/neuropixels-analysis

name: "neuropixels-analysis"
description: "Pipeline for Neuropixels extracellular electrophysiology: probe geometry (ProbeInterface), Kilosort sorting via SpikeInterface, quality metrics, unit curation (ISI, firing rate, SNR), post-sort analysis (PSTH, tuning curves, population decoding). Supports Neuropixels 1.0/2.0/Ultra in rodent/primate experiments."
license: "MIT"

Neuropixels Analysis

Overview

Neuropixels probes record extracellular voltage from 384 (NP 1.0) or 192 (NP 2.0) simultaneously recorded channels at 30 kHz. Analysis follows a canonical pipeline: raw data → spike sorting (Kilosort) → quality curation → unit analysis. **SpikeInterface** (Python) provides a unified API across 10+ spike sorters, handles data loading from multiple formats (SpikeGLX, OpenEphys, NWB), computes quality metrics, and exports sorted results. **ProbeInterface** manages probe geometry and channel maps. Post-sort analysis (PSTHs, firing rate, decoding) uses standard Python scientific stack.

When to Use

• Spike-sorting Neuropixels recordings from SpikeGLX (`.bin`) or OpenEphys (`.dat`) to extract single-unit activity
• Applying automatic quality control metrics (ISI violations, SNR, firing rate) to curate sorted units
• Computing peristimulus time histograms (PSTHs) locked to experimental events
• Analyzing population coding: decoding stimulus or behavioral variables from firing rates
• Converting sorted data to NWB (Neurodata Without Borders) format for sharing
• Comparing multiple spike sorters on the same dataset for method validation
• Visualizing unit waveforms, auto-correlograms, and spatial distribution across probe channels
• Use **SpikeInterface** instead for a unified framework that supports 10+ spike sorters with a common API and comparison tools

Prerequisites

• **Python packages**: `spikeinterface[full]`, `probeinterface`, `numpy`, `pandas`, `matplotlib`, `scipy`
• **Spike sorter**: Kilosort4 (Python, `pip install kilosort`); or MATLAB-based Kilosort2/3 (requires MATLAB license)
• **Data requirements**: raw `.bin` (SpikeGLX) or `.dat` (OpenEphys) recording; probe channel map file (`.prb` or from ProbeInterface)
• **Hardware**: NVIDIA GPU (4+ GB VRAM) required for Kilosort; CPU fallback available but ~10× slower

pip install spikeinterface[full] probeinterface kilosort
# Optional: Phy for manual curation
pip install phy

Workflow

Step 1: Load Raw Recording

import spikeinterface.full as si
from pathlib import Path

# SpikeGLX recording (most common Neuropixels format)
data_dir = Path("/data/recording/session_001")
recording = si.read_spikeglx(data_dir, stream_name="imec0.ap")

print(f"Probe type:       {recording.get_probe().name}")
print(f"Channels:         {recording.get_num_channels()}")
print(f"Sampling rate:    {recording.get_sampling_frequency()} Hz")
print(f"Duration:         {recording.get_total_duration():.1f} s")
print(f"Total samples:    {recording.get_total_samples()}")

Step 2: Apply Common Reference and Bandpass Filter

import spikeinterface.preprocessing as spre

# Common median reference (removes common noise across all channels)
recording_cmr = spre.common_reference(recording, reference="global", operator="median")

# Bandpass filter: 300–6000 Hz for spikes
recording_filt = spre.bandpass_filter(recording_cmr, freq_min=300, freq_max=6000)

# Remove bad channels automatically
recording_clean, removed_ids = spre.remove_bad_channels(recording_filt)
print(f"Removed {len(removed_ids)} bad channels: {removed_ids}")
print(f"Clean recording: {recording_clean.get_num_channels()} channels")

Step 3: Run Spike Sorting

import spikeinterface.sorters as ss
from pathlib import Path

output_dir = Path("./sorting_output")

# Kilosort4 (recommended for Neuropixels; requires GPU)
sorting = ss.run_sorter(
    "kilosort4",
    recording_clean,
    output_folder=output_dir / "kilosort4",
    remove_existing_folder=True,
    verbose=True,
    # Kilosort4 parameters
    nblocks=5,          # number of drift correction blocks
    Th_learned=8,       # threshold for learned templates
    do_correction=True, # drift correction
)

print(f"Units found: {len(sorting.get_unit_ids())}")
print(f"Spike counts (first 5): {[len(sorting.get_unit_spike_train(u, segment_index=0)) for u in sorting.unit_ids[:5]]}")

Step 4: Compute Waveforms and Quality Metrics

import spikeinterface.full as si
import spikeinterface.qualitymetrics as sqm

# Extract waveforms (snippets around each spike)
waveforms = si.extract_waveforms(
    recording_clean,
    sorting,
    folder="./waveforms",
    ms_before=1.0,      # ms before spike peak
    ms_after=2.0,       # ms after spike peak
    max_spikes_per_unit=1000,
    overwrite=True,
    n_jobs=4,
)

print(f"Waveform shape per unit: {waveforms.get_waveforms(waveforms.unit_ids[0]).shape}")
# (n_spikes, n_samples, n_channels)

# Compute quality metrics
metrics = sqm.compute_quality_metrics(
    waveforms,
    metric_names=["snr", "isi_violation", "firing_rate", "presence_ratio",
                  "amplitude_cutoff", "nearest_neighbor"],
)
print(f"\nQuality metrics summary:")
print(metrics.describe())

Step 5: Curate Units

import pandas as pd

# Curation thresholds (Allen Brain Institute defaults)
thresholds = {
    "snr":                (5.0, None),   # SNR ≥ 5
    "isi_violations_ratio": (None, 0.1), # ISI violation ratio ≤ 10%
    "firing_rate":        (0.1, None),   # firing rate ≥ 0.1 Hz
    "presence_ratio":     (0.9, None),   # present ≥ 90% of recording
    "amplitude_cutoff":   (None, 0.1),   # amplitude cutoff ≤ 10%
}

def apply_thresholds(metrics_df, thresholds):
    mask = pd.Series(True, index=metrics_df.index)
    for metric, (low, high) in thresholds.items():
        if metric not in metrics_df.columns:
            continue
        if low is not None:
            mask &= metrics_df[metric] >= low
        if high is not None:
            mask &