Neuropixels Analysis

Source: https://github.com/aipoch/medical-research-skills

When to Use

Use this skill in any of the following situations:

1. You need to load and standardize Neuropixels recordings from SpikeGLX (.ap.bin/.lf.bin/.meta), Open Ephys (.continuous/.oebin), or NWB (.nwb) into a consistent analysis pipeline.
2. You are preparing raw extracellular data for spike sorting, including high-pass filtering, phase shift correction (NP1.0), bad channel detection/removal, and common average referencing (CAR).
3. You suspect probe drift or tissue motion and need to estimate and correct motion before sorting (especially when drift is > ~10 µm).
4. You want to run spike sorting (Kilosort4 recommended; CPU alternatives supported) and then compute post-processing products (waveforms, templates, amplitudes, correlograms, unit locations).
5. You need quality control and curation using Allen/IBL-style thresholds, plus optional AI-assisted visual review for borderline units, and exports to Phy/NWB.

Key Features

• Multi-format ingestion: SpikeGLX, Open Ephys, and NWB readers via SpikeInterface.
• Neuropixels-aware preprocessing:
  • High-pass filtering for spike band
  • Phase shift correction for Neuropixels 1.0
  • Bad channel detection and removal
  • Median CAR / referencing
• Motion/drift workflow:
  • Motion estimation presets (e.g., “Kilosort-like”)
  • Optional rigid/non-rigid correction presets
  • Drift visualization outputs
• Spike sorting orchestration:
  • Kilosort4 (GPU) recommended
  • CPU alternatives (e.g., SpykingCircus2, Mountainsort5, Tridesclous2)
• Post-processing and QC:
  • SortingAnalyzer-based computation of waveforms, templates, amplitudes, correlograms, unit locations, and quality metrics
• Curation:
  • Allen/IBL-style automated labeling
  • Optional AI-assisted visual analysis for uncertain units
• Reporting and export:
  • HTML report generation
  • Export to Phy and NWB
  • Save metrics tables (CSV)

Reference guides (if present in the repository) can be used for deeper explanations:

• reference/standard_workflow.md
• reference/api_reference.md
• reference/plotting_guide.md
• reference/PREPROCESSING.md, reference/MOTION_CORRECTION.md, reference/SPIKE_SORTING.md
• reference/QUALITY_METRICS.md, reference/AUTOMATED_CURATION.md, reference/AI_CURATION.md

Dependencies

Python dependencies (typical versions known to work; adjust to your environment):

• python >= 3.9
• spikeinterface[full] >= 0.99
• probeinterface >= 0.2
• neo >= 0.13
• Spike sorters (optional, depending on what you run):
  • kilosort >= 4.0 (Kilosort4; GPU required)
  • spykingcircus >= 1.1 (SpykingCircus2; CPU)
  • mountainsort5 >= 0.5 (CPU)
• Optional (AI-assisted curation):
  • anthropic >= 0.20
• Optional (IBL tooling):
  • ibllib >= 2.0
  • ibl-neuropixel >= 1.0

Example Usage

The following example is designed to be a complete, runnable script (assuming dependencies and a valid dataset path). It loads SpikeGLX data, preprocesses, estimates/corrects motion, runs Kilosort4, computes metrics, curates units, generates a report, and exports to Phy and NWB.

import spikeinterface.full as si
import neuropixels_analysis as npa

