Spike Sorting with SpikeInterface

Import Libraries

Import the modules required for this notebook

import numpy as np
from matplotlib import pyplot as plt
import spikeinterface.full as si
import probeinterface as pi
from probeinterface.plotting import plot_probe
import scipy

%matplotlib inline

Download Data

Download the Data required for this notebook

import os
import requests
import zipfile

url = "https://uni-bonn.sciebo.de/s/qrXWd4ZHozGdJwN"
response = requests.get(url + "/download")
print("Downloading Data ...")
with open("openephys_recording.zip", "wb") as file:
    file.write(response.content)

with zipfile.ZipFile("openephys_recording.zip", "r") as zip_ref:
    zip_ref.extractall("")
os.remove("openephys_recording.zip")
print("Done!")

SpikeInterface provides a unified interface for reading, handling, and processing extracellular electrophysiological recordings from a wide variety of file formats. In this section, you will learn how to load recordings, inspect their properties, and extract data segments for further analysis. This is the first step in any spike sorting workflow, as it allows you to access and manipulate the raw voltage traces recorded from your experiment.

Run the cell below to define the data directory

Exercises

data_dir = "openephys_raw/2023-08-23_15-56-05"

Example: Load the recording stream with the id "1", which contains the auxiliary channels, assign the returned RecordingExtractor to a variable called aux_rec and print it.

rec = si.read_openephys(data_dir, stream_id="1")
rec

Exercise: Load the stream with the id "0" and assign the returned RecordingExtractor to a variable called rec.

rec = si.read_openephys(data_dir, stream_id="0")
rec

Exercise: Get the duration of the recording

rec.get_duration()

300.0

Exercise: Get the sampling frequency of the recording and assign it to a variable fs

fs = rec.get_sampling_frequency()

Example: Use the frame_slice method to select the first 10 seconds of the recording rec. Assign the returned object to a new variable seg and get it’s duration.

seg = rec.frame_slice(start_frame=0 * fs, end_frame=10 * fs)
seg.get_duration()

10.000000000000027

Exercise: Use the frame_slice method to select the segment between 0.5 and 0.7 seconds. Assign the returned object to a new variable seg and get it’s duration.

seg = rec.frame_slice(start_frame=0.5 * fs, end_frame=0.7 * fs)
seg.get_duration()

0.20000000000001594

Exercise: Use seg.get_traces() to get the raw voltage traces and assign them to a variable traces.

traces = seg.get_traces()

Example: Plot the extracted traces for the first 3 channels

plt.plot(traces[:, :3])

Exercise: Plot the extracted traces for the last 5 channels

Solution

plt.plot(traces[:, -5:])

Accurate probe information is essential for interpreting extracellular recordings, as it defines the spatial arrangement and wiring of the electrodes. ProbeInterface, integrated with SpikeInterface, allows you to load, visualize, and configure probe layouts from a wide range of manufacturers. In this section, you will learn how to select a probe, visualize its geometry, remap its wiring, and attach it to your recording for downstream analysis.

Exercises

Example: Get the "NP1010" probe from "imec".

manufacturer = "imec"
probe_name = "NP1010"
probe = pi.get_probe(manufacturer, probe_name)
probe

Probe - imec - NP1010 - 960ch - 1shanks

Example: Plot the probe and adjust the axis limits to display the probe’s tip

fig, ax = plt.subplots()
plot_probe(probe, ax=ax)
ax.set(xlim=(-50, 100), ylim=(-200, 150))

[(-50.0, 100.0), (-200.0, 150.0)]

Exercise: Get the "1S256" probe from "plexon"

manufacturer = 'plexon'
probe_name = '1S256'
probe = pi.get_probe(manufacturer, probe_name)
probe

Probe - Plexon-SiNAPS - 1S256 - 256ch - 1shanks

Exercise: Plot the probe and adjust the axis limits to display the probe’s tip

Solution

fig, ax = plt.subplots()
plot_probe(probe, ax=ax)
ax.set(xlim=(-50, 75), ylim=(-150, 150))

[(-50.0, 75.0), (-150.0, 150.0)]

Exercise: Go to the ProbeInterface library and select a new manufacturer and probe from there. Then, get the probe and plot it

Solution

probe = pi.get_probe("neuronexus", "A4x8-5mm-100-200-177")
fig, ax = plt.subplots()
plot_probe(probe, ax=ax)
ax.set(xlim=(-50, 50), ylim=(-80, 150))

