-
Notifications
You must be signed in to change notification settings - Fork 1
/
convert2nwbCaImg.m
495 lines (437 loc) · 20.1 KB
/
convert2nwbCaImg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
% Convert two-photon calcium imaging data to the NWB format in Matlab
%
% Run this script to convert two-photon calcium imaging data and associated
% intracellular physiology data generated at the University of Bristol
% (UoB) to the Neurodata Without Borders (NWB) file format. This script is
% explained in the accompanying Bristol GIN for Calcium Imaging Data
% tutorial available at
% https://dervinism.github.io/bristol-neuroscience-data-guide/tutorials/Bristol%20GIN%20for%20Calcium%20Imaging%20Data.html
%
% You can use this script to get an idea of how to convert your own
% optical imaging data to the NWB file format.
%% Record metadata
% Project (experiment) metadata
projectName = 'Intracellular Ca2+ dynamics during plateau potentials trigerred by Schaffer collateral stimulation';
experimenter = 'MU';
institution = 'University of Bristol';
publications = 'In preparation';
lab = 'Jack Mellor lab';
brainArea = 'Hippocampus CA1-2';
greenIndicator = 'Fluo5f';
redIndicator = 'Alexa594';
% Animal metadata
animalID = 'm1';
ageInDays = 100;
age = ['P' num2str(ageInDays) 'D']; % Convert to ISO8601 format: https://en.wikipedia.org/wiki/ISO_8601#Durations
strain = 'C57BL/6J';
sex = 'M';
species = 'Mus musculus';
weight = [];
description = '001'; % Animal testing order.
% Session metadata
startYear = 2020;
startMonth = 12;
startDay = 4;
startTime = datetime(startYear, startMonth, startDay);
year = num2str(startYear); year = year(3:4);
month = num2str(startMonth); if numel(month) == 1; month = ['0' month]; end
day = num2str(startDay); if numel(day) == 1; day = ['0' day]; end
sliceNumber = 2;
cellNumber = 1;
imagingRate = 1/21; % A single linescan duration is 1sec with 20sec period in between linescans
lineRate = 1000.; % A number of lines scanned in a second.
sessionID = [animalID '_' year month day '_s' num2str(sliceNumber) ...
'_c' num2str(cellNumber)]; % mouse-id_time_slice-id_cell-id
sessionDescription = 'Single cell imaging in a slice combined with somatic current clamp recordings and the stimulation of Schaffer collaterals';
sessionNotes = ['201204 - Slice 2 Imaging 3 dendrite regions roughly in the SO SR and SLM regions ' ...
'Same line scan with same intensity of stim (2.3V) at different locations along the cell ' ...
'Ephys frames to match up with linescans ' ...
'Frames ' ...
'1-8 Bottom Den ' ...
'10-19 Middle Den ' ...
'22-28 Top Dendrite ' ...
'Missed a few of the imaging with the ephys so more Ephys traces than linescans ' ...
'by the end of the experiment the top or neuron started to bleb.'];
% Generate Matlab classes from NWB core schema files
generateCore;
% Assign NWB file fields
nwb = NwbFile( ...
'session_description', sessionDescription, ...
'identifier', sessionID, ...
'session_start_time', startTime, ...
'general_experimenter', experimenter, ... % optional
'general_session_id', sessionID, ... % optional
'general_institution', institution, ... % optional
'general_related_publications', publications, ... % optional
'general_notes', sessionNotes, ... % optional
'general_lab', lab); % optional
% Create subject object
subject = types.core.Subject( ...
'subject_id', animalID, ...
'age', age, ...
'description', description, ...
'species', species, ...
'sex', sex);
nwb.general_subject = subject;
%% Convert calcium imaging data
% Create optical channels
green_optical_channel = types.core.OpticalChannel( ...
'description', ['green channel corresponding to ' greenIndicator], ...
'emission_lambda', 516.);
red_optical_channel = types.core.OpticalChannel( ...
'description', ['red channel corresponding to ' redIndicator], ...
'emission_lambda', 616.);
% Create the imaging device object
device = types.core.Device(...
'description', 'Two-photon microscope', ...
'manufacturer', 'Scientifica');
nwb.general_devices.set('2P_microscope', device);
% Create imaging plane objects
imaging_plane_name = 'green_imaging_plane';
green_imaging_plane = types.core.ImagingPlane( ...
'optical_channel', green_optical_channel, ...
'description', 'The plane for imaging calcium indicator Fluo5f.', ...
'device', types.untyped.SoftLink(device), ...
'excitation_lambda', 810., ...
'imaging_rate', imagingRate, ...
'indicator', 'Fluo5f', ...
'location', brainArea);
nwb.general_optophysiology.set(imaging_plane_name, green_imaging_plane);
imaging_plane_name = 'red_imaging_plane';
red_imaging_plane = types.core.ImagingPlane( ...
'optical_channel', red_optical_channel, ...
'description', 'The plane for imaging calcium indicator Alexa594.', ...
'device', types.untyped.SoftLink(device), ...
'excitation_lambda', 810., ...
'imaging_rate', imagingRate, ...
'indicator', 'Alexa594', ...
'location', brainArea);
nwb.general_optophysiology.set(imaging_plane_name, red_imaging_plane);
% Load data
dendrite = 'Bot';
data1 = load(['..\Analysed\' year month day '__s' num2str(sliceNumber) ...
'd' num2str(cellNumber) '_004_ED__1 ' dendrite 'den_Analysed.mat']);
dendrite = 'Mid';
data2 = load(['..\Analysed\' year month day '__s' num2str(sliceNumber) ...
'd' num2str(cellNumber) '_004_ED__1 ' dendrite 'den_Analysed.mat']);
dendrite = 'Top';
data3 = load(['..\Analysed\' year month day '__s' num2str(sliceNumber) ...
'd' num2str(cellNumber) '_004_ED__1 ' dendrite 'den_Analysed.mat']);
% Append raw linescans so they have equal widths
data1.Analysed_data.Flur5_denoised = appendLinescans(data1.Analysed_data.Flur5_denoised);
data1.Analysed_data.Alexa_denoised = appendLinescans(data1.Analysed_data.Alexa_denoised);
data2.Analysed_data.Flur5_denoised = appendLinescans(data2.Analysed_data.Flur5_denoised);
data2.Analysed_data.Alexa_denoised = appendLinescans(data2.Analysed_data.Alexa_denoised);
data3.Analysed_data.Flur5_denoised = appendLinescans(data3.Analysed_data.Flur5_denoised);
data3.Analysed_data.Alexa_denoised = appendLinescans(data3.Analysed_data.Alexa_denoised);
nFrames = [size(data1.Analysed_data.Flur5_denoised,1) ...
size(data2.Analysed_data.Flur5_denoised,1) size(data3.Analysed_data.Flur5_denoised,1)];
% Add optical physiology data: Fluorescence
input.indicator = greenIndicator;
input.imagingPlane = green_imaging_plane;
input.imagingRate = imagingRate;
input.lineRate = lineRate;
input.data = data1.Analysed_data.Flur5_denoised;
input.dendriteID = 'bottom';
nwb = setTwoPhotonSeries(nwb, input);
input.data = data2.Analysed_data.Flur5_denoised;
input.dendriteID = 'middle';
nwb = setTwoPhotonSeries(nwb, input);
input.data = data3.Analysed_data.Flur5_denoised;
input.dendriteID = 'top';
nwb = setTwoPhotonSeries(nwb, input);
input.indicator = redIndicator;
input.imagingPlane = red_imaging_plane;
input.data = data1.Analysed_data.Alexa_denoised;
input.dendriteID = 'bottom';
nwb = setTwoPhotonSeries(nwb, input);
input.data = data2.Analysed_data.Alexa_denoised;
input.dendriteID = 'middle';
nwb = setTwoPhotonSeries(nwb, input);
input.data = data3.Analysed_data.Alexa_denoised;
input.dendriteID = 'top';
nwb = setTwoPhotonSeries(nwb, input);
clear input
% Add optical physiology data: Delta fluorescence
input.indicator = redIndicator;
input.imagingPlane = red_imaging_plane;
input.imagingRate = imagingRate;
input.lineRate = lineRate;
input.data = data1.Analysed_data.Calcium_deltaF';
input.dendriteID = 'bottom';
nwb = setDeltaFSeries(nwb, input);
input.data = data2.Analysed_data.Calcium_deltaF';
input.dendriteID = 'middle';
nwb = setDeltaFSeries(nwb, input);
input.data = data3.Analysed_data.Calcium_deltaF';
input.dendriteID = 'top';
nwb = setDeltaFSeries(nwb, input);
clear input
% Add images
neuron_image = types.core.RGBImage( ...
'data', data1.Analysed_data.Neuron, ... % required: [height, width, colour]
'description', 'RGB image of the full neuron.');
bottom_dend_image = types.core.GrayscaleImage( ...
'data', data1.Analysed_data.ROI_img, ... % required: [height, width]
'description', 'Grayscale image of the bottom dendrite.');
middle_dend_image = types.core.GrayscaleImage( ...
'data', data2.Analysed_data.ROI_img, ... % required: [height, width]
'description', 'Grayscale image of the middle dendrite.');
top_dend_image = types.core.GrayscaleImage( ...
'data', data3.Analysed_data.ROI_img, ... % required: [height, width]
'description', 'Grayscale image of the top dendrite.');
image_collection = types.core.Images( ...
'description', 'A collection of neuron and dendrite images.');
image_collection.image.set('neuron_image', neuron_image);
image_collection.image.set('dendrite1_image', bottom_dend_image);
image_collection.image.set('dendrite2_image', middle_dend_image);
image_collection.image.set('dendrite3_image', top_dend_image);
nwb.acquisition.set('ImageCollection', image_collection);
%% Convert intracellular electrophysiology data
% Create the recording device object
device = types.core.Device( ...
'description', 'Amplifier for recording current clamp data.', ...
'manufacturer', 'Molecular Devices');
nwb.general_devices.set('Amplifier_Multiclamp_700A', device);
electrode = types.core.IntracellularElectrode( ...
'description', 'A patch clamp electrode', ...
'location', 'Cell soma in CA1-2 of hippocampus', ...
'slice', ['slice #' num2str(sliceNumber)], ...
'device', types.untyped.SoftLink(device));
nwb.general_intracellular_ephys.set('icephys_electrode', electrode);
% Add current clamp data
input.ephysTime = data1.Analysed_data.Ephys_Time;
input.nSweeps = nFrames;
input.data = data1.Analysed_data.Ephys_data;
input.imagingRate = imagingRate;
input.electrode = electrode;
input.dendriteID = 'bottom';
nwb = setCClampSeries(nwb, input);
input.ephysTime = data2.Analysed_data.Ephys_Time;
input.data = data2.Analysed_data.Ephys_data;
input.dendriteID = 'middle';
nwb = setCClampSeries(nwb, input);
input.ephysTime = data3.Analysed_data.Ephys_Time;
input.data = data3.Analysed_data.Ephys_data;
input.dendriteID = 'top';
nwb = setCClampSeries(nwb, input);
%% Save the converted NWB file
nwbExport(nwb, [sessionID '.nwb']);
%% Local functions
function truncatedLinescans = truncateLinescans(linescans) %#ok<*DEFNU>
% truncatedLinescans = truncateLinescans(linescans)
%
% Function takes linescans of unequal widths and truncates them to have the
% smallest width and outputs the truncated linescans as a 3D matrix with
% the first dimension corresponding to the linescan number.
% Find the smallest width
nScans = numel(linescans);
length = size(linescans{1},1);
minWidth = size(linescans{1},2);
for iScan = 2:nScans
minWidth = min([minWidth size(linescans{iScan},2)]);
end
% Truncate the linescans
truncatedLinescans = zeros(nScans, length, minWidth);
for iScan = 1:nScans
diff = size(linescans{iScan},2) - minWidth;
truncatedLinescans(iScan,:,:) = linescans{iScan}(:,1+floor(diff/2):end-ceil(diff/2));
end
end
function appendedLinescans = appendLinescans(linescans)
% appendedLinescans = appendLinescans(linescans)
%
% Function takes linescans of unequal widths and appends them to have the
% the same width and outputs the appended linescans as a 3D matrix with
% the first dimension corresponding to the linescan number.
% Find the largest width
nScans = numel(linescans);
length = size(linescans{1},1);
maxWidth = size(linescans{1},2);
for iScan = 2:nScans
maxWidth = max([maxWidth size(linescans{iScan},2)]);
end
% Append the linescans
appendedLinescans = nan(nScans, length, maxWidth);
for iScan = 1:nScans
appendedLinescans(iScan,:,1:size(linescans{iScan},2)) = linescans{iScan};
end
end
function nwb = setTwoPhotonSeries(nwb, input)
% nwb = setTwoPhotonSeries(nwb, input)
%
% Function creates, names, and adds a two-photon series data to a given NWB
% file.
%
% Input: nwb - the NWB file object.
% input - a structure variable with the following fields:
% indicator - a string variable with the calcium indicator name
% (e.g., 'Fluo5f' or 'Alexa594').
% imagingPlane - the imaging plane object corresponding to the
% indicator.
% imagingRate - a scalar with the rate at which individual full
% linescans are obtained.
% lineRate - a scalar wit the frequency of individual lines within
% a single linescan.
% data - an n-dimensional matrix containing two-photon series
% linescan data. The first dimension corresponds to time
% (that is, individual linescans). The second dimension
% corresponds to individual lines along the vertical
% denditic segment. The third dimension corresponds to the
% dendritic width. NaNs are appended at the end of linescans
% to make them equal in width.
% dendriteID - a string variable with the dendritic ID (i.e.,
% 'bottom', 'middle', and 'top').
%
% Output: nwb - the NWB file object containing the newly added two-photon
% series data.
% Create image series
image_series = types.core.TwoPhotonSeries( ...
'description', [input.indicator ' linescans of the ' input.dendriteID ' dendrite'], ...
'imaging_plane', types.untyped.SoftLink(input.imagingPlane), ...
'starting_time', 0.0, ...
'starting_time_rate', input.imagingRate, ...
'scan_line_rate', input.lineRate, ...
'data', input.data, ...
'data_unit', 'a.u.', ...
'data_continuity', 'step', ...
'comments', ['This two-photon series contains ' input.indicator ' linescans of the ', ...
input.dendriteID ' (ROI) with the first dimension corresponding to time ', ...
'(or to individual linescans). Each linescan is 1-sec in duration with ', ...
'20-sec intervals between two linescans. The second dimension corresponds ', ...
'to individual lines spanning the length of the dendrite in the ROI. ', ...
'The third dimension corresponds to the width of the dendrite. ', ...
'Some linescans may contain appended NaN values to make ', ...
'widths of different linescans be equal within the same ROI.']);
% Name and add the two-photon series to the NWB file
if strcmpi(input.indicator, 'Fluo5f')
opticalChannel = 'Green';
elseif strcmpi(input.indicator, 'Alexa594')
opticalChannel = 'Red';
end
if strcmpi(input.dendriteID, 'bottom')
dendriteID = '1';
elseif strcmpi(input.dendriteID, 'middle')
dendriteID = '2';
elseif strcmpi(input.dendriteID, 'top')
dendriteID = '3';
end
nwb.acquisition.set(['TwoPhotonSeries' opticalChannel dendriteID], image_series);
end
function nwb = setDeltaFSeries(nwb, input)
% nwb = setDeltaFSeries(nwb, input)
%
% Function creates, names, and adds a two-photon series delta fluorescence
% data to a given NWB file.
%
% Input: nwb - the NWB file object.
% input - a structure variable with the following fields:
% indicator - a string variable with the calcium indicator name
% (e.g., 'Fluo5f' or 'Alexa594').
% imagingPlane - the imaging plane object corresponding to the
% indicator.
% imagingRate - a scalar with the rate at which individual full
% linescans are obtained.
% lineRate - a scalar wit the frequency of individual lines within
% a single linescan.
% data - an 2D matrix containing two-photon series delta F data.
% The first dimension corresponds to time (that is, individual
% linescans). The second dimension corresponds to individual
% lines along the vertical dendritic segment. The data is
% averaged across the dendritic width.
% dendriteID - a string variable with the dendritic ID (i.e.,
% 'bottom', 'middle', and 'top').
%
% Output: nwb - the NWB file object containing the newly added two-photon
% delta F series data.
% Create image series
image_series = types.core.TwoPhotonSeries( ...
'description', ['Delta F data for the ' input.dendriteID, ...
' calculated based on ' input.indicator '.'], ...
'imaging_plane', types.untyped.SoftLink(input.imagingPlane), ...
'starting_time', 0.0, ...
'starting_time_rate', input.imagingRate, ...
'scan_line_rate', input.lineRate, ...
'data', input.data, ...
'data_unit', 'normalised', ...
'data_continuity', 'step', ...
'comments', ['This two-photon series contains delta F data calculated based on ', ...
input.indicator ' for the ' input.dendriteID ' (ROI) ', ...
'with the first dimension corresponding to time ', ...
'(or to individual linescans). Each linescan is 1-sec in duration with ', ...
'20-sec intervals between two linescans. The second dimension corresponds ', ...
'to individual lines spanning the length of the dendrite in the ROI. ', ...
'The data is averaged across the dendritic width.']);
% Name and add the two-photon delta F series to the NWB file
if strcmpi(input.dendriteID, 'bottom')
dendriteID = '1';
elseif strcmpi(input.dendriteID, 'middle')
dendriteID = '2';
elseif strcmpi(input.dendriteID, 'top')
dendriteID = '3';
end
nwb.acquisition.set(['TwoPhotonDeltaFSeries' dendriteID], image_series);
end
function nwb = setCClampSeries(nwb, input)
% nwb = setCClampSeries(nwb, input)
%
% Function creates, names, and adds a current clamp series data to a given
% NWB file.
%
% Input: nwb - the NWB file object.
% input - a structure variable with the following fields:
% ephysTime - a time vector for a single recording sweep.
% nSweeps - The total number of sweeps for every ROI.
% imagingRate - a scalar with the rate at which individual full
% linescans are obtained.
% data - a 2D matrix containing somatic current clamp recordings
% corresponding to the initial part of the calcium
% imaging period at dendritic ROI. The first dimension
% corresponds to individual recording sweeps. The second
% dimension corresponds to individual sweep data samples.
% electrode - an electrode object.
% dendriteID - a string variable with the dendritic ID (i.e.,
% 'bottom', 'middle', and 'top').
%
% Output: nwb - the NWB file object containing the newly added current
% clamp series data.
if strcmpi(input.dendriteID, 'bottom')
nSweeps = input.nSweeps(1);
sweeps = (1:nSweeps) - 1;
dendriteID = '1';
elseif strcmpi(input.dendriteID, 'middle')
nSweeps = input.nSweeps(2);
sweeps = (input.nSweeps(1)+1:input.nSweeps(1)+nSweeps) - 1;
dendriteID = '2';
else
nSweeps = input.nSweeps(3);
sweeps = (input.nSweeps(1)+input.nSweeps(2)+1:input.nSweeps(1)+input.nSweeps(2)+nSweeps) - 1;
dendriteID = '3';
end
%timestamps = repmat((frames./input.imagingRate)', [1 numel(input.ephysTime)]) + ...
% repmat((input.ephysTime./1000)', [nSweeps 1]); % sec
timestamps = input.ephysTime'./1000; % sec
current_clamp_series = types.core.CurrentClampSeries( ...
'description', ['Somatic current clamp recording corresponding to ', ...
'the initial part of the calcium imaging period ', ...
'at ' input.dendriteID ' dendrite.'], ...
'data', input.data(:,sweeps+1)', ...
'data_continuity', 'step', ...
'data_unit', 'millivolt', ...
'starting_time', 0., ...
'starting_time_rate', input.imagingRate, ...
'electrode', types.untyped.SoftLink(input.electrode), ...
'stimulus_description', 'N/A', ...
'timestamps', timestamps, ...
'comments', ['Somatic current clamp recording corresponding to ', ...
'the initial part of the calcium imaging period ', ...
'at ' input.dendriteID ' dendrite.', ...
'The first dimension corresponds to individual ', ...
'recording sweeps. The second dimension corresponds to ', ...
'individual sweep data samples. The associated timestamps ', ...
'variable provides timestamps for the second dimension. ', ...
'This variable has to be combined with starting_time and ', ...
'starting_time_rate variables to get absolute timestamps ', ...
'for each data point']);
nwb.acquisition.set(['CurrentClampSeries' dendriteID], current_clamp_series);
end