neuropixels-analysis

Neuropixels 分析

概述

Neuropixels 探针以 30 kHz 的采样率，同时记录 384 个（NP 1.0）或 192 个（NP 2.0）通道的细胞外电压。分析遵循标准流程：原始数据 → 尖峰分选（Kilosort）→ 质量筛选 → 单元分析。SpikeInterface（Python）为 10 多种尖峰分选器提供统一 API，支持从多种格式（SpikeGLX、OpenEphys、NWB）加载数据、计算质量指标并导出分选结果。ProbeInterface 管理探针几何结构和通道映射。分选后的分析（PSTH、发放率、解码）使用标准 Python 科学计算工具栈。

适用场景

• 对 SpikeGLX（.bin）或 OpenEphys（.dat）记录的 Neuropixels 数据进行尖峰分选，以提取单单元活动
• 应用自动质量控制指标（ISI 违规、SNR、发放率）筛选分选后的单元
• 计算与实验事件对齐的刺激后时间直方图（PSTH）
• 分析群体编码：根据发放率解码刺激或行为变量
• 将分选数据转换为 NWB（Neurodata Without Borders）格式以便共享
• 在同一数据集上比较多种尖峰分选器，以验证方法
• 可视化单元波形、自相关图以及探针通道上的空间分布
• 如需使用支持 10 多种尖峰分选器并提供统一 API 和比较工具的框架，请改用 SpikeInterface

前置条件

• Python 包：spikeinterface[full]、probeinterface、numpy、pandas、matplotlib、scipy
• 尖峰分选器：Kilosort4（Python，pip install kilosort）；或基于 MATLAB 的 Kilosort2/3（需要 MATLAB 许可证）
• 数据要求：原始 .bin（SpikeGLX）或 .dat（OpenEphys）记录；探针通道映射文件（.prb 或来自 ProbeInterface）
• 硬件：Kilosort 需要 NVIDIA GPU（4 GB 以上 VRAM）；也可使用 CPU，但速度约慢 10 倍

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

工作流

步骤 1：加载原始记录

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()}")

步骤 2：应用公共参考和带通滤波

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")

步骤 3：运行尖峰分选

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]]}")

步骤 4：计算波形和质量指标

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())

步骤 5：筛选单元

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 &= metrics_df[metric] <= high
    return mask

good_units_mask = apply_thresholds(metrics, thresholds)
good_unit_ids = metrics[good_units_mask].index.tolist()

print(f"Total units:   {len(metrics)}")
print(f"Good units:    {len(good_unit_ids)} ({100*len(good_unit_ids)/len(metrics):.0f}%)")

# Filter sorting to good units only
sorting_curated = sorting.select_units(good_unit_ids)

步骤 6：计算 PSTH 并进行可视化

import numpy as np
import matplotlib.pyplot as plt

def compute_psth(spike_times, event_times, window=(-0.5, 1.0), bin_size=0.01, fs=30000):
    """Compute peri-stimulus time histogram for one unit."""
    bins = np.arange(window[0], window[1] + bin_size, bin_size)
    spike_times_s = spike_times / fs  # samples → seconds
    counts = np.zeros(len(bins) - 1)
    for t_event in event_times:
        rel_times = spike_times_s - t_event
        in_window = rel_times[(rel_times >= window[0]) & (rel_times < window[1])]
        counts += np.histogram(in_window, bins=bins)[0]
    rate = counts / (len(event_times) * bin_size)  # convert to Hz
    return bins[:-1] + bin_size / 2, rate  # bin centers, firing rate

# Example: visual stimulus events at 1.0, 2.5, 4.0, 5.5 s
fs = 30000
event_times_s = np.array([1.0, 2.5, 4.0, 5.5, 7.0, 8.5])

# Plot PSTH for first 4 good units
fig, axes = plt.subplots(2, 2, figsize=(10, 6))
for ax, unit_id in zip(axes.flat, good_unit_ids[:4]):
    spikes = sorting_curated.get_unit_spike_train(unit_id, segment_index=0)
    times, rate = compute_psth(spikes, event_times_s, fs=fs)
    ax.bar(times, rate, width=0.01, color="steelblue", alpha=0.8)
    ax.axvline(0, color="red", lw=1.5, linestyle="--", label="Stimulus")
    ax.set_xlabel("Time from stimulus (s)")
    ax.set_ylabel("Firing rate (Hz)")
    ax.set_title(f"Unit {unit_id}")

plt.tight_layout()
plt.savefig("psth_grid.pdf", bbox_inches="tight")
print("Saved psth_grid.pdf")

步骤 7：导出为 NWB

import spikeinterface.exporters as sexp

# Export sorted data to Neurodata Without Borders format
nwb_path = "./recording_sorted.nwb"
sexp.export_to_nwb(
    sorting_curated,
    nwb_file_path=nwb_path,
    overwrite=True,
)
print(f"Exported to NWB: {nwb_path}")

# Also export waveforms for Phy manual curation
phy_dir = "./phy_export"
sexp.export_to_phy(
    waveforms,
    output_folder=phy_dir,
    compute_pc_features=True,
    copy_binary=True,
)
print(f"Phy export ready at: {phy_dir}")
print("Launch with: phy template-gui phy_export/params.py")

关键参数

