From 2fa5f49d41ebaa68912df4057c7131fc1d8ac223 Mon Sep 17 00:00:00 2001 From: chrishalcrow Date: Tue, 15 Jul 2025 07:09:21 -0400 Subject: [PATCH 1/7] non core modules --- doc/how_to/customize_a_plot.rst | 2 + doc/how_to/handle_drift.rst | 2 + doc/modules/benchmark.rst | 42 +++++------- doc/modules/comparison.rst | 27 ++++---- doc/modules/curation.rst | 19 +++--- doc/modules/extractors.rst | 21 ++++-- doc/modules/generation.rst | 2 +- doc/modules/motion_correction.rst | 105 ++++++++++++++++++------------ doc/modules/postprocessing.rst | 50 +++++++++----- doc/modules/preprocessing.rst | 35 +++++----- doc/modules/sorters.rst | 53 ++++++++------- doc/modules/sortingcomponents.rst | 69 +++++++++++--------- doc/modules/widgets.rst | 26 ++++---- 13 files changed, 250 insertions(+), 203 deletions(-) diff --git a/doc/how_to/customize_a_plot.rst b/doc/how_to/customize_a_plot.rst index 1fa9cca166..7ae9e3e79c 100644 --- a/doc/how_to/customize_a_plot.rst +++ b/doc/how_to/customize_a_plot.rst @@ -1,3 +1,5 @@ +.. _customize-a-plot: + Customize a plot ================ diff --git a/doc/how_to/handle_drift.rst b/doc/how_to/handle_drift.rst index 86c146fd02..1559accf3b 100755 --- a/doc/how_to/handle_drift.rst +++ b/doc/how_to/handle_drift.rst @@ -4,6 +4,8 @@ %load_ext autoreload %autoreload 2 +.. _handle-drift-with-spikeinterface: + Handle motion/drift with spikeinterface ======================================= diff --git a/doc/modules/benchmark.rst b/doc/modules/benchmark.rst index faf53be790..4f78cbb469 100644 --- a/doc/modules/benchmark.rst +++ b/doc/modules/benchmark.rst @@ -1,8 +1,7 @@ Benchmark module ================ -This module contains machinery to compare some sorters against ground truth in many multiple situtation. - +This module contains machinery to compare some sorters against ground truth in multiple situations. ..notes:: @@ -13,22 +12,21 @@ This module contains machinery to compare some sorters against ground truth in m This module also aims to benchmark sorting components (detection, clustering, motion, template matching) using the same base class :py:func:`~spikeinterface.benchmark.BenchmarkStudy()` but specialized to a targeted component. -By design, the main class handle the concept of "levels" : this allows to compare several complexities at the same time. +By design, the main class handles the concept of "levels" : this allows you to compare several complexities at the same time. For instance, compare kilosort4 vs kilsort2.5 (level 0) for different noises amplitudes (level 1) combined with -several motion vectors (leevel 2). +several motion vectors (level 2). **Example: compare many sorters : a ground truth study** We have a high level class to compare many sorters against ground truth: :py:func:`~spikeinterface.benchmark.SorterStudy()` - -A study is a systematic performance comparison of several ground truth recordings with several sorters or several cases +A study is a systematic performance comparison of several ground truth recordings with several sorters or several cases, like the different parameter sets. The study class proposes high-level tool functions to run many ground truth comparisons with many "cases" on many recordings and then collect and aggregate results in an easy way. -The all mechanism is based on an intrinsic organization into a "study_folder" with several subfolders: +The mechanism is based on an intrinsic organization into a "study_folder" with several subfolders: * datasets: contains ground truth datasets * sorters : contains outputs of sorters @@ -39,24 +37,20 @@ The all mechanism is based on an intrinsic organization into a "study_folder" wi .. code-block:: python - import matplotlib.pyplot as plt - import seaborn as sns - - import spikeinterface.extractors as se + import spikeinterface as si import spikeinterface.widgets as sw from spikeinterface.benchmark import SorterStudy - # generate 2 simulated datasets (could be also mearec files) - rec0, gt_sorting0 = generate_ground_truth_recording(num_channels=4, durations=[30.], seed=42) - rec1, gt_sorting1 = generate_ground_truth_recording(num_channels=4, durations=[30.], seed=91) + rec0, gt_sorting0 = si.generate_ground_truth_recording(num_channels=4, durations=[30.], seed=42) + rec1, gt_sorting1 = si.generate_ground_truth_recording(num_channels=4, durations=[30.], seed=91) datasets = { "toy0": (rec0, gt_sorting0), "toy1": (rec1, gt_sorting1), } - # define some "cases" here we want to test tridesclous2 on 2 datasets and spykingcircus2 on one dataset + # define some "cases". Here we want to test tridesclous2 on 2 datasets and spykingcircus2 on one dataset # so it is a two level study (sorter_name, dataset) # this could be more complicated like (sorter_name, dataset, params) cases = { @@ -74,20 +68,21 @@ The all mechanism is based on an intrinsic organization into a "study_folder" wi "label": "spykingcircus2 on tetrode0", "dataset": "toy0", "params": { - "sorter_name": "spykingcircus", + "sorter_name": "spykingcircus2", "docker_image": True }, }, } - # this initilizes a folder + # this initializes a folder + study_folder = "my_study_folder" study = SorterStudy.create(study_folder=study_folder, datasets=datasets, cases=cases, levels=["sorter_name", "dataset"]) - # This internally do run_sorter() for all cases in one function + # This internally does run_sorter() for all cases in one function study.run() - # Run the benchmark : this internanly do compare_sorter_to_ground_truth() for all cases + # Run the benchmark : this internally does compare_sorter_to_ground_truth() for all cases study.compute_results() # Collect comparisons one by one @@ -104,9 +99,9 @@ The all mechanism is based on an intrinsic organization into a "study_folder" wi m = comp.get_confusion_matrix() w_comp = sw.plot_agreement_matrix(sorting_comparison=comp) - # Collect synthetic dataframes and display + # Collect synthetic dataframes and display. # As shown previously, the performance is returned as a pandas dataframe. - # The spikeinterface.comparison.get_performance_by_unit() function, + # The spikeinterface.comparison.get_performance_by_unit() function # gathers all the outputs in the study folder and merges them into a single dataframe. # Same idea for spikeinterface.comparison.get_count_units() @@ -116,13 +111,10 @@ The all mechanism is based on an intrinsic organization into a "study_folder" wi # this is a dataframe unit_counts = study.get_count_units() - # Study also have several plotting methods for plotting the result + # Study also has several plotting methods for plotting the result study.plot_agreement_matrix() study.plot_unit_counts() study.plot_performances(mode="ordered") - study.plot_performances(mode="snr") - - Benchmark spike collisions diff --git a/doc/modules/comparison.rst b/doc/modules/comparison.rst index 957c974d16..abd1ae9c5a 100644 --- a/doc/modules/comparison.rst +++ b/doc/modules/comparison.rst @@ -1,7 +1,6 @@ Comparison module ================= - SpikeInterface has a :py:mod:`~spikeinterface.comparison` module, which contains functions and tools to compare spike trains and templates (useful for tracking units over multiple sessions). @@ -132,7 +131,7 @@ Given: 4. **Compute performances** With the list of matched units we can compute performance metrics. - Given : **tp** the number of true positive events, **fp** number of false + Given: **tp** the number of true positive events, **fp** number of false positive events, **fn** the number of false negative events, **num_gt** the number of events of the matched tested units, the following metrics are computed for each GT unit: @@ -142,7 +141,7 @@ Given: * false_discovery_rate = fp / (tp + fp) * miss_rate = fn / num_gt - The overall performances can be visualised with the **confusion matrix**, where + The overall performances can be visualized with the **confusion matrix**, where the last column contains the **FN** counts and the last row contains the **FP** counts. .. image:: ../images/spikecomparison_confusion.png @@ -215,15 +214,16 @@ An **over-merged** unit has a relatively high agreement (>= 0.2 by default) for .. code-block:: python - local_path = download_dataset(remote_path='mearec/mearec_test_10s.h5') - recording, sorting_true = read_mearec(local_path) + from spikeinterface.widgets import plot_agreement_matrix, plot_confusion_matrix + from spikeinterface.generation import generate_ground_truth_recording + from spikeinterface.sorters import run_sorter + recording, sorting_true = generate_ground_truth_recording() # run a sorter and compare to ground truth - sorting_HS = run_sorter(sorter_name='herdingspike', recording=recording) + sorting_HS = run_sorter(sorter_name='herdingspikes', recording=recording) cmp_gt_HS = sc.compare_sorter_to_ground_truth(sorting_true, sorting_HS, exhaustive_gt=True) - # To have an overview of the match we can use the ordered agreement matrix plot_agreement_matrix(cmp_gt_HS, ordered=True) @@ -232,7 +232,6 @@ An **over-merged** unit has a relatively high agreement (>= 0.2 by default) for # perf = cmp_gt_HS.get_performance() - # The confusion matrix is also a good summary of the score as it has # the same shape as an agreement matrix, but it contains an extra column for FN # and an extra row for FP @@ -271,20 +270,20 @@ The :py:func:`~spikeinterface.comparison.compare_two_sorters()` returns the comp .. code-block:: python import spikeinterface as si - import spikeinterface.extractors as se import spikeinterface.sorters as ss - import spikeinterface.comparisons as sc - import spikinterface.widgets as sw + import spikeinterface.comparison as scmp + import spikeinterface.widgets as sw # First, let's generate a simulated dataset recording, sorting = si.generate_ground_truth_recording() + # Then run two spike sorters and compare their outputs. sorting_HS = ss.run_sorter(sorter_name='herdingspikes', recording=recording) sorting_TDC = ss.run_sorter(sorter_name='tridesclous', recording=recording) # Run the comparison # Let's see how to inspect and access this matching. - cmp_HS_TDC = sc.compare_two_sorters( + cmp_HS_TDC = scmp.compare_two_sorters( sorting1=sorting_HS, sorting2=sorting_TDC, sorting1_name='HS', @@ -339,7 +338,7 @@ Comparison of multiple sorters uses the following procedure: sorting_TDC = ss.run_sorter(sorter_name='tridesclous', recording=recording) # Compare multiple spike sorter outputs - mcmp = sc.compare_multiple_sorters( + mcmp = scmp.compare_multiple_sorters( sorting_list=[sorting_MS4, sorting_HS, sorting_TDC], name_list=['MS4', 'HS', 'TDC'], verbose=True, @@ -354,7 +353,7 @@ Comparison of multiple sorters uses the following procedure: print(mcmp.comparisons[('MS4', 'TDC')].get_matching()) # The global multi comparison can be visualized with this graph - sw.plot_multicomp_graph(multi_comparison=mcmp) + sw.plot_multicomparison_graph(multi_comparison=mcmp) # Consensus-based method # diff --git a/doc/modules/curation.rst b/doc/modules/curation.rst index 36770086fe..d1e45c80af 100644 --- a/doc/modules/curation.rst +++ b/doc/modules/curation.rst @@ -12,7 +12,7 @@ Curation with the ``SortingAnalyzer`` ------------------------------------- The :py:class:`~spikeinterface.core.SortingAnalyzer`, as seen in previous modules, -is a powerful tool to posprocess the spike sorting output, as it can compute many +is a powerful tool to postprocess the spike sorting output, as it can compute many extensions and metrics to further characterize the spike sorting results. To facilitate the spike sorting workflow, the :py:class:`~spikeinterface.core.SortingAnalyzer` @@ -34,7 +34,7 @@ a subset of units from a spike sorting output and to perform some merges: remove_unit_ids = [1, 2] sorting_analyzer2 = sorting_analyzer.remove_units(remove_unit_ids=remove_unit_ids) - # merge some units + # merge units 4 and 5, and units 7, 8 and 12 merge_unit_groups = [[4, 5], [7, 8, 12]] sorting_analyzer3 = sorting_analyzer2.merge_units( merge_unit_groups=merge_unit_groups, @@ -98,7 +98,8 @@ The function can act both on a ``BaseSorting`` or a ``SortingAnalyzer`` object. clean_sorting = remove_redundant_units( sorting, duplicate_threshold=0.9, - remove_strategy="max_spikes" + remove_strategy="max_spikes", + align=False, ) # remove redundant units from SortingAnalyzer object @@ -106,14 +107,14 @@ The function can act both on a ``BaseSorting`` or a ``SortingAnalyzer`` object. clean_sorting = remove_redundant_units( sorting_analyzer, duplicate_threshold=0.9, - remove_strategy="min_shift" + remove_strategy="minimum_shift", ) # in order to have a SortingAnalyer with only the non-redundant units one must # select the designed units remembering to give format and folder if one wants # a persistent SortingAnalyzer. clean_sorting_analyzer = sorting_analyzer.select_units(clean_sorting.unit_ids) -We recommend using the ``SortingAnalyzer`` approach, since the ``min_shift`` strategy keeps +We recommend using the ``SortingAnalyzer`` approach, since the ``minimum_shift`` strategy keeps the unit (among the redundant ones), with a better template alignment. @@ -139,13 +140,13 @@ merges. It offers multiple "presets" and the flexibility to apply individual ste # merges is a list of unit pairs, with unit_ids to be merged. merge_unit_pairs = compute_merge_unit_groups( - analyzer=analyzer, + sorting_analyzer=analyzer, preset="similarity_correlograms", ) # with resolve_graph=True, merges_resolved is a list of merge groups, # which can contain more than two units merge_unit_groups = compute_merge_unit_groups( - analyzer=analyzer, + sorting_analyzer=analyzer, preset="similarity_correlograms", resolve_graph=True ) @@ -189,14 +190,14 @@ prone to wrong merges. To do so, you'll need to do the following: The extra keyword ``recursive`` specifies that for each presets/sequences of steps, merges are performed until no further merges are possible. The ``job_kwargs`` are the parameters for the parallelization. -**Be careful that the merges can not be reverted, so be sure to not erase your analyzer and create a new variable** +**Be careful: the merges can not be reverted, so be sure to not erase your analyzer and instead create a one** Manual curation --------------- While automatic curation tools can be very useful, manual curation is still widely used to -clean spike sorting outputs and it is sometoimes necessary to have a human in the loop. +clean spike sorting outputs and it is sometimes necessary to have a human in the loop. Curation format diff --git a/doc/modules/extractors.rst b/doc/modules/extractors.rst index ba08e45aca..44732914f5 100644 --- a/doc/modules/extractors.rst +++ b/doc/modules/extractors.rst @@ -28,13 +28,13 @@ Every format can be read with a simple function: .. code-block:: python - recording_oe = read_openephys(folder_path="open-ephys-folder") + recording_oe = se.read_openephys(folder_path="open-ephys-folder") - recording_spikeglx = read_spikeglx(folder_path="spikeglx-folder") + recording_spikeglx = se.read_spikeglx(folder_path="spikeglx-folder") - recording_blackrock = read_blackrock(folder_path="blackrock-folder") + recording_blackrock = se.read_blackrock(folder_path="blackrock-folder") - recording_mearec = read_mearec(file_path="mearec_file.h5") + recording_mearec = se.read_mearec(file_path="mearec_file.h5") Importantly, some formats directly handle the probe information: @@ -80,7 +80,7 @@ The actual reading will be done on demand using the :py:meth:`~spikeinterface.co .. code-block:: python # opening a 40GB SpikeGLX dataset is fast - recording_spikeglx = read_spikeglx(folder_path="spikeglx-folder") + recording_spikeglx = se.read_spikeglx(folder_path="spikeglx-folder") # this really does load the full 40GB into memory : not recommended!!!!! traces = recording_spikeglx.get_traces(start_frame=None, end_frame=None, return_scaled=False) @@ -110,7 +110,7 @@ You can install all extractors' dependencies with: .. code-block:: python - pip install spikeinterface[extractor] + pip install spikeinterface[extractors] Raw Data Formats @@ -119,6 +119,7 @@ Raw Data Formats For raw recording formats, we currently support: * **AlphaOmega** :py:func:`~spikeinterface.extractors.read_alphaomega()` +* **Axon** :py:func:`~spikeinterface.extractors.read_axon()` * **Axona** :py:func:`~spikeinterface.extractors.read_axona()` * **BlackRock** :py:func:`~spikeinterface.extractors.read_blackrock()` * **Binary** :py:func:`~spikeinterface.core.read_binary()` @@ -134,6 +135,7 @@ For raw recording formats, we currently support: * **Mountainsort MDA** :py:func:`~spikeinterface.extractors.read_mda_recording()` * **Neuralynx** :py:func:`~spikeinterface.extractors.read_neuralynx()` * **Neurodata Without Borders** :py:func:`~spikeinterface.extractors.read_nwb_recording()` +* **NeuroNexus** :py:func:`~spikeinterface.coextractorsre.read_neuronexus()` * **Neuroscope** :py:func:`~spikeinterface.coextractorsre.read_neuroscope_recording()` * **Neuroexplorer** :py:func:`~spikeinterface.extractors.read_neuroexplorer()` * **NIX** :py:func:`~spikeinterface.extractors.read_nix()` @@ -142,10 +144,12 @@ For raw recording formats, we currently support: * **Plexon** :py:func:`~spikeinterface.extractors.read_plexon()` * **Plexon 2** :py:func:`~spikeinterface.extractors.read_plexon2()` * **Shybrid** :py:func:`~spikeinterface.extractors.read_shybrid_recording()` +* **SpikeGadgets** :py:func:`~spikeinterface.extractors.read_spikegadgets()` * **SpikeGLX** :py:func:`~spikeinterface.extractors.read_spikeglx()` * **SpikeGLX IBL compressed** :py:func:`~spikeinterface.extractors.read_cbin_ibl()` * **SpikeGLX IBL stream** :py:func:`~spikeinterface.extractors.read_streaming_ibl()` * **Spike 2** :py:func:`~spikeinterface.extractors.read_spike2()` +* **Whitematter** :py:func:`~spikeinterface.extractors.read_whitematter()` * **TDT** :py:func:`~spikeinterface.extractors.read_tdt()` * **Zarr** :py:func:`~spikeinterface.core.read_zarr()` @@ -155,11 +159,13 @@ Sorted Data Formats For sorted data formats, we currently support: +* **ALF** :py:func:`~spikeinterface.extractors.read_alf_sorting()` * **BlackRock** :py:func:`~spikeinterface.extractors.read_blackrock_sorting()` * **Combinato** :py:func:`~spikeinterface.extractors.read_combinato()` * **Cell explorer** :py:func:`~spikeinterface.extractors.read_cellexplorer()` * **HerdingSpikes2** :py:func:`~spikeinterface.extractors.read_herdingspikes()` * **HDsort** :py:func:`~spikeinterface.extractors.read_hdsort()` +* **IBL sorter** :py:func:`~spikeinterface.extractors.read_ibl_sorting()` * **Kilosort1/2/2.5/3** :py:func:`~spikeinterface.extractors.read_kilosort()` * **Klusta** :py:func:`~spikeinterface.extractors.read_klusta()` * **MClust** :py:func:`~spikeinterface.extractors.read_mclust()` @@ -176,6 +182,7 @@ For sorted data formats, we currently support: * **Trideclous** :py:func:`~spikeinterface.extractors.read_tridesclous()` * **Wave Clus** :py:func:`~spikeinterface.extractors.read_waveclus()` * **YASS** :py:func:`~spikeinterface.extractors.read_yass()` +* **Zarr** :py:func:`~spikeinterface.extractors.read_zarr()` Dealing with Non-Supported File Formats @@ -183,4 +190,4 @@ Dealing with Non-Supported File Formats With recording and sorting objects, we hope that any user can access SpikeInterface regardless of the nature of their underlying file format. If you feel like a non-supported file format should be included in SpikeInterface as an -actual extractor, please open an issue. +actual extractor, please `open an issue `_. diff --git a/doc/modules/generation.rst b/doc/modules/generation.rst index 191cb57f30..44e5e2b7a1 100644 --- a/doc/modules/generation.rst +++ b/doc/modules/generation.rst @@ -1,7 +1,7 @@ Generation module ================= -The :py:mod:`spikeinterface.generation` provides functions to generate recordings containing spikes, +The :py:mod:`spikeinterface.generation` module provides functions to generate recordings containing spikes, which can be used as "ground-truth" for benchmarking spike sorting algorithms. There are several approaches to generating such recordings. diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index 076a560e31..be41b85f4b 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -4,21 +4,24 @@ Motion/drift correction ======================= +See a practical guide to motion correction in our How To guide: :ref:`handle-drift-with-spikeinterface`. + Overview -------- -Mechanical drift, often observed in recordings, is currently a major issue for spike sorting. This is especially striking -with the new generation of high-density devices used for in-vivo electrophysiology such as the neuropixel electrodes. -The first sorter that introduced motion/drift correction as a prepossessing step was Kilosort2.5 (see [Steinmetz2021]_ [SteinmetzDataset]_) +Mechanical drift, often observed in recordings when a probe is acutely inserted, is currently a major issue for spike sorting. This is especially striking +with the new generation of high-density devices used for in-vivo electrophysiology such as the NeuroPixels probes. +The first sorter that introduced motion/drift correction as a prepossessing step was Kilosort2.5 (see [Steinmetz2021]_ [SteinmetzDataset]_ [Pachitariu2023]_) -Long story short, the main idea is the same as the one used for non-rigid image registration, for example with calcium +The first algorithm used the same ideas as those used for non-rigid image registration, for example with calcium imaging. However, because with extracellular recording we do not have a proper image to use as a reference, the main idea of the algorithm is create an "image" via the activity profile of the cells during a given time window. Assuming this activity profile should be kept constant over time, the motion can be estimated, by blocks, along the probe's insertion axis (i.e. depth) so that we can interpolate the traces to compensate for this estimated motion. -Users with a need to handle drift were currently forced to stick to the use of Kilosort2.5 or pyKilosort (see [Pachitariu2023]_). Recently, the Paninski -group from Columbia University introduced a possibly more accurate method to estimate the drift (see [Varol2021]_ -and [Windolf2023]_), but this new method has not yet been integrated into any sorter. + +There are now several algorithms which try to correct for drift as a preprocessing step: the Paninski +group from Columbia University introduced DREDGE (see [Varol2021]_ +and [Windolf2023]_), and the Jazayeri lab introduced MEDiCINe ([Watters]_). Because motion registration is a hard topic, with numerous hypotheses and/or implementations details that might have a large impact on the spike sorting performances (see [Garcia2023]_), in SpikeInterface, we developed a full motion estimation @@ -26,7 +29,7 @@ and interpolation framework to make all these methods accessible in one place. T **the drift correction can be applied to a recording as a preprocessing step, and then used for any sorter!** In short, the motion correction is decoupled from the sorter itself. -This gives the user an incredible flexibility to check/test and correct the drift before the sorting process. +This gives the user flexibility to check/test and correct the drift before the sorting process. Here is an overview of motion correction as part of preprocessing a recording: @@ -59,16 +62,20 @@ One challenging task for motion correction is to determine the parameters. The high level :py:func:`~spikeinterface.preprocessing.correct_motion()` proposes the concept of a **"preset"** that already has predefined parameters, in order to achieve a calibrated behavior. -We currently have 3 presets: +We currently have presets: - * **"nonrigid_accurate"**: It consists of *monopolar triangulation + decentralized + inverse distance weighted* - This is the slowest combination but maybe the most accurate. The main bottleneck of this preset is the monopolar + * **"dredge"**: The official implementation of DREDGE, used in [Windolf2023]_. + * **"dredge_fast"**: A faster implementation of DREDGE, which should give similar results. + * **"nonrigid_accurate"**: A precursor to DREDGE. This consists of *monopolar triangulation + decentralized + inverse distance weighted* + It is the slowest combination, but maybe the most accurate. The main bottleneck of this preset is the monopolar triangulation for the estimation of the peaks positions. To speed it up, one could think about subsampling the space of all the detected peaks. Introduced by the Paninski group ([Varol2021]_, [Windolf2023]_) - * **"rigid_fast"**: a fast, but not very accurate method. It uses *center of mass + decentralized + inverse distance weighted* + * **"medicine"**: A wrapped version of MEDiCINe, [Watters]_. + * **"nonrigid_fast_and_accurate"**: A mixture of Kilosort and DREDGE ideas. Roughly: grid_convolution + decentralized motion estimation. + * **"rigid_fast"**: A fast, but not very accurate method. It uses *center of mass + decentralized + inverse distance weighted* To be used as check and/or control on a recording to check the presence of drift. Note that, in this case the drift is considered as "rigid" over the electrode. - * **"kilosort_like"**: It consists of *grid convolution + iterative_template + kriging*, to mimic what is done in Kilosort (see [Pachitariu2023]_). + * **"kilosort_like"**: This consists of *grid convolution + iterative_template + kriging*, to mimic what is done in Kilosort (see [Pachitariu2023]_). Note that this is not exactly 100% what Kilosort is doing, because the peak detection is done with a template matching in Kilosort, while in SpikeInterface we use a threshold-based method. However, this "preset" gives similar results to Kilosort2.5. @@ -104,21 +111,33 @@ Optionally any parameter from the preset can be overwritten: .. code-block:: python - rec_corrected = correct_motion(recording=rec, preset="nonrigid_accurate", - detect_kwargs=dict( - detect_threshold=10.), - estimate_motion_kwargs=dict( - histogram_depth_smooth_um=8., - time_horizon_s=120., - ), - correct_motion_kwargs=dict( - spatial_interpolation_method="kriging", - ) - ) + rec_corrected = correct_motion( + recording=rec, preset="nonrigid_accurate", + detect_kwargs=dict( + detect_threshold=10. + ), + estimate_motion_kwargs=dict( + histogram_depth_smooth_um=8., + time_horizon_s=120., + ), + correct_motion_kwargs=dict( + spatial_interpolation_method="kriging", + ) + ) + +Importantly, all the result and intermediate computations can returned to a motion object, for further loading, +verification and visualisation. -Importantly, all the result and intermediate computations can be saved into a folder for further loading -and verification. The folder will contain the motion vector itself of course but also detected peaks, peak location, and more. +.. code-block:: python + + motion_folder = '/somewhere/to/save/the/motion' + rec_corrected, motion = correct_motion(recording=rec, preset="nonrigid_accurate", output_motion=True) + from spikeinterface.widgets import plot_motion + plot_motion(motion) + +Alternatively, you can save the motion (and related motion info) in a folder. The folder will contain +the motion vector itself, as well as detected peaks, peak location, and more. .. code-block:: python @@ -129,7 +148,6 @@ and verification. The folder will contain the motion vector itself of course but motion_info = load_motion_info(motion_folder) - Low-level API ------------- @@ -144,7 +162,6 @@ working properly for your particular case. The high-level :py:func:`~spikeinterface.preprocessing.correct_motion()` is internally equivalent to this: - .. code-block:: python # each import is needed @@ -162,28 +179,30 @@ The high-level :py:func:`~spikeinterface.preprocessing.correct_motion()` is inte max_distance_um=150.0, **job_kwargs) # Step 2: motion inference - motion = estimate_motion(recording=rec, - peaks=peaks, - peak_locations=peak_locations, - method="decentralized", - direction="y", - bin_duration_s=2.0, - bin_um=5.0, - win_step_um=50.0, - win_sigma_um=150.0) + motion = estimate_motion( + recording=rec, + peaks=peaks, + peak_locations=peak_locations, + method="decentralized", + direction="y", + bin_um=5.0, + ) # Step 3: motion interpolation # this step is lazy - rec_corrected = interpolate_motion(recording=rec, motion=motion, - border_mode="remove_channels", - spatial_interpolation_method="kriging", - sigma_um=30.) + rec_corrected = interpolate_motion( + recording=rec, + motion=motion, + border_mode="remove_channels", + spatial_interpolation_method="kriging", + sigma_um=30. + ) Preprocessing details --------------------- -The function :py:func:`~spikeinterface.preprocessing.correct_motion()` requires an already preprocessed recording. +The function :py:func:`~spikeinterface.preprocessing.correct_motion()` requires a preprocessed recording. It is important to keep in mind that the preprocessing can have a strong impact on the motion estimation. @@ -236,3 +255,5 @@ References .. [Pachitariu2023] `Solving the spike sorting problem with Kilosort `_ .. [Garcia2023] `A modular approach to handle in-vivo drift correction for high-density extracellular recordings `_ + +.. [Watters] `MEDiCINe: Motion Correction for Neural Electrophysiology Recordings. 2025. `_ diff --git a/doc/modules/postprocessing.rst b/doc/modules/postprocessing.rst index a312633a08..c859e6d49f 100644 --- a/doc/modules/postprocessing.rst +++ b/doc/modules/postprocessing.rst @@ -1,11 +1,29 @@ Postprocessing module ===================== +SpikeInterface's postprocessing is centered around the :py:class:`~spikeinterface.core.SortingAnalyzer` +object. This combines recording and sorting information, and can be used to compute lots of information +you might be interested in: waveforms, templates, correlograms, spike amplitudes and more. We hope this +object can also serve as a basis for other tools, and is already used by SpikeInterface-GUI, SortingView +and NeuroConv. + +Learn how to make a SortingAnalyzer in `our tutorial +`_ + .. _extensions: After spike sorting, we can use the :py:mod:`~spikeinterface.postprocessing` module to further post-process -the spike sorting output. Most of the post-processing functions require a -:py:class:`~spikeinterface.core.SortingAnalyzer` as input. +the spike sorting output. First, make a :py:class:`~spikeinterface.core.SortingAnalyzer` object. We can +then compute "extensions" of the sorting analyzer. If you'd like to play with the extensions on a +simulated analyzer, we can create a simulated one as follows: + +.. code-block:: python + + import spikeinterface.core as si + + recording, sorting = si.generate_ground_truth_recording() + sorting_analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording) + .. _waveform_extensions: @@ -26,8 +44,6 @@ To check what extensions have already been calculated for a :code:`SortingAnalyz .. code-block:: python - import spikeinterface as si - available_extension_names = sorting_analyzer.get_loaded_extension_names() print(available_extension_names) @@ -43,12 +59,12 @@ To load the extension object you can run: ext = sorting_analyzer.get_extension("spike_amplitudes") ext_data = ext.get_data() -Here :code:`ext` is the extension object (in this case the :code:`SpikeAmplitudeCalculator`), and :code:`ext_data` will +Here :code:`ext` is the extension object (in this case the :code:`ComputeSpikeAmplitudes`), and :code:`ext_data` will contain the actual amplitude data. Note that different extensions might have different ways to return the extension. You can use :code:`ext.get_data?` for documentation. -To check what extensions spikeinterface can calculate, you can use the :code:`get_computable_extensions` method. +To check what extensions spikeinterface can calculate, you can use the :code:`get_computable_extensions` method: .. code-block:: python @@ -103,7 +119,7 @@ You can also compute several extensions at the same time by passing a list: .. code-block:: python - sorting_analyzer.compute(["principal_components", "templates"]) + sorting_analyzer.compute(["random_spikes", "waveforms", "principal_components"]) You might want to change the parameters. Two parameters of principal_components are :code:`n_components` and :code:`mode`. We can choose these as follows: @@ -118,8 +134,9 @@ than one thing: .. code-block:: python compute_dict = { + 'random_spikes': {}, 'principal_components': {'n_components': 3, 'mode': 'by_channel_local'}, - 'templates': {'operators': ["average"]} + 'templates': {'operators': ["average"]}, } sorting_analyzer.compute(compute_dict) @@ -146,7 +163,6 @@ Extensions are generally saved in two ways, suitable for two workflows: :code:`sorting_analyzer.compute('waveforms', save=False)`). - **NOTE**: We recommend choosing a workflow and sticking with it. Either keep everything on disk or keep everything in memory until you'd like to save. A mixture can lead to unexpected behavior. For example, consider the following code @@ -163,7 +179,7 @@ you'd like to save. A mixture can lead to unexpected behavior. For example, cons Here the random_spikes extension is **not** saved. The :code:`sorting_analyzer` is **still** saved in memory. The :code:`save_as` method only made a snapshot of the sorting analyzer which is saved in a folder. Hence :code:`compute` doesn't know about the folder -and doesn't save anything. If we wanted to save the extension we should have started with a non-memory sorting analyzer: +and doesn't save anything. If we wanted to save the extensions as we compute them, we should have started with a "binary_folder" sorting analyzer: .. code:: @@ -188,10 +204,11 @@ As an extension, this expects the :code:`Recording` as input and the computed va .. code-block:: python + # compute the noise from a recording noise = compute_noise_level(recording=recording) - - + # compute the noise from an analyzer + noise = sorting_analyzer.compute("noise_levels") principal_components @@ -219,7 +236,6 @@ For more information, see :py:func:`~spikeinterface.postprocessing.compute_princ template_similarity ^^^^^^^^^^^^^^^^^^^ - This extension computes the similarity of the templates to each other. This information could be used for automatic merging. Currently, the only available similarity method is the cosine similarity, which is the angle between the high-dimensional flattened template arrays. Note that cosine similarity does not take into account amplitude differences @@ -238,7 +254,7 @@ spike_amplitudes ^^^^^^^^^^^^^^^^ This extension computes the amplitude of each spike as the value of the traces on the extremum channel at the times of -each spike. +each spike. The extremum channel is compute from the templates. **NOTE:** computing spike amplitudes is highly recommended before calculating amplitude-based quality metrics, such as :ref:`amp_cutoff` and :ref:`amp_median`. @@ -267,7 +283,7 @@ with center of mass (:code:`method="center_of_mass"` - fast, but less accurate), input="spike_locations", ms_before=0.5, ms_after=0.5, - spike_retriever_kwargs=dict( + spike_retriver_kwargs=dict( channel_from_template=True, radius_um=50, peak_sign="neg" @@ -387,8 +403,8 @@ This extension computes the histograms of inter-spike-intervals. The computed ou .. code-block:: python - isi = sorting_analyer.compute( - input="isi_histograms" + isi = sorting_analyzer.compute( + input="isi_histograms", window_ms=50.0, bin_ms=1.0, method="auto" diff --git a/doc/modules/preprocessing.rst b/doc/modules/preprocessing.rst index 085b55c083..1769ab7e1c 100644 --- a/doc/modules/preprocessing.rst +++ b/doc/modules/preprocessing.rst @@ -32,15 +32,15 @@ These two preprocessors will not compute anything at instantiation, but the comp traces = recording_cmr.get_traces(start_frame=100_000, end_frame=200_000) -Some internal sorters (see :ref:`modules/sorters:Internal Sorters`) can work directly on these preprocessed objects so there is no need to +Most sorters (see :ref:`modules/sorters:Internal Sorters`) can work directly on these preprocessed objects so there is no need to save the object: .. code-block:: python - # here the spykingcircus2 sorter engine directly uses the lazy "recording_cmr" object + # here the spykingcircus2 sorter directly uses the lazy "recording_cmr" object sorting = run_sorter(sorter='spykingcircus2', recording=recording_cmr, sorter_name='spykingcircus2') -Most of the external sorters, however, will need a binary file as input, so we can optionally save the processed +Some sorters, however, will need a binary file as input, so we can optionally save the processed recording with the efficient SpikeInterface :code:`save()` function: .. code-block:: python @@ -51,7 +51,7 @@ In this case, the :code:`save()` function will process in parallel our original CMR, and save it to a binary file in the "/path/to/preprocessed" folder. The :code:`recording_saved` is yet another :code:`~spikeinterface.core.BaseRecording` which maps directly to the newly created binary file, for very quick access. -**NOTE:** all sorters will automatically perform the saving operation internally. +**NOTE:** some sorters will automatically perform the saving operation internally. The Preprocessing Pipeline -------------------------- @@ -64,7 +64,7 @@ bad channels step. We first make the appropriate dictionary .. code-block:: python - from spikeinterface.preprocessing import apply_pipeline, PreprocessingPipeline + from spikeinterface.preprocessing import apply_preprocessing_pipeline, PreprocessingPipeline preprocessing_dict = { 'highpass_filter': {'freq_min': 250}, @@ -72,11 +72,11 @@ bad channels step. We first make the appropriate dictionary 'detect_and_remove_bad_channels': {}, } -We can then pass this dictionary to the :code:`apply_pipeline` function to make a preprocessed recording +We can then pass this dictionary to the :code:`apply_preprocessing_pipeline` function to make a preprocessed recording .. code-block:: python - preprocessed_recording = apply_pipeline(recording, preprocessing_dict) + preprocessed_recording = apply_preprocessing_pipeline(recording, preprocessing_dict) Alternatively, we can construct a :code:`PreprocessingPipeline`, allowing us to investigate the pipeline before using it. @@ -139,7 +139,7 @@ Once we have the pipeline, we can apply it to a recording in the same way as app .. code-block:: python - preprocessed_recording_again = apply_pipeline(recording, preprocessing_pipeline) + preprocessed_recording_again = apply_preprocessing_pipeline(recording, preprocessing_pipeline) To share the pipeline you have made with another lab, you can simply share the dictionary. The dictionary can also be obtained from the pipeline object directly: @@ -211,7 +211,7 @@ Important aspects of filtering functions: common_reference() ^^^^^^^^^^^^^^^^^^ -A very common operation to remove the noise is to re-reference traces. +A very common operation to remove noise is to re-reference traces. This is implemented with the :code:`common_reference()` function. There are various options when combining :code:`operator` and :code:`reference` arguments: @@ -232,14 +232,14 @@ In fact, there is a small delay (less that a sampling period) in between channel For instance this is the case for Neuropixels devices. Applying :code:`common_reference()` on this data does not correctly remove artifacts, since we first need to compensate -for these small delays! This is exactly what :code:`phase_shift()` does. +for these small delays. This is exactly what :code:`phase_shift()` does. This function relies on an internal property of the recording called :code:`inter_sample_shift`. For Neuropixels recordings (read with the :py:func:`~spikeinterface.extractors.read_spikeglx` or the :py:func:`~spikeinterface.extractors.read_openephys` functions), the :code:`inter_sample_shift` is automatically loaded from the metadata and set. -Calling :code:`phase_shift()` alone has almost no effect, but combined with :code:`common_reference()` it makes a real +Calling :code:`phase_shift()` alone has almost no effect, but combined with :code:`common_reference()` it can make a real difference on artifact removal. @@ -326,7 +326,8 @@ The :code:`detect_bad_channels()` can be used to detect bad channels with severa approach to detect bad channels with abnormally high power and the :code:`coherence+psd` method (introduced by [IBL_spikesorting]_), which detects bad channels looking at both coherence with other channels and PSD power in the high-frequency range. -Note: The :code:`coherence+psd` method must be run on individual probes/shanks separately since it uses the coherence of the signal across the depth of the probe. See `Processing a Recording by Channel Group `_ for more information. +Note: The :code:`coherence+psd` method must be run on individual probes/shanks separately since it uses the coherence of the signal across the depth of the probe. +See `Processing a Recording by Channel Group `_ for more information. The function returns both the :code:`bad_channel_ids` and :code:`channel_labels`, which can be :code:`good`, :code:`noise`, :code:`dead`, or :code:`out` (outside of the brain). Note that the :code:`dead` and :code:`out` are only available with the :code:`coherence+psd` method. @@ -410,7 +411,7 @@ is subtracted, and the traces are finally cast to :code:`int16`: zero_channel_pad() ^^^^^^^^^^^^^^^^^^ -Pads a recording with extra channels that containing only zeros. This step can be useful when a certain shape is +Pads a recording with extra channels that contain only zeros. This step can be useful when a certain shape is required. .. code-block:: python @@ -449,7 +450,7 @@ explanation. deepinterpolation() (experimental) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The step (experimental) applies the inference step of a DeepInterpolation denoiser model [DeepInterpolation]_. +This step (experimental) applies the inference step of a DeepInterpolation denoiser model [DeepInterpolation]_. * :py:func:`~spikeinterface.preprocessing.deepinterpolation()` @@ -462,13 +463,13 @@ How to implement "IBL destriping" or "SpikeGLX CatGT" in SpikeInterface SpikeGLX has a built-in function called `CatGT `_ to apply some preprocessing on the traces to remove noise and artifacts. + IBL also has a standardized pipeline for preprocessed traces a bit similar to CatGT which is called "destriping" [IBL_spikesorting]_. -In both these cases, the traces are entirely read, processed and written back to a file. +In both these cases, the entire traces are read, processed and written back to a file. SpikeInterface can reproduce similar results without the need to write back to a file by building a *lazy* preprocessing chain. Optionally, the result can still be written to a binary (or a zarr) file. - Here is a recipe to mimic the **IBL destriping**: .. code-block:: python @@ -483,7 +484,6 @@ Here is a recipe to mimic the **IBL destriping**: rec.save(folder='clean_traces', n_jobs=10, chunk_duration='1s', progres_bar=True) - Here is a recipe to mimic the **SpikeGLX CatGT**: .. code-block:: python @@ -504,7 +504,6 @@ Of course, these pipelines can be enhanced and customized using other available Preprocessing on Snippets ------------------------- - Some preprocessing steps are available also for :py:class:`~spikeinterface.core.BaseSnippets` objects: align_snippets() diff --git a/doc/modules/sorters.rst b/doc/modules/sorters.rst index 393a357cef..b15daa76b3 100644 --- a/doc/modules/sorters.rst +++ b/doc/modules/sorters.rst @@ -4,7 +4,7 @@ Sorters module The :py:mod:`spikeinterface.sorters` module is where spike sorting happens! On one hand, SpikeInterface provides wrapper classes to many commonly used spike sorters like -Kilosort, Spyking-circus, etc. (see :ref:`compatible-sorters`). All these sorter classes inherit +Kilosort, Mountainsort, etc. (see :ref:`compatible-sorters`). All these sorter classes inherit from the :py:class:`~spikeinterface.sorters.BaseSorter` class, which provides the common tools to run spike sorters. @@ -56,17 +56,18 @@ to easily run spike sorters: # run Tridesclous sorting_TDC = run_sorter(sorter_name="tridesclous", recording=recording, folder="/folder_TDC") - # run Kilosort2.5 - sorting_KS2_5 = run_sorter(sorter_name="kilosort2_5", recording=recording, folder="/folder_KS2_5") + # run Kilosort4 + sorting_KS2_5 = run_sorter(sorter_name="kilosort4", recording=recording, folder="/folder_KS4") # run IronClust sorting_IC = run_sorter(sorter_name="ironclust", recording=recording, folder="/folder_IC") # run pyKilosort sorting_pyKS = run_sorter(sorter_name="pykilosort", recording=recording, folder="/folder_pyKS") - # run SpykingCircus - sorting_SC = run_sorter(sorter_name="spykingcircus", recording=recording, folder="/folder_SC") + # run SpykingCircus2 + sorting_SC2 = run_sorter(sorter_name="spykingcircus2", recording=recording, folder="/folder_SC2") -Then the output, which is a :py:class:`~spikeinterface.core.BaseSorting` object, can be easily +Whatever the sorter would usually save is saved in the ``folder`` you specify. In addition, each function +return a :py:class:`~spikeinterface.core.BaseSorting` object, which can be easily saved or directly post-processed: .. code-block:: python @@ -76,7 +77,7 @@ saved or directly post-processed: The :py:func:`~spikeinterface.sorters.run_sorter` function has several options: - * to remove or not the sorter working folder (:code:`output_folder/sorter_output`) + * to remove or not the sorter working folder (:code:`folder/sorter_output`) with: :code:`remove_existing_folder=True/False` (this can save lot of space because some sorters need data duplication!) * to control their verbosity: :code:`verbose=False/True` @@ -98,6 +99,8 @@ Parameters from all sorters can be retrieved with these functions: .. code-block:: python + from spikeinterface.sorters import get_default_sorter_params, get_sorter_params_description + params = get_default_sorter_params(sorter_name_or_class='spykingcircus') print("Parameters:\n", params) @@ -147,13 +150,13 @@ The containerized approach has several advantages: * Installation is much easier. * Different spike sorters with conflicting dependencies can be easily run side-by-side. * The results of the analysis are more reproducible and not dependant on the operating system -* MATLAB-based sorters can be run **without a MATLAB licence**. +* MATLAB-based sorters can be run **without a MATLAB license**. The containers can be run in Docker or Singularity, so having Docker or Singularity installed is a prerequisite. -Running spike sorting in a Docker container just requires: +Running spike sorting in a Docker container requires: 1) have docker installed 2) have docker Python SDK installed (:code:`pip install docker`) @@ -183,7 +186,9 @@ The following code creates a test recording and runs a containerized spike sorte .. code-block:: python - test_recording, _ = toy_example( + import spikeinterface.extractors as se + + test_recording, _ = se.toy_example( duration=30, seed=0, num_channels=64, @@ -273,7 +278,7 @@ Then you can build and tag the docker image with: docker build -t my-user/ks3-with-spikeinterface-test:0.1.0 . -And use the custom image whith the :code:`run_sorter` function: +And use the custom image with the :code:`run_sorter` function: .. code-block:: python @@ -296,6 +301,8 @@ an :code:`engine` that supports parallel processing (such as :code:`joblib` or : .. code-block:: python + from spikeinterface.sorters import run_sorter_jobs + # here we run 2 sorters on 2 different recordings = 4 jobs recording = ... another_recording = ... @@ -395,11 +402,12 @@ In this example, we create a 16-channel recording with 4 tetrodes: sorting = run_sorter(sorter_name='kilosort2', recording=recording, folder=f"folder_KS2_group{group}") sortings[group] = sorting +Read more about preprocessing and sorting by group in our How To, :ref:`recording-by-channel-group`. Handling multi-segment recordings --------------------------------- -In several experiments, several acquisitions are performed in sequence, for example a +In some experiments, several acquisitions are performed in sequence, for example a baseline/intervention. In these cases, since the underlying spiking activity can be assumed to be the same (or at least very similar), the recordings can be concatenated. This example shows how to concatenate the recordings before spike sorting and how to split the sorted output based @@ -413,20 +421,18 @@ do not handle multi-segment, and in that case we will use the .. code-block:: python + import spikeinterface.core as si # Let's create 4 recordings recordings_list = [] for i in range(4): - rec, _ = si.toy_example(duration=10., num_channels=4, seed=0, num_segments=1) + rec, _ = se.toy_example(duration=10., num_channels=4, seed=0, num_segments=1) recordings_list.append(rec) # Case 1: the sorter handles multi-segment objects multirecording = si.append_recordings(recordings_list) - # let's set a probe - multirecording = multirecording.set_probe(recording_single.get_probe()) - print(multirecording) # multirecording has 4 segments of 10s each # run tridesclous in multi-segment mode @@ -436,9 +442,6 @@ do not handle multi-segment, and in that case we will use the # Case 2: the sorter DOES NOT handle multi-segment objects # The `concatenate_recordings()` mimics a mono-segment object that concatenates all segments multirecording = si.concatenate_recordings(recordings_list) - # let's set a probe - multirecording = multirecording.set_probe(recording_single.get_probe()) - print(multirecording) # multirecording has 1 segment of 40s each # run mountainsort4 in mono-segment mode @@ -465,6 +468,7 @@ Here is the list of external sorters accessible using the run_sorter wrapper: * **Kilosort2** :code:`run_sorter(sorter_name='kilosort2')` * **Kilosort2.5** :code:`run_sorter(sorter_name='kilosort2_5')` * **Kilosort3** :code:`run_sorter(sorter_name='kilosort3')` +* **Kilosort4** :code:`run_sorter(sorter_name='kilosort4')` * **PyKilosort** :code:`run_sorter(sorter_name='pykilosort')` * **Klusta** :code:`run_sorter(sorter_name='klusta')` * **Mountainsort4** :code:`run_sorter(sorter_name='mountainsort4')` @@ -483,24 +487,18 @@ experimental for now: * **Spyking Circus2** :code:`run_sorter(sorter_name='spykingcircus2')` * **Tridesclous2** :code:`run_sorter(sorter_name='tridesclous2')` -In 2024, we expect to add many more sorters to this list. - - Installed Sorters ----------------- To check which sorters are useable in a given Python environment, one can print the installed -sorters list. An example is shown in a pre-defined miniconda3 environment. - - -Then you can check the installed Sorter list, +sorters list: .. code:: python from spikeinterface.sorters import installed_sorters installed_sorters() -which outputs, +which, in our case, outputs .. parsed-literal:: ['herdingspikes', @@ -551,6 +549,7 @@ From the user's perspective, they behave exactly like the external sorters: sorting = run_sorter(sorter_name="spykingcircus2", recording=recording, folder="/tmp/folder") +Read more in the :ref:`sorting-components-module` docs. Contributing ------------ diff --git a/doc/modules/sortingcomponents.rst b/doc/modules/sortingcomponents.rst index a32e111bd7..2a80ee29e7 100644 --- a/doc/modules/sortingcomponents.rst +++ b/doc/modules/sortingcomponents.rst @@ -1,3 +1,5 @@ +.. _sorting-components-module: + Sorting Components module ========================= @@ -52,11 +54,9 @@ follows: peak_sign='neg', detect_threshold=5, exclude_sweep_ms=0.2, - radius_um=100, noise_levels=None, random_chunk_kwargs={}, - outputs='numpy_compact', - engine='numpy', + gather_mode='memory', **job_kwargs, ) @@ -95,9 +95,15 @@ follows: job_kwargs = dict(chunk_duration='1s', n_jobs=8, progress_bar=True) - peak_locations = localize_peaks(recording=recording, peaks=peaks, method='center_of_mass', - radius_um=70., ms_before=0.3, ms_after=0.6, - **job_kwargs) + peak_locations = localize_peaks( + recording=recording, + peaks=peaks, + method='center_of_mass', + radius_um=70., + ms_before=0.3, + ms_after=0.6, + **job_kwargs + ) Currently, the following methods are implemented: @@ -165,19 +171,11 @@ Motion estimation ----------------- Recently, drift estimation has been added to some of the available spike sorters (Kilosort 2.5, 3) -Especially for Neuropixels-like probes, this is a crucial step. - -Several methods have been proposed to correct for drift, but only one is currently implemented in SpikeInterface. -See `Decentralized Motion Inference and Registration of Neuropixel Data `_ -for more details. - -The motion estimation step comes after peak detection and peak localization. -The idea is to divide the recording into time bins and estimate the relative motion between temporal bins. +Especially for acute Neuropixels-like probes, this is a crucial step. -This method has two options: - - * rigid drift : one motion vector is estimated for the entire probe - * non-rigid drift : one motion vector is estimated per depth bin +The motion estimation step comes after peak detection and peak localization. Read more about +it in the :ref:`_motion_correction` modules doc, and a more practical guide in the +:ref:`handle-drift-with-spikeinterface` How To. Here is an example with non-rigid motion estimation: @@ -191,20 +189,29 @@ Here is an example with non-rigid motion estimation: from spikeinterface.sortingcomponents.motion import estimate_motion - motion, temporal_bins, spatial_bins, - extra_check = estimate_motion(recording=recording, peaks=peaks, peak_locations=peak_locations, - direction='y', bin_s=10., bin_um=10., margin_um=0., - method='decentralized_registration', - rigid=False, win_shape='gaussian', win_step_um=50., win_sigma_um=150., - progress_bar=True, verbose=True) + motion = estimate_motion( + recording=recording, + peaks=peaks, + peak_locations=peak_locations, + direction='y', + bin_s=10., + bin_um=10., + margin_um=0., + method='decentralized', + rigid=False, + win_shape='gaussian', + win_step_um=50., + progress_bar=True, + verbose=True + ) In this example, because it is a non-rigid estimation, :code:`motion` is a 2d array (num_time_bins, num_spatial_bins). - +We could now check the ``motion`` object and see if we need to apply a correction. Motion interpolation -------------------- -The estimated motion can be used to interpolate traces, in other words, for drift correction. +The estimated motion can be used to interpolate traces to attempt to correct for drift. One possible way is to make an interpolation sample-by-sample to compensate for the motion. The :py:class:`~spikeinterface.sortingcomponents.motion.InterpolateMotionRecording` is a preprocessing step doing this. This preprocessing is *lazy*, so that interpolation is done on-the-fly. However, the class needs the @@ -213,14 +220,16 @@ estimation). Here is a short example that depends on the output of "Motion interpolation": - .. code-block:: python from spikeinterface.sortingcomponents.motion import InterpolateMotionRecording - recording_corrected = InterpolateMotionRecording(recording=recording_with_drift, motion=motion, temporal_bins=temporal_bins, spatial_bins=spatial_bins - spatial_interpolation_method='kriging', - border_mode='remove_channels') + recording_corrected = InterpolateMotionRecording( + recording=recording_with_drift, + motion=motion, + spatial_interpolation_method='kriging', + border_mode='remove_channels' + ) **Notes**: * :code:`spatial_interpolation_method` "kriging" or "iwd" do not play a big role. diff --git a/doc/modules/widgets.rst b/doc/modules/widgets.rst index cc6aa63d80..6a2e11d417 100644 --- a/doc/modules/widgets.rst +++ b/doc/modules/widgets.rst @@ -6,17 +6,14 @@ Widgets module The :py:mod:`spikeinterface.widgets` module includes plotting function to visualize recordings, sortings, waveforms, and more. -Since version 0.95.0, the :py:mod:`spikeinterface.widgets` module supports multiple backends: +The :py:mod:`spikeinterface.widgets` module supports multiple backends: * | :code:`matplotlib`: static rendering using the `matplotlib `_ package * | :code:`ipywidgets`: interactive rendering within a jupyter notebook using the | `ipywidgets `_ package * | :code:`sortingview`: web-based and interactive rendering using the `sortingview `_ | and `FIGURL `_ packages. - -Version 0.99.0 also comes with this new backend: - -* :code:`ephyviewer`: interactive Qt based using the `ephyviewer `_ package +* | :code:`ephyviewer`: interactive Qt based using the `ephyviewer `_ package Installing backends @@ -44,8 +41,8 @@ To install it, run: Install ipywidgets ^^^^^^^^^^^^^^^^^^ -The :code:`ipywidgets` backend allows users to interact with the plot, for example, by selecting units or -scrolling through a time series. +The :code:`ipywidgets` backend allows users to interact with the plot when you plot it in +a Jupyter Notebook. For example, in certain widgets you can select units or scroll through a time series. To install it, run: @@ -94,7 +91,6 @@ Install ephyviewer This backend is Qt based with PyQt5, PyQt6 or PySide6 support. Qt is sometimes tedious to install. - For a pip-based installation, run: .. code-block:: bash @@ -122,15 +118,17 @@ A default backend for a SpikeInterface session can be set with the .. code-block:: python + import spikeinterface.widgets as sw + # matplotlib backend - set_default_plotter_backend(backend="ipywidgets") - print(get_default_plotter_backend()) + sw.set_default_plotter_backend(backend="ipywidgets") + print(sw.get_default_plotter_backend()) # >>> "ipywidgets" All :code:`plot_*` functions return a :code:`BackendPlotter` instance. Different backend-specific plotters can expose different attributes. For example, the :code:`matplotlib` plotter has the :code:`figure`, :code:`ax`, and :code:`axes` (for multi-axes plots) attributes to enable further -customization. +customization. To learn more about customization read our how to guide, :ref:`customize-a-plot`. matplotlib @@ -186,7 +184,7 @@ sortingview The :code:`plot_*(..., backend="sortingview")` generate web-based GUIs, which are also shareable with a link (provided that :code:`kachery-cloud` is correctly setup, see :ref:`sorting_view`). -The functions have the following additional arguments: +The functions has the following additional arguments: * :code:`generate_url`: If True, the figurl URL is generated and printed. Default True * :code:`display`: If True and in jupyter notebook/lab, the widget is displayed in the cell. Default True @@ -269,14 +267,16 @@ Available plotting functions * :py:func:`~spikeinterface.widgets.plot_crosscorrelograms` (backends: :code:`matplotlib`, :code:`sortingview`) * :py:func:`~spikeinterface.widgets.plot_isi_distribution` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_motion` (backends: :code:`matplotlib`) +* :py:func:`~spikeinterface.widgets.plot_drift_raster_map` (backends: :code:`matplotlib`, :code:`ipywidgets`) * :py:func:`~spikeinterface.widgets.plot_multicomparison_agreement` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_multicomparison_agreement_by_sorter` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_multicomparison_graph` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_peak_activity` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_probe_map` (backends: :code:`matplotlib`) * :py:func:`~spikeinterface.widgets.plot_quality_metrics` (backends: :code:`matplotlib`, :code:`ipywidgets`, :code:`sortingview`) -* :py:func:`~spikeinterface.widgets.plot_rasters` (backends: :code:`matplotlib`) +* :py:func:`~spikeinterface.widgets.plot_rasters` (backends: :code:`matplotlib`, :code:`ipywidgets`) * :py:func:`~spikeinterface.widgets.plot_sorting_summary` (backends: :code:`sortingview`) +* :py:func:`~spikeinterface.widgets.plot_spike_amplitudes` (backends: :code:`matplotlib`, :code:`ipywidgets`) * :py:func:`~spikeinterface.widgets.plot_spike_locations` (backends: :code:`matplotlib`, :code:`ipywidgets`) * :py:func:`~spikeinterface.widgets.plot_spikes_on_traces` (backends: :code:`matplotlib`, :code:`ipywidgets`) * :py:func:`~spikeinterface.widgets.plot_template_metrics` (backends: :code:`matplotlib`, :code:`ipywidgets`, :code:`sortingview`) From 0e0903304635a90ccfa0dcb107e9ca5fb3a01c41 Mon Sep 17 00:00:00 2001 From: chrishalcrow Date: Mon, 11 Aug 2025 10:03:37 +0100 Subject: [PATCH 2/7] Update core.rst --- doc/how_to/physical_units.rst | 2 ++ doc/modules/core.rst | 34 +++++++++++++++++++++------------- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/doc/how_to/physical_units.rst b/doc/how_to/physical_units.rst index d2c1743930..28cae232dd 100644 --- a/doc/how_to/physical_units.rst +++ b/doc/how_to/physical_units.rst @@ -1,3 +1,5 @@ +.. _working-with-physical-units: + Working with physical units in SpikeInterface recordings ======================================================== diff --git a/doc/modules/core.rst b/doc/modules/core.rst index ef26fb14c7..54f2c59466 100644 --- a/doc/modules/core.rst +++ b/doc/modules/core.rst @@ -76,9 +76,9 @@ with 16 channels: recording_slice_frames = recording.frame_slice(start_frame=0, end_frame=int(10*sampling_frequency)) # get new recording with the first 4 channels - recording_slice_chans = recording.channel_slice(channel_ids=channel_ids[:4]) + recording_slice_chans = recording.select_channels(channel_ids=channel_ids[:4]) # remove last two channels - recording_rm_chans = recording.remove_channels(channel_ids=channel_ids[-2:]) + recording_rm_chans = recording.remove_channels(remove_channel_ids=channel_ids[-2:]) # set channel grouping (assume we have 4 groups of 4 channels, e.g. tetrodes) groups = [0] * 4 + [1] * 4 + [2] * 4 + [3] * 4 @@ -89,13 +89,14 @@ with 16 channels: # sliced recordings as values # set times (for synchronization) - assume our times start at 300 seconds + num_samples = recording.get_num_samples() timestamps = np.arange(num_samples) / sampling_frequency + 300 recording.set_times(times=timestamps, segment_index=0) **Note**: Raw data formats often store data as integer values for memory efficiency. To give these integers meaningful physical units (uV), you can apply a gain and an offset. Many devices have their own gains and offsets necessary to convert their data and these values are handled by SpikeInterface for its extractors. This -is triggered by the :code:`return_scaled` parameter in :code:`get_traces()`, (see above example), which will return the traces in uV. +is triggered by the :code:`return_in_uV` parameter in :code:`get_traces()`, (see above example), which will return the traces in uV. Read more in our how to guide, :ref:`working-with-physical-units`. Sorting @@ -118,7 +119,7 @@ with 10 units: .. code-block:: python unit_ids = sorting.unit_ids - num_channels = sorting.get_num_units() + num_units = sorting.get_num_units() sampling_frequency = sorting.sampling_frequency # retrieve spike trains for a unit (returned as sample indices) @@ -178,8 +179,8 @@ to perform further analysis, such as calculating :code:`waveforms` and :code:`te Importantly, the :py:class:`~spikeinterface.core.SortingAnalyzer` handles the *sparsity* and the physical *scaling*. Sparsity defines the channels on which waveforms and templates are calculated using, for example, a physical distance from the channel with the largest peak amplitude (see the :ref:`modules/core:Sparsity` section). Scaling, set by -the :code:`return_scaled` argument, determines whether the data is converted from integer values to :math:`\mu V` or not. -By default, :code:`return_scaled` is true and all processed data voltage values are in :math:`\mu V` (e.g., waveforms, templates, spike amplitudes, etc.). +the :code:`return_in_uV` argument, determines whether the data is converted from integer values to :math:`\mu V` or not. +By default, :code:`return_in_uV` is true and all processed data voltage values are in :math:`\mu V` (e.g., waveforms, templates, spike amplitudes, etc.). Now we will create a :code:`SortingAnalyzer` called :code:`sorting_analyzer`. @@ -217,6 +218,7 @@ is run again with one of the backends supplied. .. code-block:: python + from pathlib import Path # create a "processed" folder processed_folder = Path("processed") @@ -266,8 +268,7 @@ The :code:`sorting_analyzer` object implements convenient functions to access th num_channels = sorting_analyzer.get_num_channels() num_units = sorting_analyzer.get_num_units() - sampling_frequency = sorting_analyzer.get_sampling_frequency() - # or: sampling_frequency = sorting_analyzer.sampling_frequency + sampling_frequency = sorting_analyzer.sampling_frequency total_num_samples = sorting_analyzer.get_total_samples() total_duration = sorting_analyzer.get_total_duration() @@ -689,6 +690,7 @@ In this example, we create a recording and a sorting object from numpy objects: .. code-block:: python import numpy as np + from spikeinterface.core import NumpyRecording, NumpySorting # in-memory recording sampling_frequency = 30_000. @@ -697,7 +699,10 @@ In this example, we create a recording and a sorting object from numpy objects: num_channels = 16 random_traces = np.random.randn(num_samples, num_channels) - recording_memory = NumpyRecording(traces_list=[random_traces]) + recording_memory = NumpyRecording( + traces_list=[random_traces] + sampling_frequency=sampling_frequency, + ) # with more elements in `traces_list` we can make multi-segment objects # in-memory sorting @@ -706,7 +711,7 @@ In this example, we create a recording and a sorting object from numpy objects: spike_trains = [] labels = [] for i in range(num_units): - spike_trains_i = np.random.randint(low=0, high=num_samples, size=num_spikes_unit) + spike_trains_i = list(np.random.randint(low=0, high=num_samples, size=num_spikes_unit)) labels_i = [i] * num_spikes_unit spike_trains += spike_trains_i labels += labels_i @@ -751,7 +756,7 @@ the new objects will be *views* of the original ones. # keep one channel of every tenth channel keep_ids = recording.channel_ids[::10] - sub_recording = recording.channel_slice(channel_ids=keep_ids) + sub_recording = recording.select_channels(channel_ids=keep_ids) # keep between 5min and 12min fs = recording.sampling_frequency @@ -870,13 +875,15 @@ They are useful to make examples, tests, and small demos: .. code-block:: python + from spikeinterface.core import generate_recording, generate_sorting, generate_snippets + # recording with 2 segments and 4 channels recording = generate_recording(num_channels=4, sampling_frequency=30000., durations=[10.325, 3.5], set_probe=True) # sorting with 2 segments and 5 units sorting = generate_sorting(num_units=5, sampling_frequency=30000., durations=[10.325, 3.5], - firing_rate=15, refractory_period=1.5) + firing_rates=15, refractory_period_ms=1.5) # snippets of 60 samples on 2 channels from 5 units snippets = generate_snippets(nbefore=20, nafter=40, num_channels=2, @@ -929,7 +936,8 @@ WaveformExtractor ^^^^^^^^^^^^^^^^^ This is now a legacy object that can still be accessed through the :py:class:`MockWaveformExtractor`. It is kept -for backward compatibility. +for backward compatibility. You can convert a ``WaveformExtractor`` to a ``SortingAnalyzer`` +easily, :ref:`using this guide `. The :py:class:`~spikeinterface.core.WaveformExtractor` class is the core object to combine a :py:class:`~spikeinterface.core.BaseRecording` and a :py:class:`~spikeinterface.core.BaseSorting` object. From 35712be864afd1501d481b7a2beb563bd3e2cf44 Mon Sep 17 00:00:00 2001 From: chrishalcrow Date: Mon, 11 Aug 2025 10:12:21 +0100 Subject: [PATCH 3/7] update to handle-drift-in-your-recording --- doc/modules/motion_correction.rst | 2 +- doc/modules/sortingcomponents.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index be41b85f4b..4066fb5357 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -4,7 +4,7 @@ Motion/drift correction ======================= -See a practical guide to motion correction in our How To guide: :ref:`handle-drift-with-spikeinterface`. +See a practical guide to motion correction in our How To guide: :ref:`handle-drift-in-your-recording`. Overview -------- diff --git a/doc/modules/sortingcomponents.rst b/doc/modules/sortingcomponents.rst index 2a80ee29e7..dc80410fa8 100644 --- a/doc/modules/sortingcomponents.rst +++ b/doc/modules/sortingcomponents.rst @@ -175,7 +175,7 @@ Especially for acute Neuropixels-like probes, this is a crucial step. The motion estimation step comes after peak detection and peak localization. Read more about it in the :ref:`_motion_correction` modules doc, and a more practical guide in the -:ref:`handle-drift-with-spikeinterface` How To. +:ref:`handle-drift-in-your-recording` How To. Here is an example with non-rigid motion estimation: From a9af25da9515720bba304861931d1061fb95d9a0 Mon Sep 17 00:00:00 2001 From: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com> Date: Mon, 18 Aug 2025 09:38:47 +0100 Subject: [PATCH 4/7] Update doc/modules/motion_correction.rst Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com> --- doc/modules/motion_correction.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index 4066fb5357..2962ce6343 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -62,7 +62,7 @@ One challenging task for motion correction is to determine the parameters. The high level :py:func:`~spikeinterface.preprocessing.correct_motion()` proposes the concept of a **"preset"** that already has predefined parameters, in order to achieve a calibrated behavior. -We currently have presets: +We currently have these presets: * **"dredge"**: The official implementation of DREDGE, used in [Windolf2023]_. * **"dredge_fast"**: A faster implementation of DREDGE, which should give similar results. From 5a1b7f8b76dc4c3e81ef974a58a103ecc3f70e96 Mon Sep 17 00:00:00 2001 From: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com> Date: Mon, 18 Aug 2025 09:38:53 +0100 Subject: [PATCH 5/7] Update doc/modules/motion_correction.rst Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com> --- doc/modules/motion_correction.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index 2962ce6343..867ff93897 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -126,7 +126,7 @@ Optionally any parameter from the preset can be overwritten: ) Importantly, all the result and intermediate computations can returned to a motion object, for further loading, -verification and visualisation. +verification and visualization. .. code-block:: python From 9bd0439369c05c14a83d7bc6cd33dd7972f1f2e4 Mon Sep 17 00:00:00 2001 From: chrishalcrow Date: Mon, 18 Aug 2025 13:56:06 +0100 Subject: [PATCH 6/7] respond to zach --- doc/modules/curation.rst | 4 ++-- doc/modules/motion_correction.rst | 4 ++-- doc/modules/postprocessing.rst | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/modules/curation.rst b/doc/modules/curation.rst index d1e45c80af..e357c82124 100644 --- a/doc/modules/curation.rst +++ b/doc/modules/curation.rst @@ -34,7 +34,7 @@ a subset of units from a spike sorting output and to perform some merges: remove_unit_ids = [1, 2] sorting_analyzer2 = sorting_analyzer.remove_units(remove_unit_ids=remove_unit_ids) - # merge units 4 and 5, and units 7, 8 and 12 + # merge units 4 and 5, and separately merge units 7, 8 and 12 merge_unit_groups = [[4, 5], [7, 8, 12]] sorting_analyzer3 = sorting_analyzer2.merge_units( merge_unit_groups=merge_unit_groups, @@ -190,7 +190,7 @@ prone to wrong merges. To do so, you'll need to do the following: The extra keyword ``recursive`` specifies that for each presets/sequences of steps, merges are performed until no further merges are possible. The ``job_kwargs`` are the parameters for the parallelization. -**Be careful: the merges can not be reverted, so be sure to not erase your analyzer and instead create a one** +**Be careful: the merges can not be reverted, so be sure to not erase your analyzer and instead create a new one** Manual curation diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index 4066fb5357..7c7b93a43b 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -125,7 +125,7 @@ Optionally any parameter from the preset can be overwritten: ) ) -Importantly, all the result and intermediate computations can returned to a motion object, for further loading, +Importantly, all the results and intermediate computations can returned to a motion object, for further loading, verification and visualisation. .. code-block:: python @@ -137,7 +137,7 @@ verification and visualisation. plot_motion(motion) Alternatively, you can save the motion (and related motion info) in a folder. The folder will contain -the motion vector itself, as well as detected peaks, peak location, and more. +the motion vector itself, as well as detected peaks, peak locations, and more. .. code-block:: python diff --git a/doc/modules/postprocessing.rst b/doc/modules/postprocessing.rst index cbf799c5e3..66cd16a9c5 100644 --- a/doc/modules/postprocessing.rst +++ b/doc/modules/postprocessing.rst @@ -254,7 +254,7 @@ spike_amplitudes ^^^^^^^^^^^^^^^^ This extension computes the amplitude of each spike as the value of the traces on the extremum channel at the times of -each spike. The extremum channel is compute from the templates. +each spike. The extremum channel is computed from the templates. **NOTE:** computing spike amplitudes is highly recommended before calculating amplitude-based quality metrics, such as :ref:`amp_cutoff` and :ref:`amp_median`. From 36d04bfb78e8b9ddd2467f8c7aeee79023aa10a2 Mon Sep 17 00:00:00 2001 From: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com> Date: Fri, 22 Aug 2025 09:24:06 +0100 Subject: [PATCH 7/7] Update doc/modules/motion_correction.rst Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com> --- doc/modules/motion_correction.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/modules/motion_correction.rst b/doc/modules/motion_correction.rst index 6bf8104593..a7f687c6ec 100644 --- a/doc/modules/motion_correction.rst +++ b/doc/modules/motion_correction.rst @@ -125,7 +125,7 @@ Optionally any parameter from the preset can be overwritten: ) ) -Importantly, all the results and intermediate computations can returned to a motion object, for further loading, +Importantly, all the results and intermediate computations can be returned to a motion object, for further loading, verification and visualization. .. code-block:: python