diff --git a/+types/+core/Units.m b/+types/+core/Units.m index e6af0fb04..a09fa714e 100644 --- a/+types/+core/Units.m +++ b/+types/+core/Units.m @@ -131,7 +131,7 @@ obj.waveforms_index_index = p.Results.waveforms_index_index; obj.waveforms_sampling_rate = p.Results.waveforms_sampling_rate; obj.waveforms_unit = p.Results.waveforms_unit; - + % Only execute validation/setup code when called directly in this class's % constructor, not when invoked through superclass constructor chain if strcmp(class(obj), 'types.core.Units') %#ok diff --git a/.codespellrc b/.codespellrc index 5c4c3a73d..8820ec9f9 100644 --- a/.codespellrc +++ b/.codespellrc @@ -1,3 +1,3 @@ [codespell] skip = *.html,*.svg,*fastsearch.m,*.yaml,*testResults.xml -ignore-words-list = DNE,nd,whos +ignore-words-list = DNE,nd,whos,ans diff --git a/docs/source/_static/css/custom.css b/docs/source/_static/css/custom.css index 680bd1afe..153d9846d 100644 --- a/docs/source/_static/css/custom.css +++ b/docs/source/_static/css/custom.css @@ -72,3 +72,14 @@ button.copybtn { border: none; background: none; } + +.rst-content img.tutorial-media { + max-width: 100% !important; + height: auto; + display: block; + margin: 1rem auto; +} + +.rst-content .highlight-text button.copybtn { + display: none !important; +} diff --git a/docs/source/_static/html/tutorials/basicUsage.html b/docs/source/_static/html/tutorials/basicUsage.html index 478e73e4c..340b864bb 100644 --- a/docs/source/_static/html/tutorials/basicUsage.html +++ b/docs/source/_static/html/tutorials/basicUsage.html @@ -1,19 +1,21 @@ -Using NWB Data

Using NWB Data

last updated: January 1st, 2023
In this tutorial, we demonstrate the reading and usage of the NWB file produced in the File Conversion Tutorial. The output is a near-reproduction of Figure 1e from the Li et al publication, showing raster and peristimulus time histogram (PSTH) plots for neural recordings from anterior lateral motor cortex (ALM). This figure illustrates the main finding of the publication, showing the robustness of motor planning behavior and neural dynamics following short unilateral network silencing via optogenetic inhibition.

Reading NWB Files

NWB files can be read in using the nwbRead() function. This function returns a nwbfile object which is the in-memory representation of the NWB file structure.
nwb = nwbRead('out\ANM255201_20141124.nwb', 'ignorecache');

Constrained Sets

Analyzed data in NWB is placed under the analysis property, which is a Constrained Set. A constrained set consists of an arbitrary amount of key-value pairs similar to Map containers in MATLAB or a dictionary in Python. However, constrained sets also have the ability to validate their own properties closer to how a typed Object would.
You can get/set values in constrained sets using their respective .get()/.set() methods and retrieve all Set properties using the keys() method, like in a containers.Map.
unit_names = keys(nwb.analysis);

Dictionaries with MATLAB R2022b

As noted above, constrained sets uses Map container objects, a custom class that allows user to retrieve values using a corresponding key. For users working with MATLAB version R2022b, a notable addition of the dictionary object offers significantly faster functionality.
You can initialize a dictionary and configure its contents from relevant analysed NWB data with the same key/values noted above.
Number of entries in a dictionary can be done with the new numEntries() function. One can also return all the keys and values of a dictionary as arrays using keys(dict_name) and values(dict_name) functions.

Dynamic Tables

nwb.intervals_trials returns a unique type of table called a Dynamic Table. Dynamic tables inherit from the NWB type types.hdmf_common.DynamicTable and allow for a table-like interface in NWB. In the case below, we grab the special column start_time. Dynamic Tables allow adding your own vectors using the vectordata property, which are Constrained Sets. All columns are represented by either a types.hdmf_common.VectorData or a types.hdmf_common.VectorIndex type.

Data Stubs

The data property of the column id in nwb.units is a types.untyped.DataStub. This object is a representation of a dataset that is not loaded in memory, and is what allows MatNWB to lazily load its file data. To load the data into memory, use the .load() method which extracts all data from the NWB file. Alternatively, you can index into the DataStub directly using conventional MATLAB syntax.

Jagged Arrays in Dynamic Tables

With the new addition of addRow and getRow to Dynamic Tables, the concept of jagged arrays can be worked around and no longer require full understanding outside of specific data format concerns or low-level nwb tool development. The below paragraph is retained in its entirety from its original form as purely informational.
All data in a Dynamic Table must be aligned by row and column, but not all data fits into this paradigm neatly. In order to represent variable amounts of data that is localised to each row and column, NWB uses a concept called Jagged Arrays. These arrays consist of two column types: the familiar types.core.VectorData, and the new types.core.VectorIndex. A Vector Index holds no data, instead holding a reference to another Vector Data and a vector of indices that align to the Dynamic Table dimensions. The indices represent the last index boundary in the Vector Data object for the Vector Index row. As an example, an index of three in the first row of the Vector Index column points to the first three values in the referenced Vector Data column. Subsequently, if the next index were a five, it would indicate the fourth and fifth elements in the referenced Vector Data column.
The jagged arrays serve to represent multiple trials and spike times associated to each unit by id. A convenient way to represent these in MATLAB is to use Map containers where each unit's data is indexed directly by its unit id. Below, we utilize getRow in order to build the same Map.
Once again, for users with MATLAB R2022b versions or newer, we can utilize MATLAB dictionary objects instead of map containers to build the Map.
unit_ids = nwb.units.id.data.load(); % array of unit ids represented within this
 
% Initialize trials & times Map containers indexed by unit_ids
unit_trials = containers.Map('KeyType',class(unit_ids),'ValueType','any');
unit_times = containers.Map('KeyType',class(unit_ids),'ValueType','any');

MATLAB R2022b Versions

% unit_trials = dictionary();
% unit_times = dictionary();
Adding key-value pairs, one at a time:
last_idx = 0;
for i = 1:length(unit_ids)
unit_id = unit_ids(i);
row = nwb.units.getRow(unit_id, 'useId', true, 'columns', {'spike_times', 'trials'});
unit_trials(unit_id) = row.trials{1};
unit_times(unit_id) = row.spike_times{1};
end

Process Units

We now do the following for each Unit:
 
sorted_ids = sort(unit_ids);
Photostim = struct(...
'ind', true,... % mask into xs and ys for this photostim
'period', 'none',...
'duration', 0,... % in seconds
'ramp_offset', 0); % in seconds
 
% Initialize Map container of plotting data for each unit, stored as structure
Unit = containers.Map('KeyType',class(unit_ids),'ValueType','any');

MATLAB R2022b Versions

