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

[Julie-Fabre]

Thank you both so much @samuelgarcia and @alejoe91 !
Implementing the changes now and checking the SNR - that was a good point, I had implemented the 'bombcell' version to actually check and see how correlated they were (and hopefully drop the bombcell one and get a good threshold for the 'spikeInterface' one).
For the JSON file, I would like some option so users can load in a file if they want to. This is really nice for instance if people are sharing pipelines & parameters, if you want to reproduce a paper, etc. If not a JSON file, do you have other ideas for how we could implement that?

[Julie-Fabre]

Regarding the 'spikeInterface' vs 'bombcell' SNR:

• spikeInterace SNR: |peak(mean_template)| / MAD(recording_noise), is that right?
• Bombcell SNR: |peak(median_waveform)| / MAD(waveform_baseline). Uses median waveform and baseline noise from extracted spikes.

The do agree mostly:

• Here's a correlation plot:

• And some example units where they don't agree (top row = high spikeInterface SNR and low bombcell SNR; second row: low SI and high BC; third row: both high and fourth row: both low). I am plotting the median waveform +/- m.a.d. here.

Based on this very quick look at these examples, I would say:
(1) I think the spikeInterface way of choosing the extreme value doesn't work for non-somatic units (basically all the ones in the second row - blue ones). we should update that.
(2) for the top row (spikeInterface high, bombcell low), it's a mixed bag (some fine, others noisy) but I might tend to prefer the bombcell one (is does seem to seperate better what visually I would say are noisier units from less noisy ones), so I would like keep it at least for the bombcell pipeline.

I would really like to get your input on this! What do you think? And also what was the rationale for using the full recording to calculate noise?

[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 :)