[(-50.0, 50.0), (-80.0, 150.0)]

Example: Load the "ASSY-236-H5" probe from "cambridgeneurotech" and remap the channels for the "cambridgeneurotech_mini-amp-64" wiring pathway. Then, plot the probe with the contact and device indices and adjust the axis limit to view the tip of the probe.

probe = pi.get_probe("cambridgeneurotech", "ASSY-236-H5")
probe.wiring_to_device("cambridgeneurotech_mini-amp-64")

fig, ax = plt.subplots()
plot_probe(probe, ax=ax, with_contact_id=True, with_device_index=True)
ax.set(xlim=(-25, 50), ylim=(-50, 150))

[(-25.0, 50.0), (-50.0, 150.0)]

Exercise: Load the "ASSY-236-H5" probe from "cambridgeneurotech" and remap the channels for the "ASSY-156>RHD2164" wiring pathway. Then, plot the probe with the contact and device indices and adjust the axis limit to view the tip of the probe.

Solution

probe = pi.get_probe("cambridgeneurotech", "ASSY-236-H5")
probe.wiring_to_device("ASSY-156>RHD2164")

fig, ax = plt.subplots()
plot_probe(probe, ax=ax, with_contact_id=True, with_device_index=True)
ax.set(xlim=(-25, 50), ylim=(-50, 150))

[(-25.0, 50.0), (-50.0, 150.0)]

Exercise: List all of the available wiring pathways. Select one and apply it to the "ASSY-236-H5" probe from "cambridgeneurotech" and plot the probe (NOTE: not all wiring pathways are compatible with every device).

Solution

pi.get_available_pathways()
probe = pi.get_probe("cambridgeneurotech", "ASSY-236-H5")
probe.wiring_to_device("ASSY-77>Adpt.A64-Om32_2x-sm-NN>two_RHD2132")

fig, ax = plt.subplots()
plot_probe(probe, ax=ax, with_contact_id=True, with_device_index=True)
ax.set(xlim=(-25, 50), ylim=(-50, 150))

[(-25.0, 50.0), (-50.0, 150.0)]

Exercise: Add the probe to the recording rec using .set_probe() and then call .get_probe() method to verify the probe info was added to the recording.

rec = rec.set_probe(probe)
rec.get_probe()

Probe - 64ch - 1shanks

Spike sorting is the process of detecting and classifying action potentials (spikes) from extracellular recordings, assigning them to putative single units (neurons). SpikeInterface provides a unified API to run a variety of spike sorting algorithms (sorters) and manage their outputs. In this section, you will learn how to list available sorters, run a sorter on your data, and inspect the results, including the number of detected units and their spike trains.

Run the code below to load the recording, and add the probe information to it.

Exercises

rec = si.read_openephys(data_dir, stream_id="0")
probe = pi.get_probe("cambridgeneurotech", "ASSY-236-H5")
probe.wiring_to_device("cambridgeneurotech_mini-amp-64")
rec = rec.set_probe(probe)

Exercise: List all of the installed spike sorters

si.installed_sorters()

['simple', 'spykingcircus2', 'tridesclous2']

Example: Select the first 5 seconds of the recording rec and run the spike sorter "tridesclous2" on the recording segment seg.

t_start = 0
t_stop = 5
fs = rec.get_sampling_frequency()
seg = rec.frame_slice(start_frame=t_start * fs, end_frame=t_stop * fs)
sorting_TDC = si.run_sorter(
    sorter_name="tridesclous2",
    recording=seg,
    folder="results_TDC",
    remove_existing_folder=True,
)

Exercise: Get the number of units in sorting_TDC.

sorting_TDC.get_num_units()

29

Exercise: Get the unit ids in sorting_TDC.

sorting_TDC.get_unit_ids()

array([ 3,  8,  9, 10, 11, 12, 14, 16, 18, 19, 20, 24, 25, 30, 33, 39, 40,
       41, 43, 44, 45, 47, 48, 49, 52, 53, 65, 67, 69])

Exercise: Get the spike train for any of the units.

sorting_TDC.get_unit_spike_train(10, return_times=True)

array([353.0968    , 353.1337    , 353.6497    , 353.67243333,
       353.6785    , 353.77126667, 353.78456667, 353.78456667,
       354.0273    , 354.62606667, 354.65626667, 354.65626667,
       354.66563333, 354.68123333, 354.70546667, 354.74256667,
       355.12943333, 355.13546667, 355.14583333, 355.1731    ,
       355.2428    , 355.27963333, 355.5034    , 355.5179    ,
       355.52286667, 355.54753333, 355.577     , 355.61273333,
       356.28896667, 356.37073333, 356.66496667, 356.67043333,
       356.67443333, 356.8355    , 357.4416    , 357.4605    ])

