catalystneuro / ndx-bipolar-scheme

Structure for storing the bipolar schema of a recording in an NWB file.
BSD 3-Clause "New" or "Revised" License
1 stars 2 forks source link

without an electrical series #18

Closed nmarkowitz closed 3 years ago

nmarkowitz commented 3 years ago

Hi,

I'm trying to store information on different reference schemes on my nwb file for intracranial recordings. However I don't want to create a new electricalseries to hold the rereferenced signal so the file doesn't get extremely large. What I'd like to store information for is:

Other possible references will probably be used. For now, for each of these, I'm storing a DynamicTable in a ProcessingModule that lists which electrodes compose it, what to call the new trace and which acquisition the signal is derived from using a soft link. This last one is important because I'm not storing an electricalseries object of all the re-referenced traces in order to preserve space. Below is an example of how I do this in matnwb:

% The processing module
bipolar_proc= types.core.ProcessingModule('description', 'channel pairs to apply bipolar reference to');

% Data for table
label = {'Elec1-Elec2','Elec2-Elec3','Elec3-Elec4'}; % Names of the channels. The trace is first channel minus second channel

% Parse the name to get the two channels that compose it
nchans = length(label);
label_a = split(label,'-');
label_b = label_a(:,2);
label_a = label_a(:,1);

% Some soft links
lfp_link = types.untyped.ObjectView('/acquisition/ieeg'); % Where the raw traces of the signal is located so it can be reproduced
lfp_link_col = repmat(lfp_link, nchans, 1);
electable_link = types.untyped.ObjectView('/general/extracellular_ephys/electrodes');
electable_link_col = repmat(electable_link, nchans, 1);

% Construct DynamicTable
tbl = table(label, label_a, label_b, lfp_link_col, electable_link_col, 'VariableNames', {'label','a', 'b','lfp','electrodes'});
dyn_table = util.table2nwb(tbl);
dyn_table.description = 'channels to pair for bipolar referencing';
dyn_table.vectordata.get('label').description = 'name of bipolar pair';
dyn_table.vectordata.get('a').description = 'first channel of bipolar pair';
dyn_table.vectordata.get('b').description = 'second channel of bipolar pair';
dyn_table.vectordata.get('lfp').description = 'acquisition in which channels are found';
dyn_table.vectordata.get('electrodes').description = 'original electrode table containing all information for each channel';

% Add table to processing module and set in NWB file object
bipolar_proc.dynamictable.set('channels',dyn_table);
ecog.nwb.processing.set('bipolar_ref', bipolar_proc);

With all this in mind what would you recommend I do? Is this extension not intended for these purposes or do I just need to add more information to the created table? Please let me know when you get a chance and thank you in advance.

bendichter commented 3 years ago

@sportnoah14 thanks for this demo code and the full description of what you are trying to do.

With this extension, you can indeed store the BipolarSchema in the /general group without having to store the corresponding ElectricalSeries, so I think this extension could work for you, though it was designed specifically with bipolar schemes in mind, where there are 1+ anode(s) and 1+ cathode(s), so I'm not sure how well it would work for the other two modes. I suppose with avg signal you could store anodes and have cathodes be an empty list for each row. I'm not exactly sure how the white matter part would work, but maybe something similar.

Right now you are storing electrode info in 2 columns, the first indexes the electrodes table and the second just repeats a softlink to the electrodes table you want to target for every row of the table. The more NWB-ish way to do this, and the way it is done here, is to use a DynamicTableRegion, which stores indexes just like you have here, but instead of storing the link to the electrodes table as another column, stores it as a reference to the target table as an attribute of the dataset.

Are you able to create MATLAB code for this? If you are interested in using it we could add some MATLAB example code to the README

nmarkowitz commented 3 years ago

@bendichter , this is what I came up with to integrate what I'd like to accomplish with this extension. To fully reconstruct my referencing I have a column that specifies which acquisition the anodes/cathodes came from and their channel indices within that acquisition (ex: average reference channels are from /acquisition/ieeg channels [0,1,2,5,9]). This would work like a reverse DynamicTableRegion. Based on this code what do you suggest? And if you think it looks good then we can definitely add it as example matlab code

% soft links to use
electable_link = types.untyped.ObjectView('/general/extracellular_ephys/electrodes');
lfp_link = types.untyped.ObjectView('/acquisition/ieeg');

%% Average reference
elec_labels = {'LDa1','LDa2','LDa3','LDa4'};
elec_ids = [0:3]; % ID of the electrode in the primary electrode table
avg_ref_elecs = {'LDa1','LDa2','LDa4'}; % The electrodes that make up the average reference
avg_ref_ids = elec_ids(ismember(elec_labels,avg_ref_elecs)); % ID of the average ref electrodes
nchns = length(elec_labels);

anode_col = [];
cathode_col = [];

for chnii = 1:nchns

    anode_region = types.hdmf_common.DynamicTableRegion(...
        'data',elec_ids(chnii),...
        'description','first channel of reference',...
        'table',electable_link);

    cathode_region = types.hdmf_common.DynamicTableRegion(...
        'data',avg_ref_ids,...
        'description','mean of channels subtracted from anode',...
        'table',electable_link);

    anode_col = [anode_col; anode_region];
    cathode_col = [cathode_col; cathode_region];
