NEW: probe reader for Neuropixels 1.0 in SpikeGadgets .rec file.

src/probeinterface/__init__.py

@@ -16,6 +16,8 @@
     write_imro,
     read_BIDS_probe,
     write_BIDS_probe,
+    read_spikegadgets,
+    parse_spikegadgets_header,
     read_spikeglx,
     parse_spikeglx_meta,
     get_saved_channel_indices_from_spikeglx_meta,

src/probeinterface/io.py

@@ -19,6 +19,7 @@
 from collections import OrderedDict
 from packaging.version import Version, parse
 import numpy as np
+from xml.etree import ElementTree
 
 from . import __version__
 from .probe import Probe
@@ -1263,12 +1264,188 @@ def write_imro(file: str | Path, probe: Probe):
         f.write("".join(ret))
 
 
+def read_spikegadgets(file: str | Path) -> ProbeGroup:
+    """
+    Find active channels of the given Neuropixels probe from a SpikeGadgets .rec file.
+    SpikeGadgets headstages support up to three Neuropixels 1.0 probes (as of March 28, 2024),
+    and information for all probes will be returned in a ProbeGroup object.
+
+    Within the SpikeConfiguration header element (sconf), there is a SpikeNTrode element
+    for each electrophysiology channel that contains information relevant to scaling and
+    otherwise displaying the information from that channel, as well as the id of the electrode
+    from which it is recording ('id').
+
+    Nested within each SpikeNTrode element is a SpikeChannel element with information about
+    the electrode dynamically connected to that channel.  This contains information relevant
+    for spike sorting, i.e., its spatial location along the probe shank and the hardware channel
+    to which it is connected.
+
+    Excerpt of a sample SpikeConfiguration element:
+
+    <SpikeConfiguration chanPerChip="1889715760" device="neuropixels1" categories="">
+        <SpikeNTrode viewLFPBand="0"
+            viewStimBand="0"
+            id="1384"  # @USE: The first digit is the probe number; the last three digits are the electrode number
+            lfpScalingToUv="0.018311105685598315"
+            LFPChan="1"
+            notchFreq="60"
+            rawRefOn="0"
+            refChan="1"
+            viewSpikeBand="1"
+            rawScalingToUv="0.018311105685598315"  # For Neuropixels 1.0, raw and spike scaling are identical
+            spikeScalingToUv="0.018311105685598315"  # Extracted when reading the raw data
+            refNTrodeID="1"
+            notchBW="10"
+            color="#c83200"
+            refGroup="2"
+            filterOn="1"
+            LFPHighFilter="200"
+            moduleDataOn="0"
+            groupRefOn="0"
+            lowFilter="600"
+            refOn="0"
+            notchFilterOn="0"
+            lfpRefOn="0"
+            lfpFilterOn="0"
+            highFilter="6000"
+        >
+            <SpikeChannel thresh="60"
+                coord_dv="-480"  # @USE: dorsal-ventral coordinate in um (in pairs for staggered probe)
+                spikeSortingGroup="1782505664"
+                triggerOn="1"
+                stimCapable="0"
+                coord_ml="3192"  # @USE: medial-lateral coordinate in um
+                coord_ap="3700"  # doesn't vary, assuming the shank's flat face is along the ML axis
+                maxDisp="400"
+                hwChan="735"  # @USE: unique device channel that is reading from electrode
+            />
+        </SpikeNTrode>
+        ...
+    </SpikeConfiguration>
+
+    Parameters
+    ----------
+    file : Path or str
+        The .rec file path
+
+    Returns
+    -------
+    probe_group : ProbeGroup object
+
+    """
+    # ------------------------- #
+    #     Npix 1.0 constants    #
+    # ------------------------- #
+    TOTAL_NPIX_ELECTRODES = 960
+    MAX_ACTIVE_CHANNELS = 384
+    CONTACT_WIDTH = 16  # um
+    CONTACT_HEIGHT = 20  # um
+    # ------------------------- #
+
+    # Read the header and get Configuration elements
+    header_txt = parse_spikegadgets_header(file)
+    root = ElementTree.fromstring(header_txt)
+    hconf = root.find("HardwareConfiguration")
+    sconf = root.find("SpikeConfiguration")
+
+    # Get number of probes present (each has its own Device element)
+    probe_configs = [device for device in hconf if device.attrib['name'] == 'NeuroPixels1']
+    n_probes = len(probe_configs)
+    print(n_probes, "Neuropixels 1.0 probes found:")
+
+    # Container to store Probe objects
+    probe_group = ProbeGroup()
+
+    for curr_probe in range(1, n_probes + 1):
+        print(f'Reading probe{curr_probe}...', flush=True, end='')
+        probe_config = probe_configs[curr_probe - 1]
+
+        # Get number of active channels from probe Device element
+        active_channel_str = [option for option in probe_config if option.attrib['name'] == 'channelsOn'][0].attrib['data']
+        active_channels = [int(ch) for ch in active_channel_str.split(' ') if ch]
+        n_active_channels = sum(active_channels)
+        assert len(active_channels) == TOTAL_NPIX_ELECTRODES
+        assert n_active_channels <= MAX_ACTIVE_CHANNELS
+
+        # Find all channels/electrodes that belong to the current probe
+        contact_ids = []
+        device_channels = []
+        positions = np.zeros((n_active_channels, 2))
+
+        nt_i = 0  # Both probes are in sconf, so need an independent counter of probe electrodes while iterating through
+        for ntrode in sconf:
+            electrode_id = ntrode.attrib['id']
+            if int(electrode_id[0]) == curr_probe:  # first digit of electrode id is probe number
+                contact_ids.append(electrode_id)
+                positions[nt_i, :] = (ntrode[0].attrib['coord_ml'], ntrode[0].attrib['coord_dv'])
+                device_channels.append(ntrode[0].attrib['hwChan'])
+                nt_i += 1
+        assert len(contact_ids) == n_active_channels
+
+        # Construct Probe object
+        probe = Probe(ndim=2, si_units="um", model_name="Neuropixels 1.0", manufacturer="IMEC")
+        probe.set_contacts(
+            contact_ids=contact_ids,
+            positions=positions,
+            shapes="square",
+            shank_ids=None,
+            shape_params={"width": CONTACT_WIDTH, "height": CONTACT_HEIGHT},
+        )
+
+        # Wire it (i.e., point contact/electrode ids to corresponding hardware/channel ids)
+        probe.set_device_channel_indices(device_channels)
+
+        # Create a nice polygon background when plotting the probes
+        x_min = positions[:, 0].min()
+        x_max = positions[:, 0].max()
+        x_mid = 0.5 * (x_max + x_min)
+        y_min = positions[:, 1].min()
+        y_max = positions[:, 1].max()
+        polygon_default = [
+            (x_min - 20, y_min - CONTACT_HEIGHT / 2),
+            (x_mid, y_min - 100),
+            (x_max + 20, y_min - CONTACT_HEIGHT / 2),
+            (x_max + 20, y_max + 20),
+            (x_min - 20, y_max + 20),
+        ]
+        probe.set_planar_contour(polygon_default)
+
+        # If there are multiple probes, they must be shifted such that they don't occupy the same coordinates.
+        probe.move([250 * (curr_probe - 1), 0])
+
+        # Add the probe to the probe container
+        print(probe.get_contact_count(), 'active channels found.')
+        probe_group.add_probe(probe)
+
+    return probe_group
+
+
+def parse_spikegadgets_header(file: str | Path) -> str:
+    """
+    Parse file (SpikeGadgets .rec format) into a string until "</Configuration>",
+    which is the last tag of the header, after which the binary data begins.
+    """
+    header_size = None
+    with open(file, mode="rb") as f:
+        while True:
+            line = f.readline()
+            if b"</Configuration>" in line:
+                header_size = f.tell()
+                break
+
+        if header_size is None:
+            ValueError("SpikeGadgets: the xml header does not contain '</Configuration>'")
+
+        f.seek(0)
+        return f.read(header_size).decode("utf8")
+
+
 def read_spikeglx(file: str | Path) -> Probe:
     """
     Read probe position for the meta file generated by SpikeGLX
 
     See http://billkarsh.github.io/SpikeGLX/#metadata-guides for implementation.
-    The x_pitch/y_pitch/width are set automatically depending the NP version.
+    The x_pitch/y_pitch/width are set automatically depending on the NP version.
 
     The shape is auto generated as a shank.
 
@@ -1333,7 +1510,7 @@ def read_spikeglx(file: str | Path) -> Probe:
 def parse_spikeglx_meta(meta_file: str | Path) -> dict:
     """
     Parse the "meta" file from spikeglx into a dict.
-    All fiields are kept in txt format and must also parsed themself.
+    All fields are kept in txt format and must also be parsed themselves.
     """
     meta_file = Path(meta_file)
     with meta_file.open(mode="r") as f: