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
import spikeinterface.extractors as se
from spikeinterface.postprocessing import compute_principal_components
from spikeinterface.qualitymetrics import (
    compute_snrs,
    compute_firing_rates,
    compute_isi_violations,
    calculate_pc_metrics,
    compute_quality_metrics,
)

First, let’s download a simulated dataset from the repo ‘https://gin.g-node.org/NeuralEnsemble/ephy_testing_data’

local_path = si.download_dataset(remote_path="mearec/mearec_test_10s.h5")
recording, sorting = se.read_mearec(local_path)
print(recording)
print(sorting)

MEArecRecordingExtractor: 32 channels - 32.0kHz - 1 segments - 320,000 samples - 10.00s
                          float32 dtype - 39.06 MiB
  file_path: /home/docs/spikeinterface_datasets/ephy_testing_data/mearec/mearec_test_10s.h5
MEArecSortingExtractor: 10 units - 1 segments - 32.0kHz
  file_path: /home/docs/spikeinterface_datasets/ephy_testing_data/mearec/mearec_test_10s.h5

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, 1121.23it/s]
SortingAnalyzer: 32 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):  50%|█████     | 5/10 [00:00<00:00, 41.27it/s]
compute_waveforms (workers: 2 processes): 100%|██████████| 10/10 [00:00<00:00, 67.03it/s]

noise_level (no parallelization):   0%|          | 0/20 [00:00<?, ?it/s]
noise_level (no parallelization):  30%|███       | 6/20 [00:00<00:00, 55.35it/s]
noise_level (no parallelization):  60%|██████    | 12/20 [00:00<00:00, 55.37it/s]
noise_level (no parallelization):  90%|█████████ | 18/20 [00:00<00:00, 55.19it/s]
noise_level (no parallelization): 100%|██████████| 20/20 [00:00<00:00, 55.27it/s]
SortingAnalyzer: 32 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)

{'#0': 5.3, '#1': 5.0, '#2': 4.3, '#3': 3.0, '#4': 4.8, '#5': 3.7, '#6': 5.1, '#7': 11.1, '#8': 19.5, '#9': 12.9}
{'#0': 0.0, '#1': 0.0, '#2': 0.0, '#3': 0.0, '#4': 0.0, '#5': 0.0, '#6': 0.0, '#7': 0.0, '#8': 0.0, '#9': 0.0}
{'#0': 23.74257172594366, '#1': 25.442748383832125, '#2': 13.721017396214569, '#3': 21.774053024200047, '#4': 7.440554835744134, '#5': 7.4455556575340385, '#6': 20.958460743516735, '#7': 7.350317068466804, '#8': 8.065364978512886, '#9': 8.986552883041492}

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          5.3  23.742572               NaN
#1          5.0  25.442748               NaN
#2          4.3  13.721017               NaN
#3          3.0  21.774053               NaN
#4          4.8   7.440555               NaN
#5          3.7   7.445556               NaN
#6          5.1  20.958461               NaN
#7         11.1   7.350317               NaN
#8         19.5   8.065365               NaN
#9         12.9   8.986553               NaN

Some metrics are based on the principal component scores, so the exwtension need to 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:  90%|█████████ | 9/10 [00:00<00:00, 69.38it/s]
Fitting PCA: 100%|██████████| 10/10 [00:00<00:00, 60.75it/s]

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

calculate pc_metrics:   0%|          | 0/10 [00:00<?, ?it/s]
calculate pc_metrics:  30%|███       | 3/10 [00:00<00:00, 18.01it/s]
calculate pc_metrics: 100%|██████████| 10/10 [00:00<00:00, 42.48it/s]
    isolation_distance    d_prime        snr  firing_rate  amplitude_cutoff
#0        1.116311e+17  29.314861  23.742572          5.3               NaN
#1        2.250397e+04  27.426984  25.442748          5.0               NaN
#2        1.260674e+17  24.759110  13.721017          4.3               NaN
#3        8.952986e+17  30.909950  21.774053          3.0               NaN
#4        1.674176e+17  28.448716   7.440555          4.8               NaN
#5        8.661268e+16  20.730191   7.445556          3.7               NaN
#6        2.142995e+18  37.228947  20.958461          5.1               NaN
#7        8.878391e+03  28.968743   7.350317         11.1               NaN
#8        5.827407e+03  20.999116   8.065365         19.5               NaN
#9        7.100391e+03  30.249492   8.986553         12.9               NaN