Bombcell integration

[Julie-Fabre]

This PR ports bombcell-style unit classification to SpikeInterface.

Template metrics

• Rewrote peak/trough detection with a new get_trough_and_peak_idx() function that uses scipy.signal.find_peaks(). Since SpikeInterface stores templates based on raw data rather than the heavily smoothed templates used in template matching, the waveforms can be noisy—so you can optionally apply Savitzky-Golay smoothing before detection. The function returns dicts for troughs, peaks before, and peaks after, each containing indices, values, prominences, and widths.

from spikeinterface.postprocessing import get_trough_and_peak_idx

troughs, peaks_before, peaks_after = get_trough_and_peak_idx(
    templates,
    sampling_frequency,
    smooth=True,
    min_thresh_detect_peaks_troughs=0.4,
)

• New metrics: peak_before_to_trough_ratio, peak_after_to_trough_ratio, waveform_baseline_flatness, peak_before_width, trough_width, main_peak_to_trough_ratio.

• Renamed peak_to_valley to peak_to_trough_duration.

analyzer.compute("template_metrics", metric_names=[
    "peak_before_to_trough_ratio",
    "waveform_baseline_flatness",
    "trough_width",
])

Quality metrics

• Added snr_bombcell—peak amplitude over baseline MAD.

analyzer.compute("quality_metrics", metric_names=["snr_bombcell"])

• amplitude_cutoff now has parameters for controlling the histogram fitting:

analyzer.compute("quality_metrics", metric_names=["amplitude_cutoff"], qm_params={
    "amplitude_cutoff": {
        "num_histogram_bins": 100,
        "histogram_smoothing_value": 3,
    }
})

Unit classification

• New in spikeinterface.curation:

import spikeinterface.comparison as sc

thresholds = sc.bombcell_get_default_thresholds()
unit_type, unit_type_string = sc.bombcell_classify_units(
    quality_metrics,
    thresholds=thresholds,
    classify_non_somatic=True,
)
summary = sc.get_classification_summary(unit_type, unit_type_string)

Units get classified as NOISE → MUA → GOOD based on successive threshold checks. Optional NON_SOMA category for non-somatic waveforms.

Plots

• Added plots for classification summaries, metric histograms with threshold lines, waveform overlays by category, and UpSet plots.

from spikeinterface.widgets import (
    plot_unit_classification,
    plot_classification_histograms,
    plot_waveform_overlay,
    plot_upset,
)

plot_unit_classification(analyzer, unit_type, unit_type_string)
plot_classification_histograms(quality_metrics, thresholds=thresholds)
plot_waveform_overlay(analyzer, unit_type, unit_type_string)
plot_upset(quality_metrics, unit_type, unit_type_string)

or a wrapper for all plots:

plots = plot_unit_classification_all(
    sorting_analyzer,
    unit_type,
    unit_type_string,
    quality_metrics=quality_metrics,  # optional, will try to get from analyzer
    thresholds=thresholds,            # optional, uses defaults
    split_non_somatic=False,
    include_upset=True,
)

TODO:

• Add documentation tutorial

[samuelgarcia]

Thanks Julie for this plot this is really cool. lets have a call to discuss no ?
I suspect that some difference may also be due to the fact that we measurign the negative peak only.
lets check this.

[Julie-Fabre]

Yes, absolutely I agree! I think that's what's responsible for all the cases where the SI metric < BC metric (blue guys) but maybe not for the cases where SI metric > BC metric (red guys)? Curious to hear what you think

Great idea, let's have a call. Sent you a message on slack :)