Exercise: Select the first 5 seconds of the recording rec and run the spike sorter "spykingcircus2" with apply_motion_correction=False on the recording segment seg. Then, print the number of identified units.

t_start = 0
t_stop = 5
fs = rec.get_sampling_frequency()
seg = rec.frame_slice(start_frame=t_start * fs, end_frame=t_stop * fs)
sorting_SPC = si.run_sorter(
    sorter_name="spykingcircus2",
    recording=seg,
    folder="results_TDC",
    remove_existing_folder=True,
    apply_motion_correction=False,
)
sorting_SPC.get_num_units()

spykingcircus2 could benefit from using torch. Consider installing it

55

Many spike sorters offer a range of parameters that can be tuned to optimize detection and clustering for your specific dataset. SpikeInterface allows you to query and set these parameters in a standardized way. In this section, you will learn how to retrieve default parameters for a sorter, modify detection and filtering settings, and run sorters with custom configurations to improve sorting performance.

Exercise: Get the default parameters for the spike sorter "tridesclous2"

si.get_default_sorter_params("tridesclous2")

{'apply_preprocessing': True,
 'apply_motion_correction': False,
 'motion_correction': {'preset': 'dredge_fast'},
 'cache_preprocessing': {'mode': 'memory',
  'memory_limit': 0.5,
  'delete_cache': True},
 'waveforms': {'ms_before': 0.5, 'ms_after': 1.5, 'radius_um': 120.0},
 'filtering': {'freq_min': 150.0,
  'freq_max': 5000.0,
  'ftype': 'bessel',
  'filter_order': 2},
 'detection': {'peak_sign': 'neg',
  'detect_threshold': 5,
  'exclude_sweep_ms': 1.5,
  'radius_um': 150.0},
 'selection': {'n_peaks_per_channel': 5000, 'min_n_peaks': 20000},
 'svd': {'n_components': 4},
 'clustering': {'recursive_depth': 5,
  'split_radius_um': 40.0,
  'clusterer': 'isosplit6',
  'clusterer_kwargs': {},
  'do_merge': True,
  'merge_kwargs': {'similarity_metric': 'l1',
   'num_shifts': 4,
   'similarity_thresh': 0.75}},
 'templates': {'ms_before': 2.0,
  'ms_after': 3.0,
  'max_spikes_per_unit': 400,
  'sparsity_threshold': 1.5,
  'min_snr': 2.5},
 'matching': {'method': 'tdc-peeler',
  'method_kwargs': {},
  'gather_mode': 'memory'},
 'job_kwargs': {},
 'save_array': True}

Example: Run the spike sorter "tridesclous2" on the first 5 seconds of the recording rec and set the detection parameter to {"detect_threshold":8}. Then, print the number of detected units.

t_start = 0
t_stop = 5
fs = rec.get_sampling_frequency()
seg = rec.frame_slice(start_frame=t_start * fs, end_frame=t_stop * fs)

detection = {"detect_threshold": 8}
sorting_TDC = si.run_sorter(
    sorter_name="tridesclous2",
    recording=seg,
    folder="results_TDC",
    remove_existing_folder=True,
    detection=detection,
)
sorting_TDC.get_num_units()

29

Exercise: Run the spike sorter "tridesclous2" on the first 10 seconds of the recording rec and set the filtering parameter to {"freq_min": 1000, "freq_max":9000}. Then, print the number of detected units.

t_start = 0
t_stop = 5
fs = rec.get_sampling_frequency()
seg = rec.frame_slice(start_frame=t_start * fs, end_frame=t_stop * fs)

filtering = {"freq_min": 1000, "freq_max": 9000}
sorting_TDC = si.run_sorter(
    sorter_name="tridesclous2",
    recording=seg,
    folder="results_TDC",
    remove_existing_folder=True,
    filtering=filtering,
)
sorting_TDC.get_num_units()

18

Exercise: Run the spike sorter "tridesclous2" on the segment seg with svd={"n_components":3}

sorting_TDC = si.run_sorter(
    sorter_name="tridesclous2",
    recording=seg,
    folder="results_TDC",
    remove_existing_folder=True,
    svd={"n_components": 3},
)
sorting_TDC.get_num_units()

27