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, 374.24it/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, 16.86it/s]
compute_waveforms (workers: 2 processes): 100%|██████████| 10/10 [00:00<00:00, 56.03it/s]

noise_level (no parallelization):   0%|          | 0/20 [00:00<?, ?it/s]
noise_level (no parallelization): 100%|██████████| 20/20 [00:00<00:00, 255.72it/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'): 14.4, np.str_('1'): 14.7, np.str_('2'): 15.3, np.str_('3'): 14.4, np.str_('4'): 15.5, np.str_('5'): 15.7, np.str_('6'): 12.0, np.str_('7'): 16.4, np.str_('8'): 15.7, np.str_('9'): 13.8}
{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(25.204094112797932), np.str_('1'): np.float64(7.825129685241889), np.str_('2'): np.float64(14.764279132804175), np.str_('3'): np.float64(1.778879068792805), np.str_('4'): np.float64(9.557946185659846), np.str_('5'): np.float64(3.5978593259446545), np.str_('6'): np.float64(11.769432078201616), np.str_('7'): np.float64(37.0821118780606), np.str_('8'): np.float64(18.838303179271133), np.str_('9'): np.float64(3.1100953513372804)}

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         14.4  25.204094               NaN
1         14.7   7.825130               NaN
2         15.3  14.764279               NaN
3         14.4   1.778879               NaN
4         15.5   9.557946               NaN
5         15.7   3.597859               NaN
6         12.0  11.769432               NaN
7         16.4  37.082112               NaN
8         15.7  18.838303               NaN
9         13.8   3.110095               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, 213.00it/s]

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

calculate pc_metrics:   0%|          | 0/10 [00:00<?, ?it/s]
calculate pc_metrics: 100%|██████████| 10/10 [00:00<00:00, 348.39it/s]
   isolation_distance    d_prime  amplitude_cutoff        snr  firing_rate
0          467.768217  10.144100               NaN  25.204094         14.4
1           11.294896   0.919939               NaN   7.825130         14.7
2           87.113019   1.547497               NaN  14.764279         15.3
3            4.221376   2.225205               NaN   1.778879         14.4
4          142.627161   4.909966               NaN   9.557946         15.5
5            5.979716   1.815295               NaN   3.597859         15.7
6           46.064771   3.324997               NaN  11.769432         12.0
7         1989.268824  11.105151               NaN  37.082112         16.4
8          779.886814  12.200044               NaN  18.838303         15.7
9            8.068124   2.077843               NaN   3.110095         13.8