参数	模块/函数	默认值	范围 / 选项	效果
freq_min / freq_max	bandpass_filter	300/6000 Hz	150–500 / 3000–10000	尖峰频段；NP 1.0 的标准范围为 300–6000 Hz
nblocks	Kilosort4	5	0–10	漂移校正块数；0 表示禁用漂移校正
Th_learned	Kilosort4	8	6–12	检测阈值（× 噪声水平）；越低 = 尖峰越多，噪声也越多
ms_before / ms_after	extract_waveforms	1.0/2.0	0.5–2.0/1.0–3.0 ms	尖峰峰值附近的波形片段时间窗
max_spikes_per_unit	extract_waveforms	500	100–5000	每个单元用于提取波形的最大尖峰数
snr 阈值	质量指标	—	5–10	用于接收单元的信噪比阈值
isi_violations_ratio	质量指标	—	0.05–0.2	不应期违规率；<0.1 表示分离良好
presence_ratio	质量指标	—	0.9	单元处于活动状态的记录时长比例
bin_size	PSTH	0.01 s	0.001–0.05 s	PSTH 时间分辨率；越小 = 细节越多、噪声越大

关键概念

尖峰分选流程

原始电压 → 预处理（公共参考、带通滤波）→ 检测（阈值交叉或学习得到的模板）→ 聚类（PCA + k-means 或高斯混合）→ 模板匹配 → 质量指标 → 筛选。Kilosort4 增加了漂移校正，这对慢性 NP 记录至关重要（数小时内探针尖端漂移 5–50 µm 会降低单元分离效果）。

质量指标

• SNR：波形峰值幅度 / 噪声水平。SNR >5 表示该单元隔离良好
• ISI 违例率：短于不应期（约 1.5 ms）的脉冲间隔所占比例。比率 <0.1 表示单单元隔离
• 存在率：至少发放过一个脉冲的记录时段所占比例。<0.9 表示该单元消失或丢失
• 幅度截断率：因检测阈值而漏检的脉冲所占比例。<0.1 可确保近乎完整的检测

常用示例

示例：比较同一段记录上的两种排序器

import spikeinterface.sorters as ss
import spikeinterface.comparison as sc

# Run two sorters
sorting_ks4 = ss.run_sorter("kilosort4", recording_clean, output_folder="./ks4")
sorting_sc   = ss.run_sorter("spykingcircus2", recording_clean, output_folder="./sc2")

# Compare outputs
comparison = sc.compare_two_sorters(sorting_ks4, sorting_sc,
                                     sorting1_name="Kilosort4",
                                     sorting2_name="SpykingCircus2")
print(comparison.get_performance())
# Shows: agreement fraction, recall, precision per matched unit pair

示例：群体发放率热图

import numpy as np
import matplotlib.pyplot as plt

fs = 30000
duration_s = sorting_curated.get_total_duration()
bin_edges = np.arange(0, duration_s, 0.05)  # 50 ms bins

# Build firing rate matrix: (n_units, n_bins)
fr_matrix = np.zeros((len(good_unit_ids), len(bin_edges)-1))
for i, uid in enumerate(good_unit_ids[:50]):
    spikes_s = sorting_curated.get_unit_spike_train(uid, segment_index=0) / fs
    counts, _ = np.histogram(spikes_s, bins=bin_edges)
    fr_matrix[i] = counts / 0.05  # Hz

# Normalize each unit
fr_norm = (fr_matrix - fr_matrix.mean(1, keepdims=True)) / (fr_matrix.std(1, keepdims=True) + 1e-6)

fig, ax = plt.subplots(figsize=(12, 5))
im = ax.imshow(fr_norm, aspect="auto", cmap="RdBu_r",
               extent=[0, duration_s, len(good_unit_ids[:50]), 0], vmin=-2, vmax=2)
ax.set_xlabel("Time (s)")
ax.set_ylabel("Unit #")
ax.set_title("Population Activity Heatmap (z-scored FR)")
plt.colorbar(im, ax=ax, label="z-score")
plt.tight_layout()
plt.savefig("population_heatmap.pdf", bbox_inches="tight")

预期输出

输出	格式	典型内容
sorting/	SpikeInterface 文件夹	每个单元的脉冲时间；模板波形
waveforms/	SpikeInterface 文件夹	波形片段 (n_spikes, n_samples, n_channels)
quality_metrics.csv	CSV	每个单元的 SNR、ISI 违例、发放率、存在率
psth_grid.pdf	PDF	精选单元的 PSTH 图
recording_sorted.nwb	NWB/HDF5	用于共享的可移植脉冲排序数据
phy_export/	Phy 格式	用于手动整理 GUI 的 .npy + params.py

故障排除

问题	原因	解决方案
Kilosort GPU 内存错误	记录时长过长，或通道数量超出 VRAM 容量	使用 si.SubRecording 按时间段处理记录；使用 NP 2.0（192 ch）设置
未找到任何单元	阈值过高，或预处理移除了所有信号	将 Th_learned 降至 6；确认 recording_clean 中的通道数据非零
ISI 违例过多（>50%）	多单元活动被合并为单个单元	使用 Phy 手动整理；提高 Th_learned，使其更加保守
波形提取时内存不足	脉冲数 × 通道数 × 采样数过多	减少 max_spikes_per_unit；提高 n_jobs 以并行提取
漂移校正失败	脉冲过少或记录时间过短	设置 nblocks=1（最小校正）或 nblocks=0（禁用校正）
NWB 导出失败	缺少会话元数据（受试者、日期）	提供 NWBFile 元数据；或使用 sexp.export_to_nwb(..., metadata={...})
找不到 SpikeGLX 文件	流名称错误	列出流：si.get_neo_streams("spikeglx", data_dir)

参考资料

• SpikeInterface 文档 — 完整 API、教程和排序器比较指南
• Kilosort4 论文：Pachitariu 等（2024），Nature Methods — 漂移校正和模板学习
• SpikeInterface 论文：Buccino 等（2020），eLife — 用于细胞外电生理学的统一框架
• ProbeInterface 文档 — 探针几何结构和通道映射处理
• Neuropixels 文档 — 硬件规格、imec 参考设计
• NWB 文档 — 用于归档和共享的神经生理学数据标准