Pro Neuropixels Analysis

Analyze high-density neural recordings from Neuropixels probes using SpikeInterface and established neuroscience analysis workflows. This skill covers spike sorting, quality metrics, unit curation, probe visualization, and publication-ready neural data analysis.

When to Use This Skill

Choose Pro Neuropixels Analysis when you need to:

• Perform spike sorting on Neuropixels probe recordings
• Compute quality metrics for sorted neural units (ISI violations, SNR, drift)
• Visualize probe geometry, spike waveforms, and neural activity
• Curate spike sorting results with automated and manual quality control

Consider alternatives when:

• You need real-time neural signal processing (use Open Ephys or BrainFlow)
• You need calcium imaging analysis (use CaImAn or Suite2p)
• You need LFP spectral analysis only (use MNE-Python or Chronux)

Quick Start

Copy
# Install SpikeInterface with all sorters
pip install spikeinterface[full]

Copy
import spikeinterface.full as si
from pathlib import Path

# Load Neuropixels recording
recording = si.read_spikeglx("/path/to/recording", stream_id="imec0.ap")
print(f"Channels: {recording.get_num_channels()}")
print(f"Duration: {recording.get_total_duration():.1f} seconds")
print(f"Sampling rate: {recording.get_sampling_frequency()} Hz")

# Preprocessing
recording_filtered = si.bandpass_filter(recording, freq_min=300, freq_max=6000)
recording_cmr = si.common_reference(recording_filtered, reference="global", operator="median")

# Run spike sorting with Kilosort
sorting = si.run_sorter(
    "kilosort3",
    recording_cmr,
    output_folder="./kilosort3_output",
    verbose=True
)
print(f"Found {len(sorting.get_unit_ids())} units")

Core Concepts

Analysis Pipeline Steps

Step	Function	Description
Load data	si.read_spikeglx()	Import Neuropixels recordings
Filter	si.bandpass_filter()	300-6000 Hz bandpass for spikes
Re-reference	si.common_reference()	Remove common noise across channels
Spike sort	si.run_sorter()	Detect and cluster spikes
Quality metrics	si.compute_quality_metrics()	SNR, ISI violations, drift
Curate	Manual + auto thresholds	Keep only high-quality units
Export	si.export_to_phy()	Phy format for manual curation

Quality Metrics and Curation

Copy
import spikeinterface.full as si

def curate_sorting(recording, sorting, output_dir):
    """Compute quality metrics and auto-curate sorted units."""

    # Extract waveforms
    we = si.extract_waveforms(
        recording, sorting,
        folder=f"{output_dir}/waveforms",
        max_spikes_per_unit=500,
        n_jobs=4
    )

    # Compute quality metrics
    metrics = si.compute_quality_metrics(we, metric_names=[
        "snr", "isi_violation", "firing_rate",
        "presence_ratio", "amplitude_cutoff",
        "drift", "nearest_neighbor"
    ])

    print(metrics.describe())

    # Auto-curation thresholds
    good_units = metrics[
        (metrics["snr"] > 5) &
        (metrics["isi_violations_ratio"] < 0.01) &
        (metrics["firing_rate"] > 0.1) &
        (metrics["presence_ratio"] > 0.9) &
        (metrics["amplitude_cutoff"] < 0.1)
    ].index.tolist()

    print(f"Good units: {len(good_units)}/{len(sorting.get_unit_ids())}")

    # Create curated sorting
    curated = sorting.select_units(good_units)
    return curated, metrics

curated, metrics = curate_sorting(recording_cmr, sorting, "./results")

Visualization

Copy
import spikeinterface.full as si
import matplotlib.pyplot as plt

