Quality Metrics Tutorial

After spike sorting, you might want to validate the ‘goodness’ of the sorted units. This can be done using the qualitymetrics submodule, which computes several quality metrics of the sorted units.

import spikeinterface.core as si
from spikeinterface.qualitymetrics import (
    compute_snrs,
    compute_firing_rates,
    compute_isi_violations,
    compute_quality_metrics,
)

First, let’s generate a simulated recording and sorting

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 quality metrics we need first to create a SortingAnalyzer.

analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")
print(analyzer)

estimate_sparsity (no parallelization):   0%|          | 0/10 [00:00<?, ?it/s]
estimate_sparsity (no parallelization): 100%|██████████| 10/10 [00:00<00:00, 393.72it/s]
SortingAnalyzer: 4 channels - 10 units - 1 segments - memory - sparse - has recording
Loaded 0 extensions

Depending on which metrics we want to compute we will need first to compute some necessary extensions. (if not computed an error message will be raised)

analyzer.compute("random_spikes", method="uniform", max_spikes_per_unit=600, seed=2205)
analyzer.compute("waveforms", ms_before=1.3, ms_after=2.6, n_jobs=2)
analyzer.compute("templates", operators=["average", "median", "std"])
analyzer.compute("noise_levels")

print(analyzer)

compute_waveforms (workers: 2 processes):   0%|          | 0/10 [00:00<?, ?it/s]
compute_waveforms (workers: 2 processes):  20%|██        | 2/10 [00:00<00:00, 15.22it/s]
compute_waveforms (workers: 2 processes): 100%|██████████| 10/10 [00:00<00:00, 52.11it/s]

noise_level (no parallelization):   0%|          | 0/20 [00:00<?, ?it/s]
noise_level (no parallelization): 100%|██████████| 20/20 [00:00<00:00, 234.09it/s]
SortingAnalyzer: 4 channels - 10 units - 1 segments - memory - sparse - has recording
Loaded 4 extensions: random_spikes, waveforms, templates, noise_levels

The spikeinterface.qualitymetrics submodule has a set of functions that allow users to compute metrics in a compact and easy way. To compute a single metric, one can simply run one of the quality metric functions as shown below. Each function has a variety of adjustable parameters that can be tuned.

firing_rates = compute_firing_rates(analyzer)
print(firing_rates)
isi_violation_ratio, isi_violations_count = compute_isi_violations(analyzer)
print(isi_violation_ratio)
snrs = compute_snrs(analyzer)
print(snrs)

{np.str_('0'): 15.6, np.str_('1'): 17.0, np.str_('2'): 15.1, np.str_('3'): 13.3, np.str_('4'): 14.9, np.str_('5'): 13.9, np.str_('6'): 16.6, np.str_('7'): 13.0, np.str_('8'): 15.5, np.str_('9'): 15.3}
{np.str_('0'): np.float64(0.0), np.str_('1'): np.float64(0.0), np.str_('2'): np.float64(0.0), np.str_('3'): np.float64(0.0), np.str_('4'): np.float64(0.0), np.str_('5'): np.float64(0.0), np.str_('6'): np.float64(0.0), np.str_('7'): np.float64(0.0), np.str_('8'): np.float64(0.0), np.str_('9'): np.float64(0.0)}
{np.str_('0'): np.float64(19.842763186333524), np.str_('1'): np.float64(33.68840562009274), np.str_('2'): np.float64(10.035751713475044), np.str_('3'): np.float64(17.86917721264813), np.str_('4'): np.float64(2.6849621035396583), np.str_('5'): np.float64(38.87348761547684), np.str_('6'): np.float64(1.4315573794004424), np.str_('7'): np.float64(20.333775682383994), np.str_('8'): np.float64(14.231208226819527), np.str_('9'): np.float64(15.677888916359883)}

To compute more than one metric at once, we can use the compute_quality_metrics function and indicate which metrics we want to compute. This will return a pandas dataframe:

metrics = compute_quality_metrics(analyzer, metric_names=["firing_rate", "snr", "amplitude_cutoff"])
print(metrics)

   firing_rate        snr  amplitude_cutoff
0         15.6  19.842763               NaN
1         17.0  33.688406               NaN
2         15.1  10.035752               NaN
3         13.3  17.869177               NaN
4         14.9   2.684962               NaN
5         13.9  38.873488               NaN
6         16.6   1.431557               NaN
7         13.0  20.333776               NaN
8         15.5  14.231208               NaN
9         15.3  15.677889               NaN

Some metrics are based on the principal component scores, so the exwtension must be computed before. For instance:

analyzer.compute("principal_components", n_components=3, mode="by_channel_global", whiten=True)

metrics = compute_quality_metrics(
    analyzer,
    metric_names=[
        "isolation_distance",
        "d_prime",
    ],
)
print(metrics)

Fitting PCA:   0%|          | 0/10 [00:00<?, ?it/s]
Fitting PCA: 100%|██████████| 10/10 [00:00<00:00, 206.22it/s]

Projecting waveforms:   0%|          | 0/10 [00:00<?, ?it/s]
Projecting waveforms: 100%|██████████| 10/10 [00:00<00:00, 1230.36it/s]

calculate pc_metrics:   0%|          | 0/10 [00:00<?, ?it/s]
calculate pc_metrics: 100%|██████████| 10/10 [00:00<00:00, 336.81it/s]
   isolation_distance    d_prime  firing_rate  amplitude_cutoff        snr
0          402.497311   5.207557         15.6               NaN  19.842763
1         1410.143672  11.038832         17.0               NaN  33.688406
2           37.527535   0.792897         15.1               NaN  10.035752
3          473.263989   9.309678         13.3               NaN  17.869177
4           14.117576   2.207172         14.9               NaN   2.684962
5          126.037789   7.247826         13.9               NaN  38.873488
6           15.604535   2.521057         16.6               NaN   1.431557
7           97.521220   3.876010         13.0               NaN  20.333776
8           35.073803   1.551301         15.5               NaN  14.231208
9          156.284413   2.689595         15.3               NaN  15.677889