Getting started with MatNWB and reproducing Extracellular Electrophisiology Tutorial

getting_started_matnwb.m

% Reproducing code/results of MatNWB Extracellular Electrophisiology
% Tutorial [https://www.youtube.com/watch?v=W8t4_quIl1k]

% Create NWB object and add general information
nwb = NwbFile(...
            'session_description', 'this is my first mouse experiment',...
            'identifier', 'MouseDay1',...
            'session_start_time', datetime(2021, 07, 05, 1, 24, 33),...
            'general_experimenter', 'Vigneswaran',...
            'general_institution', 'IITM',...
            'general_session_id', 'session1',...
            'general_lab', 'CNS Lab @ IITM'...
);            

% add subject information
subject = types.core.Subject(...
    'subject_id', '001', ...
    'age', 'P1D', ...
    'description', 'This is my imagination mouse', ...
    'species', 'Mus musculus', ...
    'sex', 'M');

nwb.general_subject = subject;

% add position information
position_data = [linspace(0, 1, 100), linspace(2, 3, 100)]';
spatial_series = types.core.SpatialSeries(...
                            'data', position_data, ...
                            'reference_frame', '(0, 0) is the top-left', ...
                            'timestamps', linspace(0, 3, 300));
Position = types.core.Position('SpatialSeries', spatial_series);

% create processing module
behavior_mod = types.core.ProcessingModule(...
                    'description', 'this my processed data');
behavior_mod.nwbdatainterface.set(...
                    'Position', Position);

% add the processing module to NWB file object
nwb.processing.set('behavior', behavior_mod);

% export to a nwb file
nwbExport(nwb, 'Extracell_Elect_tutorial.nwb');

% read the nwb file with new nwb file object
read_nwb = nwbRead('Extracell_Elect_tutorial.nwb');
                        
% read and print the spatial series data
r_spatial_series = read_nwb.processing.get('behavior'). ...
                   nwbdatainterface.get('Position'). ...
                   spatialseries.get('SpatialSeries');

% add trials
trials = types.core.TimeIntervals(...
            'colnames', {'start_time', 'stop_time', 'correct'}, ...
            'description', 'This contains trials and its details', ...
            'id', types.hdmf_common.ElementIdentifiers('data', 0:2), ...
            'start_time', types.hdmf_common.VectorData('data', ...
               [.1, .2, .3], 'description', 'this is start time'), ...
            'stop_time', types.hdmf_common.VectorData('data', ...
               [0.75, .175, .275], 'description','this is stop time'), ...
            'correct', types.hdmf_common.VectorData('data', ...
                [true, false, true], 'description', 'this is custom one'));
% add to nwb object
nwb.intervals_trials = trials;

% add electrode table 
n_shanks = 4;
n_channels_per_shank = 3;
% column names in the table
columns = {'x', 'y', 'z', 'imp', 'location', 'filtering', 'group', 'label'};
tbl = cell2table(cell(0, length(columns)), 'VariableNames', columns);
% device 
device = types.core.Device('description', 'this is some device', ...
                            'manufacturer', 'Imagine Probs');
device_name = 'device_x';
% add device to nwb object
nwb.general_devices.set(device_name, device);
% create a SoftLink to device
device_link = types.untyped.SoftLink(['/general/devices/' device_name]);

% populate electrode table
for ish = 1:n_shanks

    grp_name = ['shank' num2str(ish)];
    nwb.general_extracellular_ephys.set(grp_name, ...
           types.core.ElectrodeGroup( ...
           'description', ['electrode group for shank' num2str(ish)], ...
           'location', 'some brain are', ...
           'device', device_link));
    grp_object_view = types.untyped.ObjectView( ...
           ['/general/extracellular_ephys/' grp_name]);
    for ielec = 1:n_channels_per_shank
        tbl = [tbl; {1.2, 2.3, 0.004, NaN, 'Unknown', 'Unknown', ...
                     grp_object_view, [grp_name 'elec' num2str(ielec)]}];
    end
end

% add to nwb
electrode_tbl = util.table2nwb(tbl, 'all electrodes');
nwb.general_extracellular_ephys_electrodes = electrode_tbl;

% create ElectricalSeries
electrodes_object_view = types.untyped.ObjectView( ...
    '/general/electracellular_ephys/electrodes');

electrode_tbl_region = types.hdmf_common.DynamicTableRegion( ...
    'table', electrodes_object_view, ...
    'description', 'all electrodes data', ...
    'data', (0:height(tbl)-1)');

% connect with ElectricalSeries
elect_series = types.core.ElectricalSeries( ...
    'starting_time', 0.0, ...
    'starting_time_rate', 3000, ...
    'data', randn(12, 3000), ...
    'electrodes', electrode_tbl_region, ...
    'data_unit', 'volts');
% add to nwb
nwb.acquisition.set('ElectricalSeries', elect_series);

% add LFP data in similar way
elect_series2 = types.core.ElectricalSeries( ...
    'starting_time', 0.0, ...
    'starting_time_rate', 1000, ...
    'data', randn(12, 1000), ...
    'electrodes', electrode_tbl_region, ...
    'data_unit', 'volts');
lfp = types.core.LFP('ElectricalSeries', elect_series2);

ecephys_mod = types.core.ProcessingModule( ...
     'description', 'extracellular electrophysiology');
ecephys_mod.nwbdatainterface.set('LFP', lfp);
% add nwb processing module
nwb.processing.set('ecephys', ecephys_mod);

% create spike data
num_cells = 10;
fire_rate = 20;
spikes = cell(1, num_cells);
for ish = 1:num_cells
    spikes{ish} = [];
    for iu = 1:poissrnd(fire_rate)
        spikes{ish}{end+1} = cumsum(exprnd(1/fire_rate));
    end
end
spikes = cellfun(@str2double, spikes, 'UniformOutput', false);
[spike_times_vector, spike_times_index] = util.create_indexed_column( ...
    spikes, '/units/spike_times', ...
    'description', 'spikes data');

nwb.units = types.core.Units( ...
    'colnames', {'spike_times'}, ...
    'description', 'units table', ...
    'id', types.hdmf_common.ElementIdentifiers('data', int64(0:length(spikes)-1)), ...
    'spike_times', spike_times_vector, ...
    'spike_times_index', spike_times_index);

% export to nwb file
nwbExport(nwb, 'ecephys_tutorial.nwb');