def visualize_units(we, unit_ids=None, n_units=10):
    """Generate publication-quality unit visualizations."""
    if unit_ids is None:
        unit_ids = we.sorting.get_unit_ids()[:n_units]

    fig, axes = plt.subplots(len(unit_ids), 3, figsize=(15, 4 * len(unit_ids)))

    for i, unit_id in enumerate(unit_ids):
        # Waveform
        template = we.get_template(unit_id)
        best_chan = template.ptp(axis=0).argmax()
        axes[i, 0].plot(template[:, best_chan], "k-", linewidth=1.5)
        axes[i, 0].set_title(f"Unit {unit_id} - Waveform")
        axes[i, 0].set_ylabel("µV")

        # ISI histogram
        spike_train = we.sorting.get_unit_spike_train(unit_id)
        isis = np.diff(spike_train) / we.sampling_frequency * 1000  # ms
        axes[i, 1].hist(isis[isis < 50], bins=100, color="steelblue")
        axes[i, 1].axvline(1.0, color="red", linestyle="--", label="1ms")
        axes[i, 1].set_title("ISI Distribution")
        axes[i, 1].set_xlabel("ISI (ms)")

        # Firing rate over time
        spike_times = spike_train / we.sampling_frequency
        bins = np.arange(0, spike_times.max(), 10)  # 10s bins
        rate, _ = np.histogram(spike_times, bins)
        rate = rate / 10.0  # Convert to Hz
        axes[i, 2].plot(bins[:-1], rate, color="darkgreen")
        axes[i, 2].set_title("Firing Rate")
        axes[i, 2].set_xlabel("Time (s)")
        axes[i, 2].set_ylabel("Hz")

    fig.tight_layout()
    fig.savefig("unit_summary.png", dpi=200)

import numpy as np
visualize_units(we, n_units=5)

Configuration

Parameter	Description	Default
freq_min	Highpass filter cutoff (Hz)	300
freq_max	Lowpass filter cutoff (Hz)	6000
sorter	Spike sorting algorithm	"kilosort3"
n_jobs	Parallel processing workers	4
max_spikes_per_unit	Waveforms to extract per unit	500
snr_threshold	Minimum SNR for good units	5.0
isi_threshold_ms	ISI violation window (ms)	1.0

Best Practices

1. Apply common median reference before sorting — Neuropixels recordings contain shared noise across channels. Common median referencing (reference="global", operator="median") removes this noise more effectively than per-channel bandpass filtering alone, significantly improving sorting quality.

2. Use multiple quality metrics for curation — Don't rely on a single metric. Combine SNR (signal quality), ISI violations (sorting accuracy), firing rate (physiological plausibility), presence ratio (recording stability), and amplitude cutoff (complete capture). Units should pass all criteria.

3. Inspect drift before interpreting results — Neuropixels probes often drift during long recordings. Use si.compute_drift() to quantify motion and apply drift correction before spike sorting. Uncorrected drift causes single neurons to be split across multiple units.

4. Export to Phy for final manual curation — Automated quality metrics provide a good starting point, but always perform manual curation with Phy for publication-quality data. Inspect waveform shapes, autocorrelograms, and cross-correlograms to merge over-split units and remove artifacts.

5. Save intermediate results at each pipeline stage — Spike sorting and waveform extraction are computationally expensive (hours for long recordings). Save preprocessed recordings and sorting results to disk using SpikeInterface's built-in save functions so you can resume from any stage without recomputing.

Common Issues

Kilosort not found or fails to start — Kilosort3 requires MATLAB runtime or a compiled binary. Use si.installed_sorters() to check available sorters. For a MATLAB-free workflow, use Kilosort4 (Python-native) or SpykingCircus2: si.run_sorter("spykingcircus2", recording).

Too many units with high ISI violations — ISI violation rates above 1-2% indicate sorting errors where spikes from different neurons were assigned to the same unit. Try increasing the number of clusters in the sorter parameters, or switch to a different sorting algorithm. Also verify that your bandpass filter settings aren't distorting spike shapes.

Memory errors during waveform extraction — Extracting waveforms from all spikes across 384 channels requires significant RAM. Reduce max_spikes_per_unit to 300-500, use n_jobs=1 to reduce parallel memory usage, or set dtype=np.float32 instead of float64 to halve memory requirements.