% Unit = dictionary();
Initializing unit structure:
unit_struct = struct(...
'id', [],...
'xs', [],...
'ys', [],...
'xlim', [-Inf Inf],...
'sample', 0,...
'delay', 0,...
'response', 0,...
'left_scatter', false,...
'right_scatter', false,...
'photostim', Photostim); % can have multiple photostim
for unit_id = unit_ids'
We first extract trial IDs from the Unit IDs.
unit_trial_id = unit_trials(unit_id);
Then filter out outliers from the Sample, Delay, and Response time points with which we derive a "good enough" estimate.
trial = nwb.intervals_trials.getRow(unit_trial_id, 'useId', true,...
'columns', {'PoleInTime', 'PoleOutTime', 'CueTime', 'GoodTrials'});
unit_sample = trial.PoleInTime;
unit_delay = trial.PoleOutTime;
unit_response = trial.CueTime;
unit_good_trials = trial.GoodTrials;
% Subjective parameters
delay_threshold = 0.064;
response_threshold = 0.43;
expected_delay_offset = 1.3; % determined from figure 1a
expected_response_offset = 1.3;
expected_delay = unit_sample + expected_delay_offset;
expected_response = unit_delay + expected_response_offset;
good_delay = (unit_delay > expected_delay - delay_threshold) &...
(unit_delay < expected_delay + delay_threshold);
good_response = (unit_response > expected_response - response_threshold) &...
(unit_response < expected_response + response_threshold);
avg_sample = mean(unit_sample(good_delay & good_response));
avg_delay = mean(unit_delay(good_delay & good_response));
avg_response = mean(unit_response(good_delay & good_response));
Filter the rest of the data by "good" trials.
unit_good_trials = unit_good_trials & good_delay & good_response;
unit_trial_id = unit_trial_id(unit_good_trials);
unit_spike_time = unit_times(unit_id);
unit_spike_time = unit_spike_time(unit_good_trials);
Retrieve good trial data and organize by stimulation type.
trial = nwb.intervals_trials.getRow(unit_trial_id, 'useId', true,...
'columns', {'start_time', 'HitR', 'HitL', 'StimTrials', 'PhotostimulationType'});
unit_is_photostim = logical(trial.StimTrials);
unit_stim_type = trial.PhotostimulationType;
unit_no_stim = ~unit_is_photostim & 0 == unit_stim_type;
unit_sample_stim = unit_is_photostim & 1 == unit_stim_type;
unit_early_stim = unit_is_photostim & 2 == unit_stim_type;
unit_middle_stim = unit_is_photostim & 3 == unit_stim_type;
Compose Scatter Plots and the Peristimulus Time Histogram zeroed on the Response time.
xs = unit_spike_time - trial.start_time - avg_response;
ys = unit_trial_id;
curr_unit = unit_struct;
curr_unit.xs = xs;
curr_unit.ys = ys;
curr_unit.left_scatter = logical(trial.HitL);
curr_unit.right_scatter = logical(trial.HitR);
curr_unit.sample = avg_sample - avg_response;
curr_unit.delay = avg_delay - avg_response;
curr_unit.response = 0;
% Photostim periods
curr_unit.photostim.ind = unit_no_stim;
% Sample
if any(unit_sample_stim)
SampleStim = Photostim;
SampleStim.ind = unit_sample_stim;
SampleStim.period = 'Sample';
SampleStim.duration = 0.5;
SampleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = SampleStim;
end
% Early Delay
if any(unit_early_stim)
early_stim_types = unique(unit_stim_type(unit_early_stim));
for i_early_types=1:length(early_stim_types)
early_type = early_stim_types(i_early_types);
EarlyStim = Photostim;
EarlyStim.period = 'Early Delay';
EarlyStim.ind = early_type == unit_stim_type & unit_early_stim;
if early_type == 2
EarlyStim.duration = 0.5;
EarlyStim.ramp_offset = 0.1;
else
EarlyStim.duration = 0.8;
EarlyStim.ramp_offset = 0.2;
end
curr_unit.photostim(end+1) = EarlyStim;
end
end
% Middle Delay
if any(unit_middle_stim)
MiddleStim = Photostim;
MiddleStim.ind = unit_middle_stim;
MiddleStim.period = 'Middle Delay';
MiddleStim.duration = 0.5;
MiddleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = MiddleStim;
end
Unit(unit_id) = curr_unit;
end

Plot Example Neurons

neuron_labels = [2, 3]; % neuron labels from Figure 1e
neuron_ids = [11, 2]; % neuron unit IDs corresponding to the Fig 1e labels
num_conditions = 4; % photostim conditions: nostim, sample, early, middle if applicable
num_neurons = length(neuron_ids);
 
% Initialize data structures for each summary plot of categorized neural spike data at specified stimulus condition
RasterPlot = struct(...
'xs', 0,...
'ys', 0);
ConditionPlot = struct(...
'label', '',...
'xlim', 0,...
'sample', 0,...
'delay', 0,...
'response', 0,...
'right_scatter', RasterPlot,...
'left_scatter', RasterPlot,...
'psth_bin_window', 0,...
'stim_type', '');
fig = figure;
 
% Plot neural spike data for each neuron and stimulus condition in a subplot array: num_neurons (rows) x num_conditions (columns)
for nn=1:num_neurons
Neuron = Unit(neuron_ids(nn));
% Initialize structure with neural + stimulus condition data
CurrPlot = ConditionPlot;
CurrPlot.xlim = [min(Neuron.xs) max(Neuron.xs)];
CurrPlot.sample = Neuron.sample;
CurrPlot.delay = Neuron.delay;
CurrPlot.response = Neuron.response;
% Plot each neuron/condition
plot_row = (nn - 1) * num_conditions;
for cc=1:num_conditions
ax = subplot(num_neurons, num_conditions, plot_row + cc, 'Parent', fig);
Stim = Neuron.photostim(cc);
CurrPlot.stim_type = Stim.period;
if strcmp(Stim.period, 'none')
CurrPlot.label = sprintf('Neuron %d', neuron_labels(nn));
CurrPlot.psth_bin_window = 9;
else
CurrPlot.label = Stim.period;
CurrPlot.psth_bin_window = 2;
end
stim_left_scatter_ind = Stim.ind & Neuron.left_scatter;
stim_left_scatter_trials = Neuron.ys(stim_left_scatter_ind);
CurrPlot.left_scatter.xs = Neuron.xs(stim_left_scatter_ind);
[~,CurrPlot.left_scatter.ys] = ismember(stim_left_scatter_trials,unique(stim_left_scatter_trials));
stim_right_scatter_ind = Stim.ind & Neuron.right_scatter;
stim_right_scatter_trials = Neuron.ys(stim_right_scatter_ind);
CurrPlot.right_scatter.xs = Neuron.xs(stim_right_scatter_ind);
[~,CurrPlot.right_scatter.ys] = ismember(stim_right_scatter_trials,unique(stim_right_scatter_trials));
plot_condition(ax, CurrPlot);
end
end

Helper Functions

PSTH helper function
function [psth_xs, psth_ys] = calculate_psth(xs, bin_window, bin_width)
[bin_counts, edges] = histcounts(xs, 'BinWidth', bin_width);
psth_xs = edges(1:end-1) + (bin_width / 2);
 
moving_avg_b = (1/bin_window) * ones(1,bin_window);
psth_ys = filter(moving_avg_b, 1, bin_counts);
end
Plotter function for each stimulus condition
function plot_condition(ax, ConditionPlot)
left_cdata = [1 0 0]; % red
right_cdata = [0 0 1]; % blue
hist_margin = 50;
scatter_margin = 10;
 
