Note
Go to the end to download the full example code.
Curation Tutorial
After spike sorting and computing quality metrics, you can automatically curate the spike sorting output using the quality metrics that you have calculated.
Import the modules and/or functions necessary from spikeinterface
import spikeinterface.core as si
Let’s generate a simulated dataset, and imagine that the ground-truth sorting is in fact the output of a sorter.
recording, sorting = si.generate_ground_truth_recording()
print(recording)
print(sorting)
GroundTruthRecording (InjectTemplatesRecording): 4 channels - 25.0kHz - 1 segments
250,000 samples - 10.00s - float32 dtype - 3.81 MiB
GroundTruthSorting (NumpySorting): 10 units - 1 segments - 25.0kHz
Create SortingAnalyzer
For this example, we will need a SortingAnalyzer and some extensions
to be computed first
analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")
analyzer.compute(["random_spikes", "waveforms", "templates", "noise_levels"])
analyzer.compute("principal_components", n_components=3, mode="by_channel_local")
print(analyzer)
estimate_sparsity (no parallelization): 0%| | 0/10 [00:00<?, ?it/s]
estimate_sparsity (no parallelization): 100%|██████████| 10/10 [00:00<00:00, 487.88it/s]
compute_waveforms (no parallelization): 0%| | 0/10 [00:00<?, ?it/s]
compute_waveforms (no parallelization): 100%|██████████| 10/10 [00:00<00:00, 348.71it/s]
noise_level (no parallelization): 0%| | 0/20 [00:00<?, ?it/s]
noise_level (no parallelization): 100%|██████████| 20/20 [00:00<00:00, 281.94it/s]
Fitting PCA: 0%| | 0/10 [00:00<?, ?it/s]
Fitting PCA: 100%|██████████| 10/10 [00:00<00:00, 180.54it/s]
Projecting waveforms: 0%| | 0/10 [00:00<?, ?it/s]
Projecting waveforms: 100%|██████████| 10/10 [00:00<00:00, 2623.74it/s]
SortingAnalyzer: 4 channels - 10 units - 1 segments - memory - sparse - has recording
Loaded 5 extensions: random_spikes, waveforms, templates, noise_levels, principal_components
Then we compute some quality metrics:
metrics_ext = analyzer.compute("quality_metrics", metric_names=["snr", "isi_violation", "nearest_neighbor"])
metrics = metrics_ext.get_data()
print(metrics)
snr isi_violations_ratio ... nn_hit_rate nn_miss_rate
0 9.550875 0.0 ... 0.802581 0.024819
1 6.559130 0.0 ... 0.811765 0.026853
2 16.254619 0.0 ... 0.863804 0.019124
3 38.071867 0.0 ... 0.929936 0.003779
4 27.614734 0.0 ... 0.896732 0.006087
5 5.504213 0.0 ... 0.785034 0.022367
6 6.504990 0.0 ... 0.786765 0.017752
7 0.973576 0.0 ... 0.750649 0.024801
8 10.940392 0.0 ... 0.812766 0.024282
9 12.656074 0.0 ... 0.829299 0.021076
[10 rows x 5 columns]
We can now threshold each quality metric and select units based on some rules.
The easiest and most intuitive way is to use boolean masking with a dataframe.
Then create a list of unit ids that we want to keep
keep_mask = (metrics["snr"] > 7.5) & (metrics["isi_violations_ratio"] < 0.2) & (metrics["nn_hit_rate"] > 0.80)
print(keep_mask)
keep_unit_ids = keep_mask[keep_mask].index.values
keep_unit_ids = [unit_id for unit_id in keep_unit_ids]
print(keep_unit_ids)
0 True
1 False
2 True
3 True
4 True
5 False
6 False
7 False
8 True
9 True
dtype: bool
['0', '2', '3', '4', '8', '9']
And now let’s create a sorting that contains only curated units and save it.
curated_sorting = sorting.select_units(keep_unit_ids)
print(curated_sorting)
curated_sorting.save(folder="curated_sorting", overwrite=True)
GroundTruthSorting (UnitsSelectionSorting): 6 units - 1 segments - 25.0kHz
We can also save the analyzer with only theses units
clean_analyzer = analyzer.select_units(unit_ids=keep_unit_ids, format="zarr", folder="clean_analyzer")
print(clean_analyzer)
SortingAnalyzer: 4 channels - 6 units - 1 segments - zarr - sparse - has recording
Loaded 6 extensions: random_spikes, waveforms, templates, noise_levels, principal_components, quality_metrics
Total running time of the script: (0 minutes 0.543 seconds)