k-dense-ai/scientific-agent-skills/neuropixels-analysis

Neuropixels Data Analysis

Overview

Toolkit for analyzing Neuropixels high-density neural recordings using current best practices from SpikeInterface, the Allen Institute, and the International Brain Laboratory (IBL). It covers the full workflow from raw data to reviewed, curated units. Targets SpikeInterface 0.105.0, ProbeInterface 0.4.0 and Neo 0.14.5 (reviewed 2026-10-01). Synthetic tests cover recording contracts, preprocessing, analyzers, metrics and exports; real acquisition files, native sorters, GPU execution and pretrained models remain illustrative.

All examples use the real SpikeInterface API (spikeinterface.full as si) plus the companion curation module (spikeinterface.curation as sc). The skill ships runnable scripts in scripts/ and a copy-and-edit template in assets/ that implement this workflow directly on top of SpikeInterface — there is no separate package to install beyond the dependencies listed under Installation.

When to Use This Skill

This skill should be used when:

• Working with Neuropixels recordings (.ap.bin, .lf.bin, .meta files)
• Loading data from SpikeGLX, Open Ephys, or NWB formats
• Preprocessing neural recordings (filtering, common reference, bad-channel detection)
• Detecting and correcting motion/drift
• Running spike sorting (Kilosort4, SpykingCircus2, Mountainsort5, Tridesclous2)
• Computing quality metrics (SNR, ISI violations, presence ratio, amplitude cutoff)
• Curating units (threshold-based, model-based, or AI-assisted)
• Creating visualizations and exporting to Phy or NWB

Supported Hardware & Formats

Quick Start

Import and configure parallel processing

import spikeinterface.full as si
# Global job kwargs are reused by all parallelizable steps
si.set_global_job_kwargs(n_jobs=-1, chunk_duration="1s", progress_bar=True)

Loading data

# Inspect available streams first
stream_names, stream_ids = si.get_neo_streams("spikeglx", "/path/to/run_g0/")
print(stream_names)  # e.g. ['imec0.ap', 'imec0.lf', 'nidq']
# SpikeGLX (most common) — select the AP stream by name
recording = si.read_spikeglx("/path/to/run_g0/", stream_name="imec0.ap")
# Open Ephys
recording = si.read_openephys("/path/to/Record_Node_101/")
# For quick iteration, slice the first 60 s
fs = recording.get_sampling_frequency()
recording_sub = recording.frame_slice(0, min(int(60 * fs), recording.get_num_samples()))

Full pipeline (bundled script)

The repository ships an end-to-end pipeline built on SpikeInterface:

python scripts/neuropixels_pipeline.py /path/to/spikeglx/data output/ --sorter kilosort4 --curation allen

It performs load → preprocess → drift check → optional motion correction → sorting → postprocessing → quality metrics → curation → export. Phy and the report retain all units for review with curation labels; sorting_curated/ contains selected good units. Read the steps below to run them interactively or customize the pipeline.

Standard Analysis Workflow

1. Preprocessing

Validate AP stream, calibration, channel order, probe geometry and segment boundaries first. This filter/reference chain is not full IBL destriping. Apply ADC timing correction only from valid acquisition metadata and reference each shank separately:

rec = si.highpass_filter(recording, freq_min=400.0)
bad_channel_ids, channel_labels = si.detect_bad_channels(rec)
rec = rec.remove_channels(bad_channel_ids)
rec = si.phase_shift(rec)  # Requires valid inter_sample_shift property.
rec = si.common_reference(rec, operator="median", reference="global")  # Single shank only.

For multiple shanks, use the bundled reference_by_shank helper described in PREPROCESSING.md. Trace arrays are samples × channels; get_traces(return_in_uV=True) needs calibrated gain/offset. The bundled commands reject uncalibrated, empty or multisegment input instead of guessing.

Cache preprocessed data when storage and repeated use justify it:

rec = rec.save(folder="preprocessed/", format="binary")

2. Check and correct drift

Always inspect drift before sorting:

from spikeinterface.sortingcomponents.peak_detection import detect_peaks
from spikeinterface.sortingcomponents.peak_localization import localize_peaks
noise_levels = si.get_noise_levels(rec, return_in_uV=False)
peaks = detect_peaks(rec, method='locally_exclusive', method_kwargs={'noise_levels': noise_levels, 'detect_threshold': 5, 'radius_um': 50.0})
peak_locations = localize_peaks(rec, peaks, method="center_of_mass")
# Visualize the drift raster
si.plot_drift_raster_map(peaks=peaks, peak_locations=peak_locations,
                         recording=rec, clim=(-50, 50))

Apply correction if needed (presets: rigid_fast, kilosort_like, nonrigid_accurate, nonrigid_fast_and_accurate, dredge, dredge_fast, medicine):

rec_corrected = si.correct_motion(rec, preset="nonrigid_fast_and_accurate", folder="motion/")

3. Spike sorting

The calls below are illustrative until tested on the target recording and sorter. Choose one drift-correction stage: externally corrected input uses do_correction=False for Kilosort 2.5/3/4, or apply_motion_correction=False for Spykingcircus2. These flags match SpikeInterface 0.105.0; inspect sorter parameters when using another release. Uncorrected input can use sorter defaults. A spread of peak depths across neurons is not a temporal drift estimate. The bundled pipeline estimates/corrects motion when requested; inspect its saved motion output.