% Calculate PSTH values
% moving average over 200 ms as per figure 1e
hist_bin_width = 0.2 / ConditionPlot.psth_bin_window;
[left_psth_xs, left_psth_ys] =...
calculate_psth(ConditionPlot.left_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
[right_psth_xs, right_psth_ys] =...
calculate_psth(ConditionPlot.right_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
 
right_scatter_offset = min(ConditionPlot.right_scatter.ys);
right_scatter_height = max(ConditionPlot.right_scatter.ys) - right_scatter_offset;
left_scatter_offset = min(ConditionPlot.left_scatter.ys);
left_scatter_height = max(ConditionPlot.left_scatter.ys) - left_scatter_offset;
psth_height = max([left_psth_ys right_psth_ys]);
 
left_y_offset = hist_margin...
+ psth_height...
- left_scatter_offset;
 
right_y_offset = scatter_margin...
+ left_y_offset...
+ left_scatter_offset...
+ left_scatter_height...
- right_scatter_offset;
 
subplot_height = right_y_offset...
+ right_scatter_offset...
+ right_scatter_height;
 
hold(ax, 'on');
 
% PSTH
plot(ax, left_psth_xs, left_psth_ys, 'Color', left_cdata);
plot(ax, right_psth_xs, right_psth_ys, 'Color', right_cdata);
 
% Scatter Plot
scatter(ax,...
ConditionPlot.left_scatter.xs,...
left_y_offset + ConditionPlot.left_scatter.ys,...
'Marker', '.',...
'CData', left_cdata,...
'SizeData', 1);
scatter(ax,...
ConditionPlot.right_scatter.xs,...
right_y_offset + ConditionPlot.right_scatter.ys,...
'Marker', '.',...
'CData', right_cdata,...
'SizeData', 1);
 
% sample, delay, response lines
line(ax, repmat(ConditionPlot.sample, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.delay, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.response, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
 
% blue bar for photoinhibition period
if ~strcmp(ConditionPlot.stim_type, 'none')
stim_height = subplot_height;
stim_width = 0.5; % seconds
% end time relative to 'go' cue as described in the paper.
switch ConditionPlot.stim_type
case 'Sample'
end_offset = 1.6;
case 'Early Delay'
end_offset = 0.8;
case 'Middle Delay'
end_offset = 0.3;
otherwise
error('Invalid photostim period `%s`', ConditionPlot.stim_type);
end
stim_offset = ConditionPlot.response - stim_width - end_offset;
patch_vertices = [...
stim_offset, 0;...
stim_offset, stim_height;...
stim_offset+stim_width, stim_height;...
stim_offset+stim_width, 0];
patch(ax,...
'Faces', 1:4,...
'Vertices', patch_vertices,...
'FaceColor', '#B3D3EC',... % light blue shading
'EdgeColor', 'none',...
'FaceAlpha', 0.8);
end
 
title(ax, ConditionPlot.label);
xlabel(ax, 'Time (Seconds)');
ylabel(ax, 'Spikes s^{-1}')
xticks(ax, [-2 0 2]);
yticks(ax, [0 max(10, round(psth_height, -1))]);
% legend(ax, [scatter_left_plot, scatter_right_plot], {'Left Lick', 'Right Lick'},...
% 'location', 'northwestoutside');
ax.TickDir = 'out';
ax.XLim = ConditionPlot.xlim;
ax.YLim = [0 subplot_height];
 
hold(ax, 'off');
end
+.eoOutputWrapper { width: calc(90vw - 10px) !important; } +.S12 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S13 { color: rgb(33, 33, 33); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px; }

Using NWB Data

last updated: January 1st, 2023
In this tutorial, we demonstrate the reading and usage of the NWB file produced in the File Conversion Tutorial. The output is a near-reproduction of Figure 1e from the Li et al publication, showing raster and peristimulus time histogram (PSTH) plots for neural recordings from anterior lateral motor cortex (ALM). This figure illustrates the main finding of the publication, showing the robustness of motor planning behavior and neural dynamics following short unilateral network silencing via optogenetic inhibition.

Reading NWB Files

NWB files can be read in using the nwbRead() function. This function returns a nwbfile object which is the in-memory representation of the NWB file structure.
nwb = nwbRead('out\ANM255201_20141124.nwb', 'ignorecache');

Constrained Sets

Analyzed data in NWB is placed under the analysis property, which is a Constrained Set. A constrained set consists of an arbitrary amount of key-value pairs similar to Map containers in MATLAB or a dictionary in Python. However, constrained sets also have the ability to validate their own properties closer to how a typed Object would.
You can get/set values in constrained sets using their respective .get()/.set() methods and retrieve all Set properties using the keys() method, like in a containers.Map.
unit_names = keys(nwb.analysis);

Dictionaries with MATLAB R2022b

As noted above, constrained sets uses Map container objects, a custom class that allows user to retrieve values using a corresponding key. For users working with MATLAB version R2022b, a notable addition of the dictionary object offers significantly faster functionality.
You can initialize a dictionary and configure its contents from relevant analysed NWB data with the same key/values noted above.
Number of entries in a dictionary can be done with the new numEntries() function. One can also return all the keys and values of a dictionary as arrays using keys(dict_name) and values(dict_name) functions.

Dynamic Tables

nwb.intervals_trials returns a unique type of table called a Dynamic Table. Dynamic tables inherit from the NWB type types.hdmf_common.DynamicTable and allow for a table-like interface in NWB. In the case below, we grab the special column start_time. Dynamic Tables allow adding your own vectors using the vectordata property, which are Constrained Sets. All columns are represented by either a types.hdmf_common.VectorData or a types.hdmf_common.VectorIndex type.

Data Stubs

The data property of the column id in nwb.units is a types.untyped.DataStub. This object is a representation of a dataset that is not loaded in memory, and is what allows MatNWB to lazily load its file data. To load the data into memory, use the .load() method which extracts all data from the NWB file. Alternatively, you can index into the DataStub directly using conventional MATLAB syntax.

Jagged Arrays in Dynamic Tables

With the new addition of addRow and getRow to Dynamic Tables, the concept of jagged arrays can be worked around and no longer require full understanding outside of specific data format concerns or low-level nwb tool development. The below paragraph is retained in its entirety from its original form as purely informational.
All data in a Dynamic Table must be aligned by row and column, but not all data fits into this paradigm neatly. In order to represent variable amounts of data that is localised to each row and column, NWB uses a concept called Jagged Arrays. These arrays consist of two column types: the familiar types.core.VectorData, and the new types.core.VectorIndex. A Vector Index holds no data, instead holding a reference to another Vector Data and a vector of indices that align to the Dynamic Table dimensions. The indices represent the last index boundary in the Vector Data object for the Vector Index row. As an example, an index of three in the first row of the Vector Index column points to the first three values in the referenced Vector Data column. Subsequently, if the next index were a five, it would indicate the fourth and fifth elements in the referenced Vector Data column.
The jagged arrays serve to represent multiple trials and spike times associated to each unit by id. A convenient way to represent these in MATLAB is to use Map containers where each unit's data is indexed directly by its unit id. Below, we utilize getRow in order to build the same Map.
Once again, for users with MATLAB R2022b versions or newer, we can utilize MATLAB dictionary objects instead of map containers to build the Map.
unit_ids = nwb.units.id.data.load(); % array of unit ids represented within this
 
% Initialize trials & times Map containers indexed by unit_ids
unit_trials = containers.Map('KeyType',class(unit_ids),'ValueType','any');
unit_times = containers.Map('KeyType',class(unit_ids),'ValueType','any');

MATLAB R2022b Versions

% unit_trials = dictionary();
% unit_times = dictionary();
Adding key-value pairs, one at a time:
last_idx = 0;
for i = 1:length(unit_ids)
unit_id = unit_ids(i);
row = nwb.units.getRow(unit_id, 'useId', true, 'columns', {'spike_times', 'trials'});
unit_trials(unit_id) = row.trials{1};
unit_times(unit_id) = row.spike_times{1};
end

Process Units

We now do the following for each Unit:
  • Filter out invalid trials
  • Separate datasets based on resulting mouse behavior (right/left licks).
  • Derive "sample", "delay", and "response" times for this analyzed neuron.
  • Compose a peristimulus time histogram from the data.
 
sorted_ids = sort(unit_ids);
Photostim = struct(...
'ind', true,... % mask into xs and ys for this photostim
'period', 'none',...
'duration', 0,... % in seconds
'ramp_offset', 0); % in seconds
 
% Initialize Map container of plotting data for each unit, stored as structure
Unit = containers.Map('KeyType',class(unit_ids),'ValueType','any');

MATLAB R2022b Versions

% Unit = dictionary();
Initializing unit structure:
unit_struct = struct(...
'id', [],...
'xs', [],...
'ys', [],...
'xlim', [-Inf Inf],...
'sample', 0,...
'delay', 0,...
'response', 0,...
'left_scatter', false,...
'right_scatter', false,...
'photostim', Photostim); % can have multiple photostim
for unit_id = unit_ids'
We first extract trial IDs from the Unit IDs.
unit_trial_id = unit_trials(unit_id);
Then filter out outliers from the Sample, Delay, and Response time points with which we derive a "good enough" estimate.
trial = nwb.intervals_trials.getRow(unit_trial_id, 'useId', true,...
'columns', {'PoleInTime', 'PoleOutTime', 'CueTime', 'GoodTrials'});
unit_sample = trial.PoleInTime;
unit_delay = trial.PoleOutTime;
unit_response = trial.CueTime;
unit_good_trials = trial.GoodTrials;
% Subjective parameters
delay_threshold = 0.064;
response_threshold = 0.43;
expected_delay_offset = 1.3; % determined from figure 1a
expected_response_offset = 1.3;
expected_delay = unit_sample + expected_delay_offset;
expected_response = unit_delay + expected_response_offset;
good_delay = (unit_delay > expected_delay - delay_threshold) &...
(unit_delay < expected_delay + delay_threshold);
good_response = (unit_response > expected_response - response_threshold) &...
(unit_response < expected_response + response_threshold);
avg_sample = mean(unit_sample(good_delay & good_response));
avg_delay = mean(unit_delay(good_delay & good_response));
avg_response = mean(unit_response(good_delay & good_response));
Filter the rest of the data by "good" trials.
unit_good_trials = unit_good_trials & good_delay & good_response;
unit_trial_id = unit_trial_id(unit_good_trials);
unit_spike_time = unit_times(unit_id);
unit_spike_time = unit_spike_time(unit_good_trials);
Retrieve good trial data and organize by stimulation type.
trial = nwb.intervals_trials.getRow(unit_trial_id, 'useId', true,...
'columns', {'start_time', 'HitR', 'HitL', 'StimTrials', 'PhotostimulationType'});
unit_is_photostim = logical(trial.StimTrials);
unit_stim_type = trial.PhotostimulationType;
unit_no_stim = ~unit_is_photostim & 0 == unit_stim_type;
unit_sample_stim = unit_is_photostim & 1 == unit_stim_type;
unit_early_stim = unit_is_photostim & 2 == unit_stim_type;
unit_middle_stim = unit_is_photostim & 3 == unit_stim_type;
Compose Scatter Plots and the Peristimulus Time Histogram zeroed on the Response time.
xs = unit_spike_time - trial.start_time - avg_response;
ys = unit_trial_id;
curr_unit = unit_struct;
curr_unit.xs = xs;
curr_unit.ys = ys;
curr_unit.left_scatter = logical(trial.HitL);
curr_unit.right_scatter = logical(trial.HitR);
curr_unit.sample = avg_sample - avg_response;
curr_unit.delay = avg_delay - avg_response;
curr_unit.response = 0;
% Photostim periods
curr_unit.photostim.ind = unit_no_stim;
% Sample
if any(unit_sample_stim)
SampleStim = Photostim;
SampleStim.ind = unit_sample_stim;
SampleStim.period = 'Sample';
SampleStim.duration = 0.5;
SampleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = SampleStim;
end
% Early Delay
if any(unit_early_stim)
early_stim_types = unique(unit_stim_type(unit_early_stim));
for i_early_types=1:length(early_stim_types)
early_type = early_stim_types(i_early_types);
EarlyStim = Photostim;
EarlyStim.period = 'Early Delay';
EarlyStim.ind = early_type == unit_stim_type & unit_early_stim;
if early_type == 2
EarlyStim.duration = 0.5;
EarlyStim.ramp_offset = 0.1;
else
EarlyStim.duration = 0.8;
EarlyStim.ramp_offset = 0.2;
end
curr_unit.photostim(end+1) = EarlyStim;
end
end
% Middle Delay
if any(unit_middle_stim)
MiddleStim = Photostim;
MiddleStim.ind = unit_middle_stim;
MiddleStim.period = 'Middle Delay';
MiddleStim.duration = 0.5;
MiddleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = MiddleStim;
end
Unit(unit_id) = curr_unit;
end

Plot Example Neurons

neuron_labels = [2, 3]; % neuron labels from Figure 1e
neuron_ids = [11, 2]; % neuron unit IDs corresponding to the Fig 1e labels
num_conditions = 4; % photostim conditions: nostim, sample, early, middle if applicable
num_neurons = length(neuron_ids);
 
% Initialize data structures for each summary plot of categorized neural spike data at specified stimulus condition
RasterPlot = struct(...
'xs', 0,...
'ys', 0);
ConditionPlot = struct(...
'label', '',...
'xlim', 0,...
'sample', 0,...
'delay', 0,...
'response', 0,...
'right_scatter', RasterPlot,...
'left_scatter', RasterPlot,...
'psth_bin_window', 0,...
'stim_type', '');
fig = figure;
 
% Plot neural spike data for each neuron and stimulus condition in a subplot array: num_neurons (rows) x num_conditions (columns)
for nn=1:num_neurons
Neuron = Unit(neuron_ids(nn));
% Initialize structure with neural + stimulus condition data
CurrPlot = ConditionPlot;
CurrPlot.xlim = [min(Neuron.xs) max(Neuron.xs)];
CurrPlot.sample = Neuron.sample;
CurrPlot.delay = Neuron.delay;
CurrPlot.response = Neuron.response;
% Plot each neuron/condition
plot_row = (nn - 1) * num_conditions;
for cc=1:num_conditions
ax = subplot(num_neurons, num_conditions, plot_row + cc, 'Parent', fig);
Stim = Neuron.photostim(cc);
CurrPlot.stim_type = Stim.period;
if strcmp(Stim.period, 'none')
CurrPlot.label = sprintf('Neuron %d', neuron_labels(nn));
CurrPlot.psth_bin_window = 9;
else
CurrPlot.label = Stim.period;
CurrPlot.psth_bin_window = 2;
end
stim_left_scatter_ind = Stim.ind & Neuron.left_scatter;
stim_left_scatter_trials = Neuron.ys(stim_left_scatter_ind);
CurrPlot.left_scatter.xs = Neuron.xs(stim_left_scatter_ind);
[~,CurrPlot.left_scatter.ys] = ismember(stim_left_scatter_trials,unique(stim_left_scatter_trials));
stim_right_scatter_ind = Stim.ind & Neuron.right_scatter;
stim_right_scatter_trials = Neuron.ys(stim_right_scatter_ind);
CurrPlot.right_scatter.xs = Neuron.xs(stim_right_scatter_ind);
[~,CurrPlot.right_scatter.ys] = ismember(stim_right_scatter_trials,unique(stim_right_scatter_trials));
plot_condition(ax, CurrPlot);
end
end

Helper Functions

PSTH helper function
function [psth_xs, psth_ys] = calculate_psth(xs, bin_window, bin_width)
[bin_counts, edges] = histcounts(xs, 'BinWidth', bin_width);
psth_xs = edges(1:end-1) + (bin_width / 2);
 
moving_avg_b = (1/bin_window) * ones(1,bin_window);
psth_ys = filter(moving_avg_b, 1, bin_counts);
end
Plotter function for each stimulus condition
function plot_condition(ax, ConditionPlot)
left_cdata = [1 0 0]; % red
right_cdata = [0 0 1]; % blue
hist_margin = 50;
scatter_margin = 10;
 
% Calculate PSTH values
% moving average over 200 ms as per figure 1e
hist_bin_width = 0.2 / ConditionPlot.psth_bin_window;
[left_psth_xs, left_psth_ys] =...
calculate_psth(ConditionPlot.left_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
[right_psth_xs, right_psth_ys] =...
calculate_psth(ConditionPlot.right_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
 
right_scatter_offset = min(ConditionPlot.right_scatter.ys);
right_scatter_height = max(ConditionPlot.right_scatter.ys) - right_scatter_offset;
left_scatter_offset = min(ConditionPlot.left_scatter.ys);
left_scatter_height = max(ConditionPlot.left_scatter.ys) - left_scatter_offset;
psth_height = max([left_psth_ys right_psth_ys]);
 
left_y_offset = hist_margin...
+ psth_height...
- left_scatter_offset;
 
right_y_offset = scatter_margin...
+ left_y_offset...
+ left_scatter_offset...
+ left_scatter_height...
- right_scatter_offset;
 
subplot_height = right_y_offset...
+ right_scatter_offset...
+ right_scatter_height;
 
hold(ax, 'on');
 
% PSTH
plot(ax, left_psth_xs, left_psth_ys, 'Color', left_cdata);
plot(ax, right_psth_xs, right_psth_ys, 'Color', right_cdata);
 
% Scatter Plot
scatter(ax,...
ConditionPlot.left_scatter.xs,...
left_y_offset + ConditionPlot.left_scatter.ys,...
'Marker', '.',...
'CData', left_cdata,...
'SizeData', 1);
scatter(ax,...
ConditionPlot.right_scatter.xs,...
right_y_offset + ConditionPlot.right_scatter.ys,...
'Marker', '.',...
'CData', right_cdata,...
'SizeData', 1);
 
% sample, delay, response lines
line(ax, repmat(ConditionPlot.sample, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.delay, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.response, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
 
% blue bar for photoinhibition period
if ~strcmp(ConditionPlot.stim_type, 'none')
stim_height = subplot_height;
stim_width = 0.5; % seconds
% end time relative to 'go' cue as described in the paper.
switch ConditionPlot.stim_type
case 'Sample'
end_offset = 1.6;
case 'Early Delay'
end_offset = 0.8;
case 'Middle Delay'
end_offset = 0.3;
otherwise
error('Invalid photostim period `%s`', ConditionPlot.stim_type);
end
stim_offset = ConditionPlot.response - stim_width - end_offset;
patch_vertices = [...
stim_offset, 0;...
stim_offset, stim_height;...
stim_offset+stim_width, stim_height;...
stim_offset+stim_width, 0];
patch(ax,...
'Faces', 1:4,...
'Vertices', patch_vertices,...
'FaceColor', '#B3D3EC',... % light blue shading
'EdgeColor', 'none',...
'FaceAlpha', 0.8);
end
 
title(ax, ConditionPlot.label);
xlabel(ax, 'Time (Seconds)');
ylabel(ax, 'Spikes s^{-1}')
xticks(ax, [-2 0 2]);
yticks(ax, [0 max(10, round(psth_height, -1))]);
% legend(ax, [scatter_left_plot, scatter_right_plot], {'Left Lick', 'Right Lick'},...
% 'location', 'northwestoutside');
ax.TickDir = 'out';
ax.XLim = ConditionPlot.xlim;
ax.YLim = [0 subplot_height];
 
hold(ax, 'off');
end

-
- \ No newline at end of file +
+ + \ No newline at end of file diff --git a/docs/source/_static/html/tutorials/behavior.html b/docs/source/_static/html/tutorials/behavior.html index 1b1fff997..3983ee832 100644 --- a/docs/source/_static/html/tutorials/behavior.html +++ b/docs/source/_static/html/tutorials/behavior.html @@ -43,7 +43,7 @@ .S7 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif, Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: 700; text-align: left; } .S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } .S9 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 25px; min-height: 0px; white-space: pre-wrap; color: rgb(33, 33, 33); font-family: Helvetica, Arial, sans-serif, Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: 700; text-align: left; } -.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Behavior Data

This tutorial will guide you in writing behavioral data to NWB.

Creating an NWB File

Create an NWBFile object with the required fields (session_description, identifier, and session_start_time) and additional metadata.
nwb = NwbFile( ...
'session_description', 'mouse in open exploration',...
'identifier', 'Mouse5_Day3', ...
'session_start_time', datetime(2018, 4, 25, 2, 30, 3, 'TimeZone', 'local'), ...
'general_experimenter', 'My Name', ... % optional
'general_session_id', 'session_1234', ... % optional
'general_institution', 'University of My Institution', ... % optional
'general_related_publications', 'DOI:10.1016/j.neuron.2016.12.011'); % optional
nwb
nwb =
NwbFile with properties: +.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

Behavior Data

This tutorial will guide you in writing behavioral data to NWB.

Creating an NWB File

Create an NWBFile object with the required fields (session_description, identifier, and session_start_time) and additional metadata.
nwb = NwbFile( ...
'session_description', 'mouse in open exploration',...
'identifier', 'Mouse5_Day3', ...
'session_start_time', datetime(2018, 4, 25, 2, 30, 3, 'TimeZone', 'local'), ...
'general_experimenter', 'My Name', ... % optional
'general_session_id', 'session_1234', ... % optional
'general_institution', 'University of My Institution', ... % optional
'general_related_publications', 'DOI:10.1016/j.neuron.2016.12.011'); % optional
nwb
nwb =
NwbFile with properties: nwb_version: '2.9.0' file_create_date: [] @@ -96,7 +96,7 @@ stimulus_presentation: [0×1 types.untyped.Set] stimulus_templates: [0×1 types.untyped.Set] units: [] -

Subject Information

It is also recommended to store information about the experimental subject in the file. Create a Subject object to store metadata about the subject, then assign it to nwb.general_subject.
subject = types.core.Subject( ...
'subject_id', '005', ...
'age', 'P90D', ...
'description', 'mouse 5', ...
'species', 'Mus musculus', ...
'sex', 'M' ...
);
nwb.general_subject = subject;
 

SpatialSeries: Storing continuous spatial data

SpatialSeries is a subclass of TimeSeries that represents data in space, such as the spatial direction e.g., of gaze or travel or position of an animal over time.
Create data that corresponds to x, y position over time.
position_data = [linspace(0, 10, 50); linspace(0, 8, 50)]; % 2 x nT array
In SpatialSeries data, the first dimension is always time (in seconds), the second dimension represents the x, y position. However, as described in the dimensionMapNoDataPipes tutorial, when a MATLAB array is exported to HDF5, the array is transposed. Therefore, in order to correctly export the data, in MATLAB the last dimension of an array should be time. SpatialSeries data should be stored as one continuous stream as it is acquired, not by trials as is often reshaped for analysis. Data can be trial-aligned on-the-fly using the trials table. See the trials tutorial for further information.
For position data reference_frame indicates the zero-position, e.g. the 0,0 point might be the bottom-left corner of an enclosure, as viewed from the tracking camera.
timestamps = linspace(0, 50, 50)/ 200;
position_spatial_series = types.core.SpatialSeries( ...
'description', 'Postion (x, y) in an open field.', ...
'data', position_data, ...
'timestamps', timestamps, ...
'reference_frame', '(0,0) is the bottom left corner.' ...
)
position_spatial_series =
SpatialSeries with properties: +

Subject Information

It is also recommended to store information about the experimental subject in the file. Create a Subject object to store metadata about the subject, then assign it to nwb.general_subject.
subject = types.core.Subject( ...
'subject_id', '005', ...
'age', 'P90D', ...
'description', 'mouse 5', ...
'species', 'Mus musculus', ...
'sex', 'M' ...
);
nwb.general_subject = subject;
 

SpatialSeries: Storing continuous spatial data

SpatialSeries is a subclass of TimeSeries that represents data in space, such as the spatial direction e.g., of gaze or travel or position of an animal over time.
Create data that corresponds to x, y position over time.
position_data = [linspace(0, 10, 50); linspace(0, 8, 50)]; % 2 x nT array
In SpatialSeries data, the first dimension is always time (in seconds), the second dimension represents the x, y position. However, as described in the dimensionMapNoDataPipes tutorial, when a MATLAB array is exported to HDF5, the array is transposed. Therefore, in order to correctly export the data, in MATLAB the last dimension of an array should be time. SpatialSeries data should be stored as one continuous stream as it is acquired, not by trials as is often reshaped for analysis. Data can be trial-aligned on-the-fly using the trials table. See the trials tutorial for further information.
For position data reference_frame indicates the zero-position, e.g. the 0,0 point might be the bottom-left corner of an enclosure, as viewed from the tracking camera.
timestamps = linspace(0, 50, 50)/ 200;
position_spatial_series = types.core.SpatialSeries( ...
'description', 'Postion (x, y) in an open field.', ...
'data', position_data, ...
'timestamps', timestamps, ...
'reference_frame', '(0,0) is the bottom left corner.' ...
)
position_spatial_series =
SpatialSeries with properties: reference_frame: '(0,0) is the bottom left corner.' starting_time_unit: 'seconds' @@ -115,7 +115,7 @@ starting_time: [] starting_time_rate: [] timestamps: [0 0.0051 0.0102 0.0153 0.0204 0.0255 0.0306 0.0357 0.0408 0.0459 0.0510 0.0561 0.0612 0.0663 0.0714 0.0765 0.0816 0.0867 0.0918 0.0969 0.1020 0.1071 0.1122 0.1173 0.1224 0.1276 0.1327 0.1378 0.1429 0.1480 0.1531 … ] (1×50 double) -

Position: Storing position measured over time

To help data analysis and visualization tools know that this SpatialSeries object represents the position of the subject, store the SpatialSeries object inside a Position object, which can hold one or more SpatialSeries objects.
position = types.core.Position();
position.spatialseries.set('SpatialSeries', position_spatial_series);

Create a Behavior Processing Module

Create a processing module called "behavior" for storing behavioral data in the NWBFile, then add the Position object to the processing module.
behavior_processing_module = types.core.ProcessingModule('description', 'stores behavioral data.');
behavior_processing_module.nwbdatainterface.set("Position", position);
nwb.processing.set("behavior", behavior_processing_module);

CompassDirection: Storing view angle measured over time

Analogous to how position can be stored, we can create a SpatialSeries object for representing the view angle of the subject.
For direction data reference_frame indicates the zero direction, for instance in this case "straight ahead" is 0 radians.
view_angle_data = linspace(0, 4, 50);
direction_spatial_series = types.core.SpatialSeries( ...
'description', 'View angle of the subject measured in radians.', ...
'data', view_angle_data, ...
'timestamps', timestamps, ...
'reference_frame', 'straight ahead', ...
'data_unit', 'radians' ...
);
direction = types.core.CompassDirection();
direction.spatialseries.set('spatial_series', direction_spatial_series);
We can add a CompassDirection object to the behavior processing module the same way we have added the position data.
behavior_processing_module.nwbdatainterface.set('CompassDirection', direction);

BehaviorTimeSeries: Storing continuous behavior data

BehavioralTimeSeries is an interface for storing continuous behavior data, such as the speed of a subject.
speed_data = linspace(0, 0.4, 50);
 
speed_time_series = types.core.TimeSeries( ...
'data', speed_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 10.0, ... % Hz
'description', 'he speed of the subject measured over time.', ...
'data_unit', 'm/s' ...
);
 
behavioral_time_series = types.core.BehavioralTimeSeries();
behavioral_time_series.timeseries.set('speed', speed_time_series);
 
% Add behavioral_time_series to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralTimeSeries', behavioral_time_series);

BehavioralEvents: Storing behavioral events

BehavioralEvents is an interface for storing behavioral events. We can use it for storing the timing and amount of rewards (e.g. water amount) or lever press times.
reward_amount = [1.0, 1.5, 1.0, 1.5];
event_timestamps = [1.0, 2.0, 5.0, 6.0];
 
time_series = types.core.TimeSeries( ...
'data', reward_amount, ...
'timestamps', event_timestamps, ...
'description', 'The water amount the subject received as a reward.', ...
'data_unit', 'ml' ...
);
 
behavioral_events = types.core.BehavioralEvents();
behavioral_events.timeseries.set('lever_presses', time_series);
 
% Add behavioral_events to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralEvents', behavioral_events);
Storing only the timestamps of the events is possible with the ndx-events NWB extension. You can also add labels associated with the events with this extension. You can find information about installation and example usage here.

BehavioralEpochs: Storing intervals of behavior data

BehavioralEpochs is for storing intervals of behavior data. BehavioralEpochs uses IntervalSeries to represent the time intervals. Create an IntervalSeries object that represents the time intervals when the animal was running. IntervalSeries uses 1 to indicate the beginning of an interval and -1 to indicate the end.
run_intervals = types.core.IntervalSeries( ...
'description', 'Intervals when the animal was running.', ...
'data', [1, -1, 1, -1, 1, -1], ...
'timestamps', [0.5, 1.5, 3.5, 4.0, 7.0, 7.3] ...
);
 
behavioral_epochs = types.core.BehavioralEpochs();
behavioral_epochs.intervalseries.set('running', run_intervals);
You can add more than one IntervalSeries to a BehavioralEpochs object.
sleep_intervals = types.core.IntervalSeries( ...
'description', 'Intervals when the animal was sleeping', ...
'data', [1, -1, 1, -1], ...
'timestamps', [15.0, 30.0, 60.0, 95.0] ...
);
behavioral_epochs.intervalseries.set('sleeping', sleep_intervals);
 
% Add behavioral_epochs to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralEpochs', behavioral_epochs);

Another approach: TimeIntervals

Using TimeIntervals to represent time intervals is often preferred over BehavioralEpochs and IntervalSeries. TimeIntervals is a subclass of DynamicTable, which offers flexibility for tabular data by allowing the addition of optional columns which are not defined in the standard DynamicTable class.
sleep_intervals = types.core.TimeIntervals( ...
'description', 'Intervals when the animal was sleeping.', ...
'colnames', {'start_time', 'stop_time', 'stage'} ...
);
 
sleep_intervals.addRow('start_time', 0.3, 'stop_time', 0.35, 'stage', 1);
sleep_intervals.addRow('start_time', 0.7, 'stop_time', 0.9, 'stage', 2);
sleep_intervals.addRow('start_time', 1.3, 'stop_time', 3.0, 'stage', 3);
 
nwb.intervals.set('sleep_intervals', sleep_intervals);

EyeTracking: Storing continuous eye-tracking data of gaze direction

EyeTracking is for storing eye-tracking data which represents direction of gaze as measured by an eye tracking algorithm. An EyeTracking object holds one or more SpatialSeries objects that represent the gaze direction over time extracted from a video.
eye_position_data = [linspace(-20, 30, 50); linspace(30, -20, 50)];
 
right_eye_position = types.core.SpatialSeries( ...
'description', 'The position of the right eye measured in degrees.', ...
'data', eye_position_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 50.0, ... % Hz
'reference_frame', '(0,0) is middle', ...
'data_unit', 'degrees' ...
);
 
left_eye_position = types.core.SpatialSeries( ...
'description', 'The position of the right eye measured in degrees.', ...
'data', eye_position_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 50.0, ... % Hz
'reference_frame', '(0,0) is middle', ...
'data_unit', 'degrees' ...
);
 
eye_tracking = types.core.EyeTracking();
eye_tracking.spatialseries.set('right_eye_position', right_eye_position);
eye_tracking.spatialseries.set('left_eye_position', left_eye_position);
 
behavior_processing_module.nwbdatainterface.set('EyeTracking', eye_tracking);

PupilTracking: Storing continuous eye-tracking data of pupil size

PupilTracking is for storing eye-tracking data which represents pupil size. PupilTracking holds one or more TimeSeries objects that can represent different features such as the dilation of the pupil measured over time by a pupil tracking algorithm.
pupil_diameter = types.core.TimeSeries( ...
'description', 'Pupil diameter extracted from the video of the right eye.', ...
'data', linspace(0.001, 0.002, 50), ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 20.0, ... % Hz
'data_unit', 'meters' ...
);
 
pupil_tracking = types.core.PupilTracking();
pupil_tracking.timeseries.set('pupil_diameter', pupil_diameter);
 
behavior_processing_module.nwbdatainterface.set('PupilTracking', pupil_tracking);

Writing the behavior data to an NWB file

All of the above commands build an NWBFile object in-memory. To write this file, use nwbExport.
% Save to tutorials/tutorial_nwb_files folder
nwbFilePath = misc.getTutorialNwbFilePath('behavior_tutorial.nwb');
nwbExport(nwb, nwbFilePath);
fprintf('Exported NWB file to "%s"\n', 'behavior_tutorial.nwb')
Exported NWB file to "behavior_tutorial.nwb"
+

Position: Storing position measured over time

To help data analysis and visualization tools know that this SpatialSeries object represents the position of the subject, store the SpatialSeries object inside a Position object, which can hold one or more SpatialSeries objects.
position = types.core.Position();
position.spatialseries.set('SpatialSeries', position_spatial_series);

Create a Behavior Processing Module

Create a processing module called "behavior" for storing behavioral data in the NWBFile, then add the Position object to the processing module.
behavior_processing_module = types.core.ProcessingModule('description', 'stores behavioral data.');
behavior_processing_module.nwbdatainterface.set("Position", position);
nwb.processing.set("behavior", behavior_processing_module);

CompassDirection: Storing view angle measured over time

Analogous to how position can be stored, we can create a SpatialSeries object for representing the view angle of the subject.
For direction data reference_frame indicates the zero direction, for instance in this case "straight ahead" is 0 radians.
view_angle_data = linspace(0, 4, 50);
direction_spatial_series = types.core.SpatialSeries( ...
'description', 'View angle of the subject measured in radians.', ...
'data', view_angle_data, ...
'timestamps', timestamps, ...
'reference_frame', 'straight ahead', ...
'data_unit', 'radians' ...
);
direction = types.core.CompassDirection();
direction.spatialseries.set('spatial_series', direction_spatial_series);
We can add a CompassDirection object to the behavior processing module the same way we have added the position data.
behavior_processing_module.nwbdatainterface.set('CompassDirection', direction);

BehaviorTimeSeries: Storing continuous behavior data

BehavioralTimeSeries is an interface for storing continuous behavior data, such as the speed of a subject.
speed_data = linspace(0, 0.4, 50);
 
speed_time_series = types.core.TimeSeries( ...
'data', speed_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 10.0, ... % Hz
'description', 'he speed of the subject measured over time.', ...
'data_unit', 'm/s' ...
);
 
behavioral_time_series = types.core.BehavioralTimeSeries();
behavioral_time_series.timeseries.set('speed', speed_time_series);
 
% Add behavioral_time_series to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralTimeSeries', behavioral_time_series);

BehavioralEvents: Storing behavioral events

BehavioralEvents is an interface for storing behavioral events. We can use it for storing the timing and amount of rewards (e.g. water amount) or lever press times.
reward_amount = [1.0, 1.5, 1.0, 1.5];
event_timestamps = [1.0, 2.0, 5.0, 6.0];
 
time_series = types.core.TimeSeries( ...
'data', reward_amount, ...
'timestamps', event_timestamps, ...
'description', 'The water amount the subject received as a reward.', ...
'data_unit', 'ml' ...
);
 
behavioral_events = types.core.BehavioralEvents();
behavioral_events.timeseries.set('lever_presses', time_series);
 
% Add behavioral_events to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralEvents', behavioral_events);
Storing only the timestamps of the events is possible with the ndx-events NWB extension. You can also add labels associated with the events with this extension. You can find information about installation and example usage here.

BehavioralEpochs: Storing intervals of behavior data

BehavioralEpochs is for storing intervals of behavior data. BehavioralEpochs uses IntervalSeries to represent the time intervals. Create an IntervalSeries object that represents the time intervals when the animal was running. IntervalSeries uses 1 to indicate the beginning of an interval and -1 to indicate the end.
run_intervals = types.core.IntervalSeries( ...
'description', 'Intervals when the animal was running.', ...
'data', [1, -1, 1, -1, 1, -1], ...
'timestamps', [0.5, 1.5, 3.5, 4.0, 7.0, 7.3] ...
);
 
behavioral_epochs = types.core.BehavioralEpochs();
behavioral_epochs.intervalseries.set('running', run_intervals);
You can add more than one IntervalSeries to a BehavioralEpochs object.
sleep_intervals = types.core.IntervalSeries( ...
'description', 'Intervals when the animal was sleeping', ...
'data', [1, -1, 1, -1], ...
'timestamps', [15.0, 30.0, 60.0, 95.0] ...
);
behavioral_epochs.intervalseries.set('sleeping', sleep_intervals);
 
% Add behavioral_epochs to the processing module
behavior_processing_module.nwbdatainterface.set('BehavioralEpochs', behavioral_epochs);

Another approach: TimeIntervals

Using TimeIntervals to represent time intervals is often preferred over BehavioralEpochs and IntervalSeries. TimeIntervals is a subclass of DynamicTable, which offers flexibility for tabular data by allowing the addition of optional columns which are not defined in the standard DynamicTable class.
sleep_intervals = types.core.TimeIntervals( ...
'description', 'Intervals when the animal was sleeping.', ...
'colnames', {'start_time', 'stop_time', 'stage'} ...
);
 
sleep_intervals.addRow('start_time', 0.3, 'stop_time', 0.35, 'stage', 1);
sleep_intervals.addRow('start_time', 0.7, 'stop_time', 0.9, 'stage', 2);
sleep_intervals.addRow('start_time', 1.3, 'stop_time', 3.0, 'stage', 3);
 
nwb.intervals.set('sleep_intervals', sleep_intervals);

EyeTracking: Storing continuous eye-tracking data of gaze direction

EyeTracking is for storing eye-tracking data which represents direction of gaze as measured by an eye tracking algorithm. An EyeTracking object holds one or more SpatialSeries objects that represent the gaze direction over time extracted from a video.
eye_position_data = [linspace(-20, 30, 50); linspace(30, -20, 50)];
 
right_eye_position = types.core.SpatialSeries( ...
'description', 'The position of the right eye measured in degrees.', ...
'data', eye_position_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 50.0, ... % Hz
'reference_frame', '(0,0) is middle', ...
'data_unit', 'degrees' ...
);
 
left_eye_position = types.core.SpatialSeries( ...
'description', 'The position of the right eye measured in degrees.', ...
'data', eye_position_data, ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 50.0, ... % Hz
'reference_frame', '(0,0) is middle', ...
'data_unit', 'degrees' ...
);
 
eye_tracking = types.core.EyeTracking();
eye_tracking.spatialseries.set('right_eye_position', right_eye_position);
eye_tracking.spatialseries.set('left_eye_position', left_eye_position);
 
behavior_processing_module.nwbdatainterface.set('EyeTracking', eye_tracking);

PupilTracking: Storing continuous eye-tracking data of pupil size

PupilTracking is for storing eye-tracking data which represents pupil size. PupilTracking holds one or more TimeSeries objects that can represent different features such as the dilation of the pupil measured over time by a pupil tracking algorithm.
pupil_diameter = types.core.TimeSeries( ...
'description', 'Pupil diameter extracted from the video of the right eye.', ...
'data', linspace(0.001, 0.002, 50), ...
'starting_time', 1.0, ... % NB: Important to set starting_time when using starting_time_rate
'starting_time_rate', 20.0, ... % Hz
'data_unit', 'meters' ...
);
 
pupil_tracking = types.core.PupilTracking();
pupil_tracking.timeseries.set('pupil_diameter', pupil_diameter);
 
behavior_processing_module.nwbdatainterface.set('PupilTracking', pupil_tracking);

Writing the behavior data to an NWB file

All of the above commands build an NWBFile object in-memory. To write this file, use nwbExport.
% Save to tutorials/tutorial_nwb_files folder
nwbFilePath = misc.getTutorialNwbFilePath('behavior_tutorial.nwb');
nwbExport(nwb, nwbFilePath);
fprintf('Exported NWB file to "%s"\n', 'behavior_tutorial.nwb')
Exported NWB file to "behavior_tutorial.nwb"

diff --git a/docs/source/_static/html/tutorials/read_demo.html b/docs/source/_static/html/tutorials/read_demo.html index 2f1c4e6d0..9be65ca66 100644 --- a/docs/source/_static/html/tutorials/read_demo.html +++ b/docs/source/_static/html/tutorials/read_demo.html @@ -88,7 +88,7 @@ Stimulus Quick PSTH and raster Plot results -Conclusion

Download the data

First, let's download an NWB data file from the DANDI neurophysiology data archive.
A NWB file represents a single session of an experiment. It contains all the data of that session and the metadata required to understand the data.
We will use data from one session of an experiment by Chandravadia et al. (2020), where the authors recorded single neuron electrophysiological activity from the medial temporal lobes of human subjects while they performed a visual recognition memory task.
  1. Go to the DANDI page for this dataset: https://dandiarchive.org/dandiset/000004/draft
  2. Toward the top middle of the page, click the "Files" button.
demo_dandi_view_files_in_dataset.png
3. Click on the folder "sub-P11MHM" (click the folder name, not the checkbox).
demo_dandi_select_folder.png
4. Then click on the download symbol to the right of the filename "sub-P11HMH_ses-20061101_ecephys+image.nwb" to download the data file (69 MB) to your computer.
demo_dandi_download_data.png

Installing matnwb

Use the code below to install MatNWB from source using git. Ensure git is on your path before running this line.
if ~exist('nwbRead', 'file') % Skip if MatNWB is on MATLAB's search path
!git clone https://github.com/NeurodataWithoutBorders/matnwb.git
% add the path to matnwb and generate the core classes
addpath('matnwb');
end
MatNWB works by automatically creating API classes based on the schema. For most NWB files, the classes are generated automatically by calling nwbRead farther down. This particular NWB file was created before this feature was supported, so we must ensure that these classes for the correct schema versions are properly generated before attempting to read from the file.
% Reminder: YOU DO NOT NORMALLY NEED TO CALL THIS FUNCTION. Only attempt this
% method if you encounter read errors.
generateCore(util.getSchemaVersion('sub-P11HMH_ses-20061101_ecephys+image.nwb'));

Read the NWB file

You can read any NWB file using nwbRead. You will find that the print out for this shows a summary of the data within.
% ignorecache informs the `nwbRead` call to not generate files by default. Since
% we have already done this, we can skip this automated step when reading. If
% you are reading the file before generating, you can omit this argument flag.
nwb = nwbRead('sub-P11HMH_ses-20061101_ecephys+image.nwb', 'ignorecache')
nwb =
NwbFile with properties: +Conclusion

Download the data

First, let's download an NWB data file from the DANDI neurophysiology data archive.
A NWB file represents a single session of an experiment. It contains all the data of that session and the metadata required to understand the data.
We will use data from one session of an experiment by Chandravadia et al. (2020), where the authors recorded single neuron electrophysiological activity from the medial temporal lobes of human subjects while they performed a visual recognition memory task.
  1. Go to the DANDI page for this dataset: https://dandiarchive.org/dandiset/000004/draft
  2. Toward the top middle of the page, click the "Files" button.
demo_dandi_view_files_in_dataset.png
3. Click on the folder "sub-P11MHM" (click the folder name, not the checkbox).
demo_dandi_select_folder.png
4. Then click on the download symbol to the right of the filename "sub-P11HMH_ses-20061101_ecephys+image.nwb" to download the data file (69 MB) to your computer.
demo_dandi_download_data.png

Installing matnwb

Use the code below to install MatNWB from source using git. Ensure git is on your path before running this line.
if ~exist('nwbRead', 'file') % Skip if MatNWB is on MATLAB's search path
!git clone https://github.com/NeurodataWithoutBorders/matnwb.git
% add the path to matnwb and generate the core classes
addpath('matnwb');
end
MatNWB works by automatically creating API classes based on the schema. For most NWB files, the classes are generated automatically by calling nwbRead farther down. This particular NWB file was created before this feature was supported, so we must ensure that these classes for the correct schema versions are properly generated before attempting to read from the file.
% Reminder: YOU DO NOT NORMALLY NEED TO CALL THIS FUNCTION. Only attempt this
% method if you encounter read errors.
generateCore(util.getSchemaVersion('sub-P11HMH_ses-20061101_ecephys+image.nwb'));

Read the NWB file

You can read any NWB file using nwbRead. You will find that the print out for this shows a summary of the data within.
% ignorecache informs the `nwbRead` call to not generate files by default. Since
% we have already done this, we can skip this automated step when reading. If
% you are reading the file before generating, you can omit this argument flag.
nwb = nwbRead('sub-P11HMH_ses-20061101_ecephys+image.nwb', 'ignorecache')
nwb =
NwbFile with properties: nwb_version: '2.1.0' file_create_date: [1×1 types.untyped.DataStub] @@ -134,10 +134,10 @@ stimulus_presentation: [1×1 types.untyped.Set] stimulus_templates: [0×1 types.untyped.Set] units: [1×1 types.core.Units] -
You can also use util.nwbTree to actively explore the NWB file.
util.nwbTree(nwb);

Stimulus

Now lets take a look at the visual stimuli presented to the subject. They will be in nwb.stimulus_presentation
nwb.stimulus_presentation
ans =
Set with properties: +
You can also use util.nwbTree to actively explore the NWB file.
util.nwbTree(nwb);

Stimulus

Now lets take a look at the visual stimuli presented to the subject. They will be in nwb.stimulus_presentation
nwb.stimulus_presentation
ans =
Set with properties: StimulusPresentation: [types.core.OpticalSeries] -
This results shows us that nwb.stimulus_presentation is a Set object that contains a single data object called StimulusPresentation, which is an OpticalSeries neurodata type. Use the get method to return this OpticalSeries. Set objects store a collection of other NWB objects.
nwb.stimulus_presentation.get('StimulusPresentation')
ans =
OpticalSeries with properties: +
This results shows us that nwb.stimulus_presentation is a Set object that contains a single data object called StimulusPresentation, which is an OpticalSeries neurodata type. Use the get method to return this OpticalSeries. Set objects store a collection of other NWB objects.
nwb.stimulus_presentation.get('StimulusPresentation')
ans =
OpticalSeries with properties: distance: 0.7000 field_of_view: [1×1 types.untyped.DataStub] @@ -160,15 +160,15 @@ starting_time: [] starting_time_rate: [] timestamps: [1×1 types.untyped.DataStub] -
OpticalSeries is a neurodata type that stores information about visual stimuli presented to subjects. This print out shows all of the attributes in the OpticalSeries object named StimulusPresentation. The images are stored in StimulusPresentation.data
StimulusImageData = nwb.stimulus_presentation.get('StimulusPresentation').data
StimulusImageData =
DataStub with properties: +
OpticalSeries is a neurodata type that stores information about visual stimuli presented to subjects. This print out shows all of the attributes in the OpticalSeries object named StimulusPresentation. The images are stored in StimulusPresentation.data
StimulusImageData = nwb.stimulus_presentation.get('StimulusPresentation').data
StimulusImageData =
DataStub with properties: filename: 'sub-P11HMH_ses-20061101_ecephys+image.nwb' path: '/stimulus/presentation/StimulusPresentation/data' dims: [3 300 400 200] ndims: 4 dataType: 'uint8' -
When calling a data object directly, the data is not read but instead a DataStub is returned. This is because data is read "lazily" in MatNWB. Instead of reading the entire dataset into memory, this provides a "window" into the data stored on disk that allows you to read only a section of the data. In this case, the last dimension indexes over images. You can index into any DataStub as you would any MATLAB matrix.
% get the image and display it
% the dimension order is provided as follows:
% [rgb, y, x, image index]
img = StimulusImageData(1:3, 1:300, 1:400, 32);
A bit of manipulation allows us to display the image using MATLAB's imshow.
img = permute(img,[3, 2, 1]); % fix orientation
img = flip(img, 3); % reverse color order
F = figure();
imshow(img, 'InitialMagnification', 'fit');
daspect([3, 5, 5]);
To read an entire dataset, use the DataStub.load method without any input arguments. We will use this approach to read all of the image display timestamps into memory.
stimulus_times = nwb.stimulus_presentation.get('StimulusPresentation').timestamps.load();

Quick PSTH and raster

Here, I will pull out spike times of a particular unit, align them to the image display times, and finally display the results.
First, let us show the first row of the NWB Units table representing the first unit.
nwb.units.getRow(1)
ans = 1×8 table
 origClusterIDwaveform_mean_encodingwaveform_mean_recognitionIsolationDistSNRwaveform_mean_sampling_ratespike_timeselectrodes
11102256×1 double256×1 double11.29171.440798400373×1 double0
Let us specify some parameters for creating a cell array of spike times aligned to each stimulus time.
%% Align spikes by stimulus presentations
 
unit_ind =8;
before =1;
after =3;
getRow provides a convenient method for reading this data out.
unit_spikes = nwb.units.getRow(unit_ind, 'columns', {'spike_times'}).spike_times{1}
unit_spikes = 2116×1
103 ×
5.9338 +
When calling a data object directly, the data is not read but instead a DataStub is returned. This is because data is read "lazily" in MatNWB. Instead of reading the entire dataset into memory, this provides a "window" into the data stored on disk that allows you to read only a section of the data. In this case, the last dimension indexes over images. You can index into any DataStub as you would any MATLAB matrix.
% get the image and display it
% the dimension order is provided as follows:
% [rgb, y, x, image index]
img = StimulusImageData(1:3, 1:300, 1:400, 32);
A bit of manipulation allows us to display the image using MATLAB's imshow.
img = permute(img,[3, 2, 1]); % fix orientation
img = flip(img, 3); % reverse color order
F = figure();
imshow(img, 'InitialMagnification', 'fit');
daspect([3, 5, 5]);
To read an entire dataset, use the DataStub.load method without any input arguments. We will use this approach to read all of the image display timestamps into memory.
stimulus_times = nwb.stimulus_presentation.get('StimulusPresentation').timestamps.load();

Quick PSTH and raster

Here, I will pull out spike times of a particular unit, align them to the image display times, and finally display the results.
First, let us show the first row of the NWB Units table representing the first unit.
nwb.units.getRow(1)
ans = 1×8 table
 origClusterIDwaveform_mean_encodingwaveform_mean_recognitionIsolationDistSNRwaveform_mean_sampling_ratespike_timeselectrodes
11102256×1 double256×1 double11.29171.440798400373×1 double0
Let us specify some parameters for creating a cell array of spike times aligned to each stimulus time.
%% Align spikes by stimulus presentations
 
unit_ind =8;
before =1;
after =3;
getRow provides a convenient method for reading this data out.
unit_spikes = nwb.units.getRow(unit_ind, 'columns', {'spike_times'}).spike_times{1}
unit_spikes = 2116×1
103 ×
5.9338 5.9343 5.9346 5.9358 @@ -187,7 +187,7 @@ 5.9358 5.9364 5.9375 - " style="white-space: normal; font-style: normal; color: rgb(33, 33, 33); font-size: 12px;">
Spike times from this unit are aligned to each stimulus time and compiled in a cell array
results = cell(1, length(stimulus_times));
for itime = 1:length(stimulus_times)
stimulus_time = stimulus_times(itime);
spikes = unit_spikes - stimulus_time;
spikes = spikes(spikes > -before);
spikes = spikes(spikes < after);
results{itime} = spikes;
end

Plot results

Finally, here is a (slightly sloppy) peri-stimulus time histogram
figure();
hold on
for i = 1:length(results)
spikes = results{i};
yy = ones(length(spikes)) * i;
 
plot(spikes, yy, 'k.');
end
hold off
ylabel('trial');
xlabel('time (s)');
axis('tight')
figure();
all_spikes = cat(1, results{:});
histogram(all_spikes, 30);
ylabel('count')
xlabel('time (s)');
axis('tight')

Conclusion

This is an example of how to get started with understanding and analyzing public NWB datasets. This particular dataset was published with an extensive open analysis conducted in both MATLAB and Python, which you can find here. For more datasets, or to publish your own NWB data for free, check out the DANDI archive here. Also, make sure to check out the DANDI breakout session later in this event.
+ " style="white-space: normal; font-style: normal; color: rgb(33, 33, 33); font-size: 12px;">
Spike times from this unit are aligned to each stimulus time and compiled in a cell array
results = cell(1, length(stimulus_times));
for itime = 1:length(stimulus_times)
stimulus_time = stimulus_times(itime);
spikes = unit_spikes - stimulus_time;
spikes = spikes(spikes > -before);
spikes = spikes(spikes < after);
results{itime} = spikes;
end

Plot results

Finally, here is a (slightly sloppy) peri-stimulus time histogram
figure();
hold on
for i = 1:length(results)
spikes = results{i};
yy = ones(length(spikes)) * i;
 
plot(spikes, yy, 'k.');
end
hold off
ylabel('trial');
xlabel('time (s)');
axis('tight')
figure();
all_spikes = cat(1, results{:});
histogram(all_spikes, 30);
ylabel('count')
xlabel('time (s)');
axis('tight')

Conclusion

This is an example of how to get started with understanding and analyzing public NWB datasets. This particular dataset was published with an extensive open analysis conducted in both MATLAB and Python, which you can find here. For more datasets, or to publish your own NWB data for free, check out the DANDI archive here. Also, make sure to check out the DANDI breakout session later in this event.