end

id = types.hdmf_common.VectorIndex('data', [1:nchns]-1);
anode_acq_col = repmat(lfp_link,nchns,1);
cathode_acq_col = repmat(lfp_link,nchns,1);
anode_acq_chns = elec_ids;
cathode_acq_chns = repmat(avg_ref_ids,nchns,1);
label_col = elec_labels;

avg_ref_table = types.ndx_bipolar_scheme.BipolarSchemeTable('anodes',anode_col,'cathodes',cathode_col,'anodes_index',id,'cathodes_index',id);
types.util.dynamictable.addRawData(avg_ref_table,'anode_acquisition',anode_acq_col);
avg_ref_table.vectordata.get('anode_acquisition').description = 'original acquisition anode channel data is derived from';
types.util.dynamictable.addRawData(avg_ref_table,'anode_acquisition_channels',anode_acq_chns);
avg_ref_table.vectordata.get('anode_acquisition_channels').description = 'channel indices in the acquisition data is derived from';
types.util.dynamictable.addRawData(avg_ref_table,'cathode_acquisition',cathode_acq_col);
avg_ref_table.vectordata.get('cathode_acquisition').description = 'original acquisition cathode channel data is derived from';
types.util.dynamictable.addRawData(avg_ref_table,'cathode_acquisition_channels',cathode_acq_chns);
avg_ref_table.vectordata.get('cathode_acquisition_channels').description = 'channel indices in the acquisition data is derived from';
types.util.dynamictable.addRawData(avg_ref_table,'label',label_col);
avg_ref_table.vectordata.get('label').description = 'name to assign the channel';

%% Bipolar reference
elec_labels = {'LDa1-LDa2','LDa2-LDa3','LDa3-LDa4'};
primary_elec_labels = {'LDa1','LDa2','LDa3','LDa4'};
elec_ids = [0:3]; % ID of the electrode in the primary electrode table
nchns = length(elec_labels);

anode_col = [];
cathode_col = [];

for chnii = 1:nchns

    this_label = elec_labels{chnii};
    components = split(this_label,'-');

    chn1_idx = find( strcmp( primary_elec_labels,components{1} ) ) -1;
    chn2_idx = find( strcmp( primary_elec_labels,components{2} ) ) -1;

    anode_region = types.hdmf_common.DynamicTableRegion(...
        'data',chn1_idx,...
        'description','first channel',...
        'table',electable_link);

    cathode_region = types.hdmf_common.DynamicTableRegion(...
        'data',chn2_idx,...
        'description','second channel',...
        'table',electable_link);

    anode_col = [anode_col; anode_region];
    cathode_col = [cathode_col; cathode_region];
end

id = types.hdmf_common.VectorIndex('data', [1:nchns]-1);
anode_acq_col = repmat(lfp_link,nchns,1);
cathode_acq_col = repmat(lfp_link,nchns,1);
anode_acq_chns = [anode_col.data];
cathode_acq_chns = [cathode_col.data];
label_col = elec_labels;

bipolar_ref_table = types.ndx_bipolar_scheme.BipolarSchemeTable('anodes',anode_col,'cathodes',cathode_col,'anodes_index',id,'cathodes_index',id);
types.util.dynamictable.addRawData(bipolar_ref_table,'anode_acquisition',anode_acq_col);
bipolar_ref_table.vectordata.get('anode_acquisition').description = 'original acquisition anode channel data is derived from';
types.util.dynamictable.addRawData(bipolar_ref_table,'anode_acquisition_channels',anode_acq_chns);
bipolar_ref_table.vectordata.get('anode_acquisition_channels').description = 'channel indices in the acquisition data is derived from';
types.util.dynamictable.addRawData(bipolar_ref_table,'cathode_acquisition',cathode_acq_col);
bipolar_ref_table.vectordata.get('cathode_acquisition').description = 'original acquisition cathode channel data is derived from';
types.util.dynamictable.addRawData(bipolar_ref_table,'cathode_acquisition_channels',cathode_acq_chns);
bipolar_ref_table.vectordata.get('cathode_acquisition_channels').description = 'channel indices in the acquisition data is derived from';
types.util.dynamictable.addRawData(bipolar_ref_table,'label',label_col);
bipolar_ref_table.vectordata.get('label').description = 'name to assign the channel';

Thank you again

bendichter commented 3 years ago

@sportnoah14 I have added example usage in MATLAB to the readme. You should be able to get this to work by copy/pasting and changing just a few variables. The one thing we are missing is an explicit link to the source data. I think I might just add this as an optional attribute attached to the bipolar scheme object.

bendichter commented 3 years ago

@sportnoah14 I have updated the schema in a few minor ways. You may need to make some minor adjustments to your code, but you can now link to a source electrical series and there is now example code to show you exactly how to use the extension in MatNWB

nmarkowitz commented 3 years ago

@bendichter this looks great. With this I can do precisely what I want. The only difference I may make that's different from your example code is placing the NdxBipolarScheme object in /processing as I'm using it as a derivative of the raw trace and do not want to store a whole new ElectricalSeries. Let me know if you have an opinion on this. If this sounds ok to you then you can close this issue. Thanks again!

bendichter commented 3 years ago

@sportnoah14 yes, I think it's fine to put this in /processing/