# Kilosort4 (external install; CUDA recommended, CPU mode also supported)
sorting = si.run_sorter("kilosort4", rec_corrected, folder="ks4_output", do_correction=False)
# CPU alternatives (SC2/TDC2 need SI optional dependencies; MS5 needs mountainsort5)
sorting = si.run_sorter("spykingcircus2", rec_corrected, folder="sc2_output", apply_motion_correction=False)
sorting = si.run_sorter("tridesclous2", rec_corrected, folder="tdc2_output")
sorting = si.run_sorter("mountainsort5", rec_corrected, folder="ms5_output")
# External sorters can run in containers without local install
sorting = si.run_sorter("kilosort2_5", rec_corrected, folder="ks25_output", docker_image=True, do_correction=False)
print(si.installed_sorters())

Note: run_sorter uses the folder= argument. The older output_folder= is deprecated.

4. Postprocessing

analyzer = si.create_sorting_analyzer(sorting, rec_corrected, sparse=True,
                                      format="binary_folder", folder="analyzer/")
analyzer.compute("random_spikes", method="uniform", max_spikes_per_unit=500)
analyzer.compute("waveforms", ms_before=1.0, ms_after=2.0)
analyzer.compute("templates", operators=["average", "std"])
analyzer.compute("noise_levels")
analyzer.compute("spike_amplitudes")
analyzer.compute("amplitude_scalings")
analyzer.compute("correlograms", window_ms=50.0, bin_ms=1.0)
analyzer.compute("unit_locations", method="monopolar_triangulation")
analyzer.compute("template_similarity")
metric_names = ["firing_rate", "presence_ratio", "snr", "isi_violation", "amplitude_cutoff"]
analyzer.compute("quality_metrics", metric_names=metric_names)
metrics = analyzer.get_extension("quality_metrics").get_data()

5. Curation by metric thresholds

# Example screen, not a guarantee of single-neuron isolation.
query = "(amplitude_cutoff < 0.1) & (isi_violations_ratio < 0.5) & (presence_ratio > 0.9)"
good_unit_ids = metrics.query(query).index.values

For reusable local screening with allen / legacy ibl / strict presets, use the bundled scripts/compute_metrics.py. See references/AUTOMATED_CURATION.md for details and the Bombcell / UnitMatch tools. The legacy ibl preset is not the IBL classifier. Missing/nonfinite metrics remain unsorted, and boundary values fail the strict thresholds. All bundled entry points now use the same screening criteria.

6. Model-based curation (UnitRefine)

SpikeInterface can apply pretrained machine-learning classifiers from Hugging Face via the spikeinterface.curation module. The UnitRefine models were trained on real Neuropixels data (V1, SC, ALM). The public model metadata currently requests SI 0.102.0 and scikit-learn 1.4.2, with empty metric-parameter metadata; compatibility with this 0.105.0 environment is untested. Inspect model requirements/features first:

import spikeinterface.curation as sc
# 1) noise vs neural
noise_labels = sc.model_based_label_units(
    sorting_analyzer=analyzer,
    repo_id="SpikeInterface/UnitRefine_noise_neural_classifier",
    trust_model=True,
    enforce_metric_params=True,
)
neural = analyzer.remove_units(noise_labels[noise_labels["prediction"] == "noise"].index)
# 2) single-unit (sua) vs multi-unit (mua) on the surviving units
sua_mua_labels = sc.model_based_label_units(
    sorting_analyzer=neural,
    repo_id="SpikeInterface/UnitRefine_sua_mua_classifier",
    trust_model=True,
    enforce_metric_params=True,
)

Each call returns a DataFrame with prediction and probability (confidence) per unit. trust_model=True (or an explicit trusted=[...] list) is required to load the .skops model — only load models from sources you trust. Parameter enforcement cannot validate training settings absent from model metadata. Models trained on other brain areas/datasets may not transfer; validate against a manually labelled subset.

7. AI-assisted curation (for uncertain units)

When running inside an agent such as Cursor or Claude Code, the agent can directly inspect waveform/correlogram plots and suggest review questions — no API setup required. Generate plots and ask the agent to assess isolation quality.

For programmatic vision-model access, read API keys from the environment — never hardcode credentials in analysis scripts (they leak into version control and logs):

import os
from anthropic import Anthropic
client = Anthropic(api_key=os.environ["ANTHROPIC_API_KEY"])  # set this in your shell, not in code

See references/AI_CURATION.md for the full pattern (rendering a unit summary image, building the prompt, and retaining the response as advisory evidence).

8. Export results

# Keep only good units, then export
analyzer_clean = analyzer.select_units(good_unit_ids, folder="analyzer_clean/", format="binary_folder")
# Phy for manual review
si.export_to_phy(analyzer_clean, output_folder="phy_export/",
                 compute_pc_features=True, compute_amplitudes=True)
# Figures report
si.export_report(analyzer_clean, "report/", format="png")
# Metrics table
metrics.to_csv("quality_metrics.csv")

SpikeInterface 0.105.0 has no export_to_nwb exporter. Use the source-specific NeuroConv NWBConverter workflow, with session metadata, electrodes, calibration and aligned unit times; validate the resulting NWB file. This optional conversion was not executed here.