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,
)

[samuelgarcia]

Salut Julie,
I read this super quickly. This is super impressive what you did during the hackahton!
I was not aware that you also did the widgets stuff. Waou.

I will be back with more carefully reading.

But some main stuff:

• we avoid to push ipynb in the repo because in saturate the history so we use jupytext instead and push only the resulting generated rst,if the notebook is fast to generate (with simulate data) we also have the tutorial way to push doc through notebooks which is a py file run and generated by the documentaion build.
• I would prefer to not have json directly in the code to handle parameters. I think simple python file with the same contents. lets discuss more
• I would be courious to see the correlation between the basic snr and the one median based you did. I will try to make some plot on this.

[alejoe91]

@Julie-Fabre massive effort! Thanks!

I did a first round of reviewing and I'm happy to discuss some details and also work on it :)

[alejoe91]

I'm thinking it could be useful to add a deprecated_column_names, so we could automate backward compatibility :)

[Julie-Fabre]

That's a great idea, go for it! :)

[alejoe91]

I would be consistent and output everything in the same unit. For now we have been doing seconds for the durations. The bombcell curation could still accept thresholds in us and do the conversion on the fly.

Alternatively, we could add a unit field to the BaseMetric, to specify units for each column. I think I would go with this, but it requires an additional refactoring. @chrishalcrow what do you think?

[Julie-Fabre]

Converted to seconds, I think that makes sense. You don't have to inherit some of bombcell's early decisions 😂

[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?