nwbextractors: handling offsets in the electrode table

[miketrumpis]

Hi, these lines are leading to a TypeError for me when multiplying the offsets by a scalar (1e6). I think when grabbing offset = self._nwbfile.electrodes["offset"] you actually need to do a slice syntax [:] to get a regular ndarray. I.e. the fix would be

offset = self._nwbfile.electrodes["offset"][:]

This is with HDMF 3.4.6 and PyNWB 2.1.0

https://github.com/SpikeInterface/spikeinterface/blob/a21fac0b91c2502c6f50e609b85b753be26ff551/spikeinterface/extractors/nwbextractors.py#L169-L172

[miketrumpis]

After thinking some more, I'm not sure it's safe to assume that the main electrode table has metadata that applies to all derived ElectricalSeries. As I understand, there is only one electrodes table in the NWB file, and any ElectricalSeries can refer to a range of rows in that table (that may well be my misunderstanding!). But if there's an offset meant to apply to the raw ADC counts, it would certainly not apply to a filtered version of that signal.

[alejoe91]

@CodyCBakerPhD @bendichter @h-mayorquin

What do you guys think?

[CodyCBakerPhD]

Offsets were added to all TimeSeries children a while back as an attribute to the object alongside the conversion factor (and channel-wise conversion factor specific to ElectricalSeries). The whole offsets added to the electrode table thing was the pre-2.1.0 PyNWB workaround, so if you're pinning SI to >=2.1.0 you can just do away with the extra column here.

[CodyCBakerPhD]

That said, if channel-wise offsets are also requested for ElectricalSeries, I'd submit an issue for that on https://github.com/NeurodataWithoutBorders/nwb-schema/issues

[miketrumpis]

@CodyCBakerPhD I have wished for that exact thing, and hadn't gotten around to finding out the right venue to request it. TY!

[h-mayorquin]

Can you show the trace @miketrumpis (the type of error) or even better a minimal example? I was the last one that tamper with that function so I might be able to patch it quickly.

[CodyCBakerPhD]

I have wished for that exact thing, and hadn't gotten around to finding out the right venue to request it. TY!

Cool - out of curiosity (and since I'm probably the one who will make the changes upstream), what format are you using that has different offsets per channel and would you be willing to share some files that showcase that property in their headers?

[h-mayorquin]

After thinking some more, I'm not sure it's safe to assume that the main electrode table has metadata that applies to all derived ElectricalSeries. As I understand, there is only one electrodes table in the NWB file, and any ElectricalSeries can refer to a range of rows in that table (that may well be my misunderstanding!). But if there's an offset meant to apply to the raw ADC counts, it would certainly not apply to a filtered version of that signal.

You are correct that we are making assumption here that every the ElectricalSeries corresponds channel to channel with the electrode table. As this is not the case, we should use the electrodes property of the ElectricalSeries object (a region mapped to the dynamical table) to disambiguate and map properties correctly: https://pynwb.readthedocs.io/en/stable/pynwb.ecephys.html#pynwb.ecephys.ElectricalSeries

[h-mayorquin]

I am linking here the previous issue for backwards following and tracking: https://github.com/SpikeInterface/spikeinterface/pull/846

[miketrumpis]

Can you show the trace @miketrumpis (the type of error) or even better a minimal example? I was the last one that tamper with that function so I might be able to patch it quickly.

"Minimal" is hard with NWB haha.. still, this probably has more than necessary. I just spun up a new virtualenv with

pip install numpy spikeinterface pynwb

import numpy as np
import os
from datetime import datetime
from pynwb import NWBHDF5IO, NWBFile, get_manager
from pynwb.ecephys import ElectricalSeries, FilteredEphys

electrode_ids = np.arange(1, 50)
np.random.shuffle(electrode_ids)

if os.path.exists('nwb-demo.nwb'):
    os.unlink('nwb-demo.nwb')

nwb = NWBFile('demo', 'demo', datetime.now().astimezone())
dev = nwb.create_device(name='recorder')
grp = nwb.create_electrode_group('electrode', device=dev, location='brain', description='fake')
# Minimum necessary columns
info = dict(x=np.nan, y=np.nan, z=np.nan, imp=np.nan, group=grp, location='brain', filtering='none')
offsets = np.random.rand(len(electrode_ids)) * 50
# Extra "offset" column
nwb.add_electrode_column(name='offset', index=False, description='Level offset')
for id in electrode_ids:
    nwb.add_electrode(id=id, **info, offset=offsets[id - 1])
# Make an arbitrary region in the table to apply to the ElectricalSeries
electrode_region = nwb.create_electrode_table_region(list(range(15, 20)), 'record electrodes')
ephys = ElectricalSeries('raw_acq', np.random.randn(1000, 5), electrode_region, starting_time=0.0, rate=100.0)
nwb.add_acquisition(ephys)
with NWBHDF5IO('nwb-demo.nwb', mode='w') as io:
    io.write(nwb)

from spikeinterface.extractors.nwbextractors import NwbRecordingExtractor

ex = NwbRecordingExtractor('nwb-demo.nwb', electrical_series_name='raw_acq')

ex.get_channel_offsets()

(si-mwe) denali:notebooks(main)$ python mwe.py 
Traceback (most recent call last):
  File "/Users/mike/work//notebooks/mwe.py", line 33, in <module>
    ex = NwbRecordingExtractor('nwb-demo.nwb', electrical_series_name='raw_acq')
  File "/Users/mike/.pyenv/versions/si-mwe/lib/python3.9/site-packages/spikeinterface/extractors/nwbextractors.py", line 173, in __init__
    self.set_channel_offsets(offset * 1e6)
TypeError: unsupported operand type(s) for *: 'VectorData' and 'float'

[miketrumpis]

I have wished for that exact thing, and hadn't gotten around to finding out the right venue to request it. TY!

Cool - out of curiosity (and since I'm probably the one who will make the changes upstream), what format are you using that has different offsets per channel and would you be willing to share some files that showcase that property in their headers?

I'm afraid it is an experimental system. What specifically would be helpful to you?

[CodyCBakerPhD]

I'm afraid it is an experimental system. What specifically would be helpful to you?

Ah, I see. What would have been helpful is the structure of the header information that describes all these channel properties but if that isn't set in stone yet I'll just emulate it at the SI Recording level.