def main():
    # Parallelization / chunking settings used by SpikeInterface functions
    job_kwargs = dict(n_jobs=-1, chunk_duration="1s", progress_bar=True)

    # 1) Load data (SpikeGLX example)
    # For Open Ephys: si.read_openephys("/path/to/Record_Node_101/")
    # For NWB:        si.read_nwb("/path/to/file.nwb")
    recording = si.read_spikeglx("/path/to/spikeglx_folder", stream_id="imec0.ap")

    # Optional: slice first 60 seconds for a quick test
    fs = recording.get_sampling_frequency()
    recording = recording.frame_slice(0, int(60 * fs))

    # 2) Preprocess (recommended chain; wrapper may include the same steps)
    # Note: phase_shift is mandatory for Neuropixels 1.0 and not needed for 2.0.
    rec = npa.preprocess(recording)

    # 3) Estimate drift/motion and correct if needed
    motion_info = npa.estimate_motion(rec, preset="kilosort_like", **job_kwargs)
    npa.plot_drift(rec, motion_info, output="drift_map.png")

    # Example threshold: correct if max drift exceeds 10 µm
    if float(motion_info["motion"].max()) > 10.0:
        rec = npa.correct_motion(rec, preset="nonrigid_accurate", **job_kwargs)

    # 4) Spike sorting (Kilosort4 recommended; requires GPU)
    sorting = si.run_sorter("kilosort4", rec, folder="ks4_output", **job_kwargs)

    # 5) Post-processing + metrics
    analyzer = si.create_sorting_analyzer(sorting, rec, sparse=True)

    analyzer.compute("random_spikes", max_spikes_per_unit=500, **job_kwargs)
    analyzer.compute("waveforms", ms_before=1.0, ms_after=2.0, **job_kwargs)
    analyzer.compute("templates", operators=["average", "std"], **job_kwargs)
    analyzer.compute("spike_amplitudes", **job_kwargs)
    analyzer.compute("correlograms", window_ms=50.0, bin_ms=1.0, **job_kwargs)
    analyzer.compute("unit_locations", method="monopolar_triangulation", **job_kwargs)
    analyzer.compute("quality_metrics", **job_kwargs)

    metrics = analyzer.get_extension("quality_metrics").get_data()
    metrics.to_csv("quality_metrics.csv")

    # 6) Automated curation (Allen/IBL-style)
    labels = npa.curate(metrics, method="allen")  # e.g., "allen", "ibl", "strict"

    # 7) Report
    results = {"sorting": sorting, "metrics": metrics, "labels": labels, "analyzer": analyzer}
    npa.generate_analysis_report(results, "output_report/")
    npa.print_analysis_summary(results)

    # 8) Export
    si.export_to_phy(
        analyzer,
        output_folder="phy_export/",
        compute_pc_features=True,
        compute_amplitudes=True,
    )

    from spikeinterface.exporters import export_to_nwb
    export_to_nwb(rec, sorting, "output.nwb")

if __name__ == "__main__":
    main()

Implementation Details

1) Data I/O and supported formats

• SpikeGLX: si.read_spikeglx(path, stream_id="imec0.ap")
• Open Ephys: si.read_openephys(path)
• NWB: si.read_nwb(path)

Neuropixels probe types commonly encountered:

• Neuropixels 1.0: requires phase shift correction to align channels.
• Neuropixels 2.0: denser geometries; phase shift correction typically not required.

2) Preprocessing chain (typical)

A standard spike-band preprocessing sequence is:

1. High-pass filter (commonly 300–400 Hz) to isolate spikes.
2. Phase shift correction (si.phase_shift) for NP1.0.
3. Bad channel detection (si.detect_bad_channels) and removal.
4. Common reference (often median CAR) to reduce shared noise.

Key parameters:

• freq_min (high-pass cutoff): typical 300–400 Hz
• bad channel detection sensitivity (implementation-dependent; often exposed as thresholds/presets)

3) Motion estimation and correction

• Motion/drift can strongly degrade sorting quality; a practical rule is to inspect drift before sorting.
• Presets:
  • preset="kilosort_like": faster estimation aligned with common sorter assumptions
  • preset="nonrigid_accurate": more robust correction for severe drift

Operational threshold often used in practice:

• If estimated drift exceeds ~10 µm, apply correction before sorting.

4) Spike sorting

• Kilosort4 is recommended for Neuropixels due to speed and quality, but requires a GPU.
• CPU alternatives can be used when GPU is unavailable (at the cost of runtime and sometimes quality).

Sorter parameters to tune (Kilosort4 examples):

• batch_size: samples per batch (often ~30000 by default)
• nblocks: number of drift blocks (increase for long recordings)
• Th_learned: detection threshold (lower → more spikes, potentially more false positives)

5) Post-processing and quality metrics

Using SortingAnalyzer, the pipeline typically computes:

• waveforms (window: ms_before, ms_after)
• templates (average/std)
• spike amplitudes
• correlograms (e.g., window_ms=50, bin_ms=1)
• unit locations (e.g., monopolar_triangulation)
• quality metrics (e.g., SNR, ISI violations, presence ratio, amplitude cutoff)

Common QC thresholds (dataset-dependent; document your choices):

• snr_threshold: often 3–5
• isi_violations_ratio: often 0.01–0.5
• presence_ratio: often 0.5–0.95

6) Curation logic (Allen/IBL-style)

A conservative “good unit” selection often combines:

• high presence ratio (stable across recording)
• low ISI violations (refractory period respected)
• low amplitude cutoff (less truncation / missed spikes)

Example rule (illustrative):

• presence_ratio > 0.9
• isi_violations_ratio < 0.5
• amplitude_cutoff < 0.1

7) AI-assisted visual analysis (optional)

For borderline units (e.g., moderate SNR), AI-assisted review can be used to interpret:

• waveform shape consistency
• refractory period evidence in autocorrelograms
• amplitude stability and drift effects
• multi-unit contamination indicators

If your repository provides npa.analyze_unit_visually(...), it can be integrated with an API client (e.g., anthropic) to generate structured curation suggestions.