In this chapter, we will use SpikeInterface to implement a standard spike sorting pipeline from start to finish.
Here we will go through the key steps at a high level, and develop the theory while adding more advanced methods and use cases in the following chapters.
In this walkthrough, we will use a short example dataset available in the data/pipeline-walkthrough of the course repository. Alterantively, the data can be downloaded from here ADD ZENODO.
This recording is only 2.5 seconds long, from 384 channels on a 4-shank Neuropixels 2 probe acquired in SpikeGLX. We use two of the four shanks, with 192 recording channels on each.
We will load and preprocess the data, run spike sorting and compute quality metrics for assessing the sorting output.
Letβs get started!
First, we will load the data with SpikeInterface.
SpikeInterface can load data from many different of acquisition systems, using functions within its extractors module. As our data is from SpikeGLX, we will use the []read_spikeglx](https://github.com/neuroinformatics-unit/course-ephys-osss/tree/main/book/data/pipeline-walkthrough) method:
import spikeinterface.extractors as si_extractors
# This should be the path to the folder containing
# the data. On Windows, you may need to prepend
# r to the string e.g. r"C:\my\path"
path_to_data_folder = "data/pipeline-walkthrough"
recording = si_extractors.read_spikeglx(
path_to_data_folder,
stream_id="imec0.ap"
)
print(recording)SpikeGLXRecordingExtractor: 384 channels - 30.0kHz - 1 segments - 75,000 samples - 2.50s
int16 dtype - 54.93 MiBWe need to tell SpikeInterface where the data is on our machine, and what data streamπ‘ we want to load. If you are unsure of the stream, run the function without the stream_id argument and SpikeInterface will prompt you with the available streams.
Because we want to load the AP data stream for spike sorting, we pass "imec0.ap".
SpikeInterface uses βlazy loadingβ, which means at this stage no data is actually loaded into memory. Data is only loaded when we request it (e.g. for visualisation, or spike sorting). This makes operations very fast.
SpikeInterface extractor functions return the recording as a Python object. This object has many useful methods for exploring the recording as well as accessing the underlying data. In our case, recording is a SpikeGLXRecordingExtractor recording object. See the inheritanceπ‘ for details on inheritence in SpikeInterface.
You can also use Pythonβs dir() function to see all attributes and methods associated with any Python object.
print(
[method for method in dir(recording) if method.startswith("get")]
)['get_annotation', 'get_annotation_keys', 'get_binary_description', 'get_channel_gains', 'get_channel_groups', 'get_channel_ids', 'get_channel_locations', 'get_channel_offsets', 'get_channel_property', 'get_dtype', 'get_duration', 'get_end_time', 'get_memory_size', 'get_neo_io_reader', 'get_num_blocks', 'get_num_channels', 'get_num_frames', 'get_num_samples', 'get_num_segments', 'get_parent', 'get_preferred_mp_context', 'get_probe', 'get_probegroup', 'get_probes', 'get_property', 'get_property_keys', 'get_sampling_frequency', 'get_start_time', 'get_streams', 'get_time_info', 'get_times', 'get_total_duration', 'get_total_memory_size', 'get_total_samples', 'get_traces']This is a good time to explore the recording object using the API docs. Have a go at the following:
read_spikeglx function without stream_id to see what streams are on the recording.get_traces() to retrieve a segment of raw data as a NumPy array. What are its dimensions?dir() call above. Try out any functions that look interesting!There are many brands of recording probes used in extracellular electrophysiology (e.g. Neuropixels, Cambridge Neurotech or NeuroNexus ) which makes loading probes into a common framework complex.
Under the hood, SpikeInterface uses the ProbeInterface package to load many different types of probe. In our case, it automatically identifies and loads the Neuropixels 2 (NP2) probe.
It is important to check that the probe has loaded as expected, and that the geometry matches how it was set up in your acquisition software. We can plot the probe to confirm it loaded correctly:
from probeinterface.plotting import plot_probe
import matplotlib.pyplot as plt
probe = recording.get_probe()
plot_probe(probe)
plt.show()
As expected, we have a 4-shank Neurpixels probe in which all the recording channels are located in the middle two shanks.
Electrophysiological recordings are susceptible to contamination from multiple noise sources: electrical interference, mechanical noise caused by probe movement, and artefacts arising from hardware defects or poor contact.
As such, it is important to visualise the data, especially after preprocessing, to check the data quality. First, letsβs plot one second of data from our raw recording:
import spikeinterface.widgets as si_widgets
start_time = recording.get_times()[0]
si_widgets.plot_traces(
recording,
time_range=(start_time, start_time + 0.1),
return_in_uV=True,
order_channel_by_depth=True
)
plt.show()
return_in_uV to ensure the data is returned in microvolts, as by default most electrophysiological recordings are stored as unitless int16 valuesπ‘.
order_channel_by_depth is used as, depending on the probe, the default ordering of channels is not always by depth, meaning the order may be scrambled on the plot.
In general when working with raw data, you can use recording.get_channel_locations() to inspect the default channel ordering of the channels, and order_channels_by_depth() function can be used to get the indicies for reordering.
As we can see - the raw data looks terrible! It is severely contaminated by electrical noise. We will need to preprocess it to tidy it up.
Below, we will apply three key preprocessing steps:
phase_shift : corrects a hardware artefact in the Neuropixels probe which causes small shifts between the recording time across nearby channels.
bandpass_filter : removes high and low frequency oscillations from the data. It is applied to each channel separately.
commmon_reference : removes noise spikes that appear as verticle stripes across the data. It is applied across channels separately at each timepoint.
To preprocess the data, we will pass our recording object into functions that return a new, recording object. e.g. filtered_recording = bandpass_filter(recording). We can do this for every step in the preprocessing chain:
import spikeinterface.preprocessing as si_prepro
phase_shift_recording = si_prepro.phase_shift(recording)
filtered_recording = si_prepro.bandpass_filter(
phase_shift_recording, freq_min=300, freq_max=6000
)
cmr_recording = si_prepro.common_reference(
filtered_recording, operator="median"
)In SpikeInterface, the preprocessing is applied almost instantly due to the βlazyβ loading mentioned earlier. When we apply the preprocessing function, we only really just change the metadata associated with the recording, the data is not actually changed at this stage.
Only when we get the data (e.g. in plot_traces, or with get_traces() or when saving the data) is the data preprocessed in the requested chunk. This makes operating very fast!
Now we are ready to plot the data:
si_widgets.plot_traces(
cmr_recording,
time_range=(start_time, start_time + 0.1),
return_in_uV=True,
order_channel_by_depth=True
)
plt.show()
As you can see, it is much cleaner and we can see action potentials firing! However, there does seem to be an artefact, small gaps in the action potential waveform. We will investigate these in the next section.
You might notice the data in the above image looks slightly strange. All spikes have predictable gaps between them, we can see this clearly if we zoom into a single action potential waveform:
si_widgets.plot_traces(
cmr_recording,
time_range=(23.923, 23.928),
return_in_uV=True,
order_channel_by_depth=True,
)
ax = plt.gca()
ax.set_ylim(200, 300)
plt.show()
Our recording channels are located on two shanks, however we have plot all recording channels together. The channels on different shanks have the same y-position, so they are plot next to eachother, but different x-positions.
This means that the AP recorded on one shanks channels is interleaved with blank channels from the other shank.
We can solve this by spliting the recording by shank with the groups property. When SpikeInterface loads the data, it will automatically assign different shanks to different groups:
split_recording = cmr_recording.split_by("group")
print(
split_recording
){0: ChannelSliceRecording: 192 channels - 30.0kHz - 1 segments - 75,000 samples - 2.50s - int16 dtype
27.47 MiB, 1: ChannelSliceRecording: 192 channels - 30.0kHz - 1 segments - 75,000 samples - 2.50s - int16 dtype
27.47 MiB}We can plot both shanks separately:
fig, ax = plt.subplots(1, 2)
for shank_idx, shank_recording in split_recording.items():
si_widgets.plot_traces(
shank_recording,
time_range=(start_time, start_time + 0.1),
return_in_uV=True,
order_channel_by_depth=True,
ax=ax[shank_idx]
)
ax[shank_idx].set_title(f"Shank: {shank_idx}")
plt.show()
si_widgets.plot_traces(
split_recording[0],
time_range=(23.923, 23.926),
return_in_uV=True,
order_channel_by_depth=True,
mode="map"
)
ax = plt.gca()
ax.set_ylim(200, 300)
plt.show()
Much better!
Now we have our preprocessed data, and have visually confirmed the data quality is good, we are ready to perform spike sorting!
median vs. mean for common_average?return_in_uV=False. Does the data look much different?.save attributeget_traces() with and without return_in_uV as True and False to inspect the data type.get_traces()Now we are ready to move onto spike sorting. There are many possible spike sorters to choose from, today we will use Kilosort 4.
All we need to do is pass our preprocessed recording to the sorter of choice, along with any parameters. In our case, we will skip the filtering and common average referencing that Kilosort 4 can perform, as we have already run these steps in SpikeInterface (we cannot turn filtering off, so just set it to a very low cutoff).
Kilosort will then detect spikes and assign each spike to a unit (i.e. a putative neuron). It will output the time of each spike, what unit it belongs to.
It is often easiest to sort each shank separately. From now on, we will only sort the first shank of the recording.
In SpikeInterface, running a sorter is quite simple:
from pathlib import Path
import spikeinterface.sorters as si_sorters
sorter_folder = Path("kilosort4_output")
sorter_output_folder = sorter_folder / "sorter_output"
first_shank_recording = split_recording[0]
sorting = si_sorters.run_sorter(
"kilosort4",
first_shank_recording,
folder=sorter_folder,
remove_existing_folder=True,
do_CAR=False,
highpass_cutoff=0.1,
)
print(sorting)
0%| | 0/2 [00:00<?, ?it/s]
50%|βββββ | 1/2 [00:07<00:07, 7.14s/it]
100%|ββββββββββ| 2/2 [00:13<00:00, 6.56s/it]
100%|ββββββββββ| 2/2 [00:13<00:00, 6.65s/it]
0%| | 0/2 [00:00<?, ?it/s]
50%|βββββ | 1/2 [00:07<00:07, 7.52s/it]
100%|ββββββββββ| 2/2 [00:13<00:00, 6.76s/it]
100%|ββββββββββ| 2/2 [00:13<00:00, 6.88s/it]
0%| | 0/1 [00:00<?, ?it/s]
100%|ββββββββββ| 1/1 [00:00<00:00, 21.83it/s]
100%|ββββββββββ| 1/1 [00:00<00:00, 21.44it/s]
0%| | 0/2 [00:00<?, ?it/s]
50%|βββββ | 1/2 [00:00<00:00, 1.17it/s]
100%|ββββββββββ| 2/2 [00:01<00:00, 1.28it/s]
100%|ββββββββββ| 2/2 [00:01<00:00, 1.27it/s]
0%| | 0/1 [00:00<?, ?it/s]
100%|ββββββββββ| 1/1 [00:00<00:00, 45.36it/s]
100%|ββββββββββ| 1/1 [00:00<00:00, 44.11it/s]KiloSortSortingExtractor: 40 units - 1 segments - 30.0kHzThe sorting output is written to disk. The output contains all of the output files of the sorter, and its organisation is sorter-dependent. As we as spike times and unit assignments, it contains other outputs from the sorter e.g. representative waveforms for each units (βtemplatesβ).
In SpikeInterface, we can access the spike times and unit assignments in a sorter-independent manner, using the sorting object. For example, to get the spike times and unit assignments:
spikes = sorting.to_spike_vector()
spike_times = spikes["sample_index"] / sorting.get_sampling_frequency()
# show the first 5 spikes
print(
f"`spikes` is a structured number array with fields`: {spikes.dtype.names}\n"
)
print(
f"First 5 spike times (s): {spike_times[:5]}\n"
)
print(
f"First 5 unit assignments: {spikes['unit_index'][:5]}"
)`spikes` is a structured number array with fields`: ('sample_index', 'unit_index', 'segment_index')
First 5 spike times (s): [6.66666667e-05 6.66666667e-05 1.00000000e-04 1.00000000e-04
1.33333333e-04]
First 5 unit assignments: [ 3 11 7 19 39]TODO! Depending on timing and content of the other sections, we may want to reduce this down to a single command where we compute all of the extensions, and then go through them in more detail in the SpikeInterface GUI
Spike sorting is a complex task undertaken in the precense of significant noise. Therefore, many units will not cleanly represent a single neuron. As well as good units (i.. where each spike is from the same putative neuron) units may be noise (XXX) or multi-unit activity (MUA).
To determine whether units are good or not, we must look at the individual spikes that belong to a unit, and compute metrics. Cut then out, waveforms. In this section we will go over some key concepts at a very high level - we will develop this much further in the Postprocessing chapter.
SortingAnalzyerThe SortingAnalzyer is the key object for postprocessing in SpikeInterface. We can create the SortingAnalzyer by passing the sorting output and preprocessed recording:
from spikeinterface import create_sorting_analyzer
analyzer = create_sorting_analyzer(
sorting=sorting,
recording=first_shank_recording
)Using this sorting analyzer interface, we will compute all of the postprocessing results needed to assess the quality of our sorting.
Some sorters will output their own waveforms, templates, amplitudes ETC (for example, Kilosort has a templates.npy file containnig the templates it generated during sorting).
These are not loaded by SpikeInterface. Instead, the sorting analzyer will always recompute these features directly from the passed recording.
In our sorting, we can have many many detected spikes. This may be in the order of millios of spikes for an hour long recording.
In general, we will compute summary statistics over metrics computed from individual spikes. Therefore, including millions of spikes is computationally intensive while providing marginal gains. Instead, we can subselect spikes to include in postprocessing computations (typically, 3-500 spikes per unit is used a default).
analyzer.compute(
"random_spikes",
method="uniform",
max_spikes_per_unit=500,
)<spikeinterface.core.analyzer_extension_core.ComputeRandomSpikes at 0x7f5e1fa5abd0>To access the data assicated with a computed extension, we can use the get_extension function as below. We can see the indices (and compute the spike times) for the selected spikes:
spike_indices = analyzer.get_extension("random_spikes").data["random_spikes_indices"]
spike_times = spike_indices / recording.get_sampling_frequency()
print(spike_times)[0.00000000e+00 3.33333333e-05 6.66666667e-05 ... 5.10000000e-02
5.10333333e-02 5.10666667e-02]The action potential waveforms are important for assessing unit quality. For example, we expect that the waveforms of spikes which belong to a single unit are all similar.
To compute waveforms in SpikeInterface, we can use the waveforms extension. This will βchop outβ the spike waveform from the data, with a time window specified by the ms_before and ms_after (the spike peak) arguments. For example in the schematic below, we could cut out the blue spike in a window defined by the red blox:
analyzer.compute("waveforms", ms_before=1.0, ms_after=2.0)<spikeinterface.core.analyzer_extension_core.ComputeWaveforms at 0x7f5e1f6c5990>The waveforms array is three dimensional `(num waveforms, num samples, num channels)
waveforms = analyzer.get_extension("waveforms").data["waveforms"]
print(waveforms.shape)
# This is the waveform of the 19th spike (chosen as it looks nice)
plt.imshow(waveforms[18, :, :].T, aspect="auto")
plt.show()(1533, 90, 26)
The template is a representative wavecform from a putative neuron. Often, it is computed as the average waveform from all of the spikes that were assigned to a unit.
We can compute the extension and access the templates with:
analyzer.compute("templates")
templates = analyzer.get_extension("templates").data["average"]
print(
templates.shape
)(40, 90, 192)SpikeInterface also includes many useful post-processing plots in the widget smodule. For example we can plot a unit template as below:
# Plot the template waveform of the first unit
si_widgets.plot_unit_templates(analyzer, unit_ids=[0])
plt.show()
Finally, spike amplitudes are a useful measure for assessing unit quality. We expect that the amplitude of a spike is relatively stable througghout a recording. Changes in spike amplitudes may indicate problems with hte neuron health, drift, or other issues.
To compute the spike amplitudes, we can run:
analyzer.compute("spike_amplitudes")
amplitudes = analyzer.get_extension("spike_amplitudes").data["amplitudes"]
print(amplitudes)[ 0.76293945 -57.9834 10.681152 ... -87.73804 -76.293945
-61.798096 ]This gives a list of amplitudes, one for every spike. We will explore how visualisng spike amplitude across time is a useful measure for quality asessment later.
From this, we can compute quality metrics used to filter the units e.g. bombcell. We cover this more in the postprocessing chapter.
For now, see SpikeInterface GUI video?
https://www.youtube.com/watch?v=OGlQtxzip-Q