Skip to content

Fix Neuropixels handling in SpikeGadgetsRawIO #1847

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
h-mayorquin wants to merge 1 commit into
base: master
Choose a base branch
from
+66 −41
Open

Fix Neuropixels handling in SpikeGadgetsRawIO #1847

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 40 additions & 41 deletions neo/rawio/spikegadgetsrawio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,12 +140,6 @@ def _parse_header(self):
"than the number of channels in the hardware configuration"
)

# as spikegadgets change we should follow this
try:
num_chan_per_chip = int(sconf.attrib["chanPerChip"])
except KeyError:
num_chan_per_chip = 32 # default value for Intan chips

# explore sub stream and count packet size
# first bytes is 0x55
packet_size = 1
Expand Down Expand Up @@ -219,56 +213,61 @@ def _parse_header(self):
signal_streams.append((stream_name, stream_id, buffer_id))
self._mask_channels_bytes[stream_id] = []

# we can only produce these channels for a subset of spikegadgets setup. If this criteria isn't
# true then we should just use the raw_channel_ids and let the end user sort everything out
if num_ephy_channels % num_chan_per_chip == 0:
all_hw_chans = [int(schan.attrib["hwChan"]) for trode in sconf for schan in trode]
missing_hw_chans = set(range(num_ephy_channels)) - set(all_hw_chans)
channel_ids = self._produce_ephys_channel_ids(
num_ephy_channels_xml, num_chan_per_chip, missing_hw_chans
)
raw_channel_ids = False
else:
raw_channel_ids = True
# Channel id resolution.
# - Default: take each spike channel's id from its hwChan attribute in the XML.
# That's the identifier produced by the acquisition hardware, and it's the right
# thing for Neuropixels and most modern setups.
# - Intan special case (`device` is "intan" or absent for legacy files): recompute
# the ids from the chip layout instead. Intan recordings lay out the binary
# stream in a chip-interleaved order (all chips' channel 0, then all chips'
# channel 1, ...), and the synthesised labels reflect that ordering. See
# `_produce_ephys_channel_ids` and issue #1215 for the rationale. Only valid
# when chanPerChip divides evenly into the total channel count.
# Default: use hwChan directly from the XML.
channel_id_and_name = [
(schan.attrib["hwChan"], f"trode{trode.attrib['id']}chan{schan.attrib['hwChan']}")
for trode in sconf
for schan in trode
]
# Intan special case: replace with synthesised chip-interleaved ids when applicable.
spike_device = sconf.attrib.get("device")
if spike_device in (None, "intan"):
intan_chans_per_chip = int(sconf.attrib.get("chanPerChip", 32)) # RHD2132 default for legacy
if intan_chans_per_chip > 0 and num_ephy_channels % intan_chans_per_chip == 0:
hw_chans_in_xml = {int(schan.attrib["hwChan"]) for trode in sconf for schan in trode}
missing_hw_chans = set(range(num_ephy_channels)) - hw_chans_in_xml
synthesised_ids = self._produce_ephys_channel_ids(
num_ephy_channels_xml, intan_chans_per_chip, missing_hw_chans
)
channel_id_and_name = [(str(cid), f"hwChan{cid}") for cid in synthesised_ids]

chan_ind = 0
self.is_scaleable = all("spikeScalingToUv" in trode.attrib for trode in sconf)
if not self.is_scaleable:
self.logger.warning(
"Unable to read channel gain scaling (to uV) from .rec header. Data has no physical units!"
)

for trode in sconf:
trode_per_channel = [trode for trode in sconf for _ in trode]
for chan_ind, trode in enumerate(trode_per_channel):
if "spikeScalingToUv" in trode.attrib:
gain = float(trode.attrib["spikeScalingToUv"])
units = "uV"
else:
gain = 1 # revert to hardcoded gain
units = ""

for schan in trode:
# Here we use raw ids if necessary for parsing (for some neuropixel recordings)
# otherwise we default back to the raw hwChan IDs
if raw_channel_ids:
name = "trode" + trode.attrib["id"] + "chan" + schan.attrib["hwChan"]
chan_id = schan.attrib["hwChan"]
else:
chan_id = str(channel_ids[chan_ind])
name = "hwChan" + chan_id

offset = 0.0
buffer_id = ""
signal_channels.append(
(name, chan_id, self._sampling_rate, "int16", units, gain, offset, stream_id, buffer_id)
)

chan_mask = np.zeros(packet_size, dtype="bool")
num_bytes = packet_size - 2 * num_ephy_channels + 2 * chan_ind
chan_mask[num_bytes] = True
chan_mask[num_bytes + 1] = True
self._mask_channels_bytes[stream_id].append(chan_mask)
chan_id, name = channel_id_and_name[chan_ind]
offset = 0.0
buffer_id = ""
signal_channels.append(
(name, chan_id, self._sampling_rate, "int16", units, gain, offset, stream_id, buffer_id)
)

chan_ind += 1
chan_mask = np.zeros(packet_size, dtype="bool")
num_bytes = packet_size - 2 * num_ephy_channels + 2 * chan_ind
chan_mask[num_bytes] = True
chan_mask[num_bytes + 1] = True
self._mask_channels_bytes[stream_id].append(chan_mask)

# make mask as array (used in _get_analogsignal_chunk(...))
self._mask_streams = {}
Expand Down
26 changes: 26 additions & 0 deletions neo/test/rawiotest/test_spikegadgetsrawio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ class TestSpikeGadgetsRawIO(
"spikegadgets/20210225_em8_minirec2_ac.rec",
"spikegadgets/W122_06_09_2019_1_fromSD.rec",
"spikegadgets/SpikeGadgets_test_data_2xNpix1.0_20240318_173658.rec",
"spikegadgets/neuropixels2_4shank/20260122_134412_merged_cropped_1min_NP2.rec",
]

def test_parse_header_missing_channels(self):
Expand Down Expand Up @@ -53,6 +54,31 @@ def test_parse_header_missing_channels(self):
# fmt: on
)

def test_neuropixels_uses_hwchan_ids(self):
# Regression test: Neuropixels recordings must expose hwChan ids, not synthesised
# Intan chip ids. The original reports of the related ZeroDivisionError on NP2
# (chanPerChip="0") are in issues #1844 and #1810 and are exercised through the
# NP2 fixture in `entities_to_test`. This test guards against the separate NP1
# silent-misrouting case that surfaced while investigating those issues: NP1
# files like this one have chanPerChip="32" despite the probe not using Intan
# chips, which made the old divisibility check pass and silently produced
# synthesised ids. The fix relies on SpikeConfiguration.device to detect
# non-Intan setups regardless of chanPerChip.
file_path = Path(
self.get_local_path("spikegadgets/SpikeGadgets_test_data_2xNpix1.0_20240318_173658.rec")
)
reader = SpikeGadgetsRawIO(filename=file_path)
reader.parse_header()

trodes_mask = reader.header["signal_channels"]["stream_id"] == "trodes"
trodes_ids = list(reader.header["signal_channels"]["id"][trodes_mask])

# Real hwChan values from the XML, in document order. The synthesised path with
# chanPerChip=32 would produce ['0', '32', '64', '96'] instead.
self.assertEqual(trodes_ids[:4], ["735", "734", "671", "670"])
self.assertEqual(trodes_ids[-4:], ["97", "96", "33", "32"])
self.assertEqual(len(trodes_ids), 768)

def test_opening_gibberish_file(self):
"""Test that parsing a file without </Configuration> raises ValueError instead of infinite loop."""
# Create a temporary file with gibberish content that doesn't have the required tag
Expand Down
Loading