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, 481.48it/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):  70%|███████   | 7/10 [00:00<00:00, 61.27it/s]
compute_waveforms (workers: 2 processes): 100%|██████████| 10/10 [00:00<00:00, 80.16it/s]

noise_level (no parallelization):   0%|          | 0/20 [00:00<?, ?it/s]
noise_level (no parallelization): 100%|██████████| 20/20 [00:00<00:00, 291.63it/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'): 16.7, np.str_('1'): 14.6, np.str_('2'): 15.1, np.str_('3'): 16.0, np.str_('4'): 15.6, np.str_('5'): 15.8, np.str_('6'): 12.8, np.str_('7'): 15.2, np.str_('8'): 14.6, np.str_('9'): 13.5}
{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(4.463681918845905), np.str_('1'): np.float64(2.2697719474698657), np.str_('2'): np.float64(7.627989798451524), np.str_('3'): np.float64(0.8658098832612761), np.str_('4'): np.float64(23.89764692044111), np.str_('5'): np.float64(18.637072776422684), np.str_('6'): np.float64(46.803744160788874), np.str_('7'): np.float64(7.602967293973333), np.str_('8'): np.float64(1.553112139081377), np.str_('9'): np.float64(14.764471190284995)}

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         16.7   4.463682               NaN
1         14.6   2.269772               NaN
2         15.1   7.627990               NaN
3         16.0   0.865810               NaN
4         15.6  23.897647               NaN
5         15.8  18.637073               NaN
6         12.8  46.803744               NaN
7         15.2   7.602967               NaN
8         14.6   1.553112               NaN
9         13.5  14.764471               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, 183.27it/s]

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

calculate pc_metrics:   0%|          | 0/10 [00:00<?, ?it/s]
calculate pc_metrics: 100%|██████████| 10/10 [00:00<00:00, 524.43it/s]
   isolation_distance    d_prime  firing_rate        snr  amplitude_cutoff
0           11.036483   1.507127         16.7   4.463682               NaN
1            7.191853   1.384724         14.6   2.269772               NaN
2          108.827310   2.099603         15.1   7.627990               NaN
3            7.468967   1.715334         16.0   0.865810               NaN
4         1558.911586  18.952799         15.6  23.897647               NaN
5          337.247538   8.365758         15.8  18.637073               NaN
6         9761.434491  21.351398         12.8  46.803744               NaN
7           29.040292   2.069659         15.2   7.602967               NaN
8            3.338606   1.326974         14.6   1.553112               NaN
9           91.239511   7.220374         13.5  14.764471               NaN