spikeinterface-electrophysiology

SpikeInterface — 统一的细胞外电生理框架

概述

SpikeInterface 提供通用的 Python API，可读取 20 多种文件格式的细胞外记录、预处理原始电压轨迹、运行 10 多种 spike sorter、对分选后的单元进行后处理和质量控制，并导出结果，全程无需编写特定格式的代码。其模块化设计让用户可以更换 sorter、格式和预处理步骤，无需重写工作流。SpikeInterface 基于延迟执行、可链式调用的对象构建：Recording 保存原始数据，Sorting 保存 spike 时间，SortingAnalyzer 则将二者关联起来，用于计算波形和指标。

适用场景

• 使用统一 API 加载来自多种采集系统（SpikeGLX、OpenEphys、Intan、NWB、Maxwell MEA、Blackrock）的记录，而不是使用特定格式的解析器
• 在不同硬件上记录的实验中运行相同的预处理和分选工作流
• 对同一段记录比较两个或更多 spike sorter，评估一致性并选择最佳结果
• 通过 Docker 或 Singularity 运行容器化 sorter（Kilosort、IronClust、MountainSort5），无需本地安装
• 计算标准质量指标（SNR、ISI 违规、放电率、存在率、振幅截断），并根据阈值进行筛选
• 使用合成或混合的真实值记录验证 spike 分选准确性
• 将分选结果导出为 NWB 以共享数据，或导出到 Phy 进行手动筛选
• 如果需要完整的、针对 Neuropixels 的 Kilosort4 工作流（包括 PSTH 计算、调谐曲线和群体解码），请改用 neuropixels-analysis
• 如果要处理 EEG、ECG 或其他生物信号（而非 spike 分选），请改用 neurokit2

前置条件

• Python 包：spikeinterface[full]>=0.101、probeinterface、numpy、matplotlib
• 可选 sorter 依赖：kilosort（pip），或用于容器化 sorter 的 Docker/Singularity
• 数据要求：原始二进制记录文件，以及探针几何信息（.prb、.json，或从格式中自动检测）
• 硬件：Kilosort4 需要 GPU；其他 sorter 均可在 CPU 上运行

pip install "spikeinterface[full]>=0.101" probeinterface
# Optional: Kilosort4 Python package
pip install kilosort
# Optional: Phy for manual curation
pip install phy

快速入门

import spikeinterface.full as si
import spikeinterface.preprocessing as spre
import spikeinterface.sorters as ss
import spikeinterface.qualitymetrics as sqm

# Load, preprocess, sort, and inspect quality metrics in 10 lines
recording = si.read_openephys("/data/session_001", stream_name="Signals CH")
recording_pp = spre.bandpass_filter(
    spre.common_reference(recording, reference="global", operator="median"),
    freq_min=300, freq_max=6000,
)
sorting = ss.run_sorter("spykingcircus2", recording_pp, output_folder="./sc2_out")
analyzer = si.create_sorting_analyzer(sorting, recording_pp, folder="./analyzer")
analyzer.compute(["random_spikes", "waveforms", "templates", "noise_levels"])
metrics = sqm.compute_quality_metrics(analyzer, metric_names=["snr", "firing_rate", "isi_violation"])
print(metrics.describe())

核心 API

模块 1：记录 I/O

SpikeInterface 将所有采集格式封装在统一的 BaseRecording 接口之后。加载后，无论源格式为何，所有对象都提供相同的方法。

import spikeinterface.full as si

# SpikeGLX (.bin + .meta)
recording_sglx = si.read_spikeglx("/data/session_001", stream_name="imec0.ap")

# OpenEphys (binary or classic)
recording_oe = si.read_openephys("/data/oe_session", stream_name="Signals CH")

# NWB file
recording_nwb = si.read_nwb_recording("/data/recording.nwb",
                                       electrical_series_name="ElectricalSeries")

# Intan RHD/RHS
recording_intan = si.read_intan("/data/session.rhd", stream_name="RHn")

# Inspect any recording with the same API
print(f"Format:       {type(recording_sglx).__name__}")
print(f"Channels:     {recording_sglx.get_num_channels()}")
print(f"Sampling rate:{recording_sglx.get_sampling_frequency()} Hz")
print(f"Duration:     {recording_sglx.get_total_duration():.1f} s")
print(f"Probe:        {recording_sglx.get_probe().name}")

# List available streams before loading (useful when a file has multiple streams)
streams = si.get_neo_streams("spikeglx", "/data/session_001")
print("Available streams:", streams)
# e.g. ['imec0.ap', 'imec0.lf', 'nidq']

# Select a time slice (lazy, no data loaded until get_traces() is called)
recording_slice = recording_sglx.frame_slice(
    start_frame=0,
    end_frame=int(60 * recording_sglx.get_sampling_frequency()),  # first 60 s
)
print(f"Sliced duration: {recording_slice.get_total_duration():.1f} s")

模块 2：预处理

预处理函数会返回新的 Recording 对象，用于封装原始对象；读取数据时才会延迟应用处理链。即使处理数小时的记录，也能保持较低的内存占用。

import spikeinterface.preprocessing as spre

# 1. Common median reference — removes shared noise across all channels
recording_cmr = spre.common_reference(recording_sglx,
                                       reference="global",
                                       operator="median")

# 2. Bandpass filter for action potentials (300–6000 Hz typical)
recording_filt = spre.bandpass_filter(recording_cmr,
                                       freq_min=300,
                                       freq_max=6000)

# 3. Remove bad channels automatically (coherence-based detection)
recording_clean, removed_ids = spre.remove_bad_channels(recording_filt,
                                                          method="coherence+psd")
print(f"Removed {len(removed_ids)} bad channels: {removed_ids}")
print(f"Clean channels: {recording_clean.get_num_channels()}")

# Whitening — decorrelates channels; recommended before template-matching sorters
recording_white = spre.whiten(recording_clean, mode="local")

# Phase shift correction for Neuropixels (samples acquired with small time offsets)
recording_shifted = spre.phase_shift(recording_clean)

# Inspect a short snippet of preprocessed data
traces = recording_white.get_traces(start_frame=0, end_frame=3000, segment_index=0)
print(f"Trace snippet shape: {traces.shape}")   # (3000, n_channels)
print(f"Trace range: [{traces.min():.2f}, {traces.max():.2f}] µV")

模块 3：Spike 分选

ss.run_sorter() 以统一的调用签名封装所有受支持的 sorter。Sorter 特定参数通过关键字参数传入；其他所有工作流步骤都相同。

import spikeinterface.sorters as ss
from pathlib import Path

# List all sorters available in the current environment
available = ss.available_sorters()
print("Available sorters:", available)

# List sorters that can run without local installation (via container)
installed = ss.installed_sorters()
print("Installed locally:", installed)

# Run SpykingCircus2 (CPU, no external deps)
sorting_sc2 = ss.run_sorter(
    "spykingcircus2",
    recording_clean,
    output_folder=Path("./sorter_output/sc2"),
    remove_existing_folder=True,
    verbose=True,
)
print(f"SpykingCircus2 units: {len(sorting_sc2.get_unit_ids())}")

# Run Kilosort4 via Docker container (no local GPU/MATLAB required)
sorting_ks4 = ss.run_sorter(
    "kilosort4",
    recording_clean,
    output_folder=Path("./sorter_output/ks4"),
    singularity_image=False,   # use Docker; set True for Singularity
    docker_image=True,
    remove_existing_folder=True,
    # Kilosort4-specific parameters
    nblocks=5,
    Th_learned=8,
    do_correction=True,
)
print(f"Kilosort4 units: {len(sorting_ks4.get_unit_ids())}")

# Run MountainSort5 (CPU, fast, good for tetrode/low-channel-count probes)
sorting_ms5 = ss.run_sorter(
    "mountainsort5",
    recording_clean,
    output_folder=Path("./sorter_output/ms5"),
    scheme="2",          # scheme 2 is recommended for high-density probes
    detect_threshold=5.5,
)
print(f"MountainSort5 units: {len(sorting_ms5.get_unit_ids())}")

模块 4：后处理（SortingAnalyzer）

SortingAnalyzer 是 SpikeInterface >= 0.101 中的核心后处理对象。它取代了旧版 WaveformExtractor，并为波形、模板、PCA 和下游指标提供统一接口。

import spikeinterface.full as si
import spikeinterface.postprocessing as spost

# Create analyzer (saves to disk; use format="memory" for in-RAM only)
analyzer = si.create_sorting_analyzer(
    sorting_sc2,
    recording_clean,
    folder="./analyzer_sc2",
    format="binary_folder",
    overwrite=True,
    sparse=True,           # sparse=True: only nearby channels per unit
    ms_before=1.0,
    ms_after=2.0,
)

# Compute extensions in dependency order
analyzer.compute([
    "random_spikes",       # subsample spike indices for waveform extraction
    "waveforms",           # raw waveform snippets per unit
    "templates",           # mean/std template per unit
    "noise_levels",        # per-channel noise estimate
])

# Retrieve templates
templates = analyzer.get_extension("templates").get_data(outputs="Templates")
print(f"Templates object: {templates}")
print(f"Unit 0 template shape: {templates.get_one_template_dense(0).shape}")
# (n_samples, n_channels)

# Compute amplitude and PCA extensions (needed for quality metrics)
analyzer.compute([
    "spike_amplitudes",          # amplitude at peak channel per spike
    "principal_components",      # PCA scores (n_components x n_spikes)
    "template_similarity",       # pairwise template correlation matrix
    "correlograms",              # auto- and cross-correlograms
    "unit_locations",            # estimated unit position on probe (center of mass)
])

# Access spike amplitudes for first unit
ext_amp = analyzer.get_extension("spike_amplitudes")
unit_ids = analyzer.unit_ids
amps = ext_amp.get_data()[analyzer.sorting.ids_to_indices([unit_ids[0]])]
print(f"Unit {unit_ids[0]} — median amplitude: {abs(amps).median():.1f} µV")

模块 5：质量指标

质量指标用于概括单元分离质量。只需要 spike 时间的指标（ISI 违规、放电率）计算很快；需要波形的指标（SNR、振幅截断）则需要先填充 SortingAnalyzer。

import spikeinterface.qualitymetrics as sqm

# Compute a standard panel of quality metrics
metrics = sqm.compute_quality_metrics(
    analyzer,
    metric_names=[
        "snr",                    # signal-to-noise ratio of template peak
        "isi_violation",          # fraction of ISIs < refractory period
        "firing_rate",            # mean firing rate (Hz) over recording
        "presence_ratio",         # fraction of time windows with ≥1 spike
        "amplitude_cutoff",       # estimated fraction of spikes below threshold
        "nearest_neighbor",       # isolation distance in PCA space
        "silhouette_score",       # cluster separation in PCA space
    ],
)
print(metrics.head())
print(f"\nShape: {metrics.shape}")  # (n_units, n_metrics)

import pandas as pd

# Apply threshold-based curation (Allen Brain Institute defaults)
thresholds = {
    "snr":                   (">=", 5.0),
    "isi_violations_ratio":  ("<=", 0.1),
    "firing_rate":           (">=", 0.1),
    "presence_ratio":        (">=", 0.9),
    "amplitude_cutoff":      ("<=", 0.1),
}

keep = pd.Series(True, index=metrics.index)
for col, (op, val) in thresholds.items():
    if col not in metrics.columns:
        continue
    if op == ">=":
        keep &= metrics[col] >= val
    else:
        keep &= metrics[col] <= val

good_unit_ids = metrics[keep].index.tolist()
print(f"Total units:   {len(metrics)}")
print(f"Curated units: {len(good_unit_ids)} ({100*len(good_unit_ids)/len(metrics):.0f}%)")

# Filter analyzer to good units
sorting_curated = sorting_sc2.select_units(good_unit_ids)

模块 6：比较与导出

将不同 sorter 相互比较，或与真实值比较，然后以便于共享的格式导出结果。

import spikeinterface.comparison as sc

# Compare two sorters — matches units by spike train overlap
comparison = sc.compare_two_sorters(
    sorting_sc2,
    sorting_ks4,
    sorting1_name="SpykingCircus2",
    sorting2_name="Kilosort4",
    match_score=0.5,        # minimum overlap to count as a match
    delta_time=0.4,         # coincidence window (ms)
)

# Performance summary per matched unit pair
perf = comparison.get_performance(method="by_unit")
print(perf.head(10))
# Columns: accuracy, recall, precision, false_discovery_rate, miss_rate

# Agreement score matrix (fraction overlap between all unit pairs)
agreement_matrix = comparison.get_agreement_fraction_table()
print(f"Agreement matrix shape: {agreement_matrix.shape}")

import spikeinterface.exporters as sexp

# Export curated sorting to NWB (Neurodata Without Borders)
sexp.export_to_nwb(
    sorting_curated,
    nwb_file_path="./session_sorted.nwb",
    overwrite=True,
)
print("Exported to NWB: session_sorted.nwb")

# Export to Phy for manual curation
sexp.export_to_phy(
    analyzer,
    output_folder="./phy_export",
    compute_pc_features=True,
    copy_binary=True,
    remove_if_exists=True,
)
print("Phy export ready at: ./phy_export")
print("Launch Phy with: phy template-gui phy_export/params.py")

常见工作流

工作流 1：比较 OpenEphys 数据上的多个 Sorter

目标：加载 OpenEphys 记录、进行预处理、运行两个 sorter、比较二者的一致性、筛选产出较高的结果，并导出为 NWB。

import spikeinterface.full as si
import spikeinterface.preprocessing as spre
import spikeinterface.sorters as ss
import spikeinterface.comparison as sc
import spikeinterface.qualitymetrics as sqm
import spikeinterface.exporters as sexp
from pathlib import Path

# --- Step 1: Load ---
data_dir = Path("/data/oe_recording")
streams = si.get_neo_streams("openephys", data_dir)
print("Streams:", streams)

recording = si.read_openephys(data_dir, stream_name="Signals CH")
print(f"Loaded: {recording.get_num_channels()} ch, "
      f"{recording.get_sampling_frequency()} Hz, "
      f"{recording.get_total_duration():.1f} s")

# --- Step 2: Preprocess ---
rec = spre.bandpass_filter(recording, freq_min=300, freq_max=6000)
rec = spre.common_reference(rec, reference="global", operator="median")
rec, bad_ids = spre.remove_bad_channels(rec, method="coherence+psd")
print(f"Preprocessing complete. Removed channels: {bad_ids}")

# --- Step 3: Run two sorters ---
out = Path("./sorting_outputs")
sorting_sc2 = ss.run_sorter("spykingcircus2", rec,
                              output_folder=out / "sc2",
                              remove_existing_folder=True)
sorting_tdc = ss.run_sorter("tridesclous2", rec,
                              output_folder=out / "tdc",
                              remove_existing_folder=True)
print(f"SC2 units: {len(sorting_sc2.unit_ids)}, "
      f"TDC units: {len(sorting_tdc.unit_ids)}")

# --- Step 4: Compare ---
cmp = sc.compare_two_sorters(sorting_sc2, sorting_tdc,
                               sorting1_name="SC2",
                               sorting2_name="Tridesclous2",
                               match_score=0.5)
perf = cmp.get_performance(method="pooled_with_average")
print(f"\nAgreement performance:\n{perf}")

# --- Step 5: Quality metrics on SC2 (higher yield) ---
analyzer = si.create_sorting_analyzer(sorting_sc2, rec,
                                        folder="./analyzer_sc2",
                                        overwrite=True, sparse=True)
analyzer.compute(["random_spikes", "waveforms", "templates",
                  "noise_levels", "spike_amplitudes"])
metrics = sqm.compute_quality_metrics(
    analyzer,
    metric_names=["snr", "firing_rate", "isi_violation",
                  "presence_ratio", "amplitude_cutoff"],
)

keep = (metrics["snr"] >= 5) & (metrics["isi_violations_ratio"] <= 0.1) \
     & (metrics["firing_rate"] >= 0.1) & (metrics["presence_ratio"] >= 0.9)
sorting_curated = sorting_sc2.select_units(metrics[keep].index.tolist())
print(f"\nCurated: {len(sorting_curated.unit_ids)} / {len(sorting_sc2.unit_ids)} units")

# --- Step 6: Export winner to NWB ---
sexp.export_to_nwb(sorting_curated,
                    nwb_file_path="./session_sorted.nwb",
                    overwrite=True)
print("Saved: session_sorted.nwb")

工作流 2：使用合成录制进行真实值验证

目标：生成包含已知尖峰序列的合成录制，运行排序器，并根据真实值测量实际准确度（召回率、精确率）——用于对排序器进行基准测试或测试预处理流水线。

import spikeinterface.full as si
import spikeinterface.preprocessing as spre
import spikeinterface.sorters as ss
import spikeinterface.comparison as sc
import numpy as np

# --- Step 1: Generate ground-truth synthetic recording ---
# Uses a Marsaglia noise model with realistic waveform templates
recording_gt, sorting_gt = si.generate_ground_truth_recording(
    durations=[120.0],             # 120 s recording
    sampling_frequency=30000.0,
    num_channels=32,
    num_units=10,
    noise_kwargs={"noise_level": 10.0, "dtype": "float32"},
    seed=42,
)
print(f"GT recording: {recording_gt.get_num_channels()} ch, "
      f"{recording_gt.get_total_duration():.0f} s")
print(f"GT units: {len(sorting_gt.unit_ids)}")
print(f"GT firing rates: "
      f"{[round(len(sorting_gt.get_unit_spike_train(u, 0))/120, 1) for u in sorting_gt.unit_ids]} Hz")

# --- Step 2: Preprocess ---
rec_pp = spre.bandpass_filter(recording_gt, freq_min=300, freq_max=6000)
rec_pp = spre.common_reference(rec_pp, reference="global", operator="median")

# --- Step 3: Sort with two sorters ---
sorting_sc2 = ss.run_sorter("spykingcircus2", rec_pp,
                              output_folder="./gt_sc2",
                              remove_existing_folder=True)
sorting_ms5 = ss.run_sorter("mountainsort5", rec_pp,
                              output_folder="./gt_ms5",
                              remove_existing_folder=True,
                              scheme="2")

# --- Step 4: Compare each sorter against ground truth ---
for name, sorting_test in [("SC2", sorting_sc2), ("MS5", sorting_ms5)]:
    cmp = sc.compare_sorter_to_ground_truth(sorting_gt, sorting_test,
                                              exhaustive_gt=True)
    perf = cmp.get_performance(method="pooled_with_average")
    print(f"\n{name} vs Ground Truth:")
    print(f"  Accuracy:  {perf['accuracy']:.3f}")
    print(f"  Recall:    {perf['recall']:.3f}")
    print(f"  Precision: {perf['precision']:.3f}")
    print(f"  Well-detected units: {cmp.get_well_detected_units(well_detected_score=0.8)}")

工作流 3：批量处理多个会话

目标：将相同的预处理 + 排序流水线应用于多个录制会话，并收集所有会话的质量指标。

import spikeinterface.full as si
import spikeinterface.preprocessing as spre
import spikeinterface.sorters as ss
import spikeinterface.qualitymetrics as sqm
import pandas as pd
from pathlib import Path

sessions = list(Path("/data/experiment").glob("session_*/"))
all_metrics = []

for session_dir in sessions:
    print(f"Processing {session_dir.name} ...")
    try:
        streams = si.get_neo_streams("spikeglx", session_dir)
        ap_stream = [s for s in streams if "ap" in s][0]
        rec = si.read_spikeglx(session_dir, stream_name=ap_stream)

        # Preprocess
        rec = spre.bandpass_filter(
            spre.common_reference(rec, reference="global", operator="median"),
            freq_min=300, freq_max=6000,
        )
        rec, _ = spre.remove_bad_channels(rec)

        # Sort
        out_dir = session_dir / "sorting"
        sorting = ss.run_sorter("spykingcircus2", rec,
                                 output_folder=out_dir,
                                 remove_existing_folder=True)

        # Compute metrics
        analyzer = si.create_sorting_analyzer(
            sorting, rec, folder=session_dir / "analyzer", overwrite=True, sparse=True
        )
        analyzer.compute(["random_spikes", "waveforms", "templates",
                          "noise_levels", "spike_amplitudes"])
        m = sqm.compute_quality_metrics(
            analyzer, metric_names=["snr", "firing_rate", "isi_violation"]
        )
        m["session"] = session_dir.name
        all_metrics.append(m)

    except Exception as e:
        print(f"  FAILED: {e}")
        continue

# Combine across sessions
combined = pd.concat(all_metrics)
combined.to_csv("all_sessions_metrics.csv")
print(f"\nSaved metrics: {combined.shape[0]} units across {len(all_metrics)} sessions")
print(combined.groupby("session")[["snr", "firing_rate"]].median())

关键参数

参数	模块 / 函数	默认值	范围 / 选项	影响
freq_min / freq_max	spre.bandpass_filter	300 / 6000 Hz	150–500 / 3000–10000 Hz	尖峰频段；AP 活动使用 300–6000 Hz
reference	spre.common_reference	"global"	"global"、"local"、"single"	用于中位数参考减法的通道子集
method	spre.remove_bad_channels	"coherence+psd"	"coherence+psd"、"std"、"mad"	检测不良通道的算法
scheme	ss.run_sorter("mountainsort5")	"2"	"1"、"2"、"3"	排序方案；高密度探针推荐使用方案 2
nblocks	ss.run_sorter("kilosort4")	5	0–10	漂移校正块数；0 表示禁用漂移校正
Th_learned	ss.run_sorter("kilosort4")	8	6–12	检测阈值（× 噪声）；值越低 = 单元越多、噪声越多
match_score	sc.compare_two_sorters	0.5	0.1–0.9	判定单元匹配所需的最小尖峰序列重叠度
sparse	si.create_sorting_analyzer	True	True、False	将波形提取限制在每个单元附近的通道；减少内存占用
ms_before / ms_after	si.create_sorting_analyzer	1.0 / 2.0 ms	0.5–2.0 / 1.0–3.0 ms	相对于检测到的尖峰峰值的波形片段窗口
snr threshold	sqm.compute_quality_metrics	—	推荐 5–10	振幅 / 噪声比；> 5 表示单元隔离良好
isi_violations_ratio	sqm.compute_quality_metrics	—	推荐 ≤ 0.1	ISI 小于不应期（1.5 ms）的比例；< 0.1 = 单单元
presence_ratio	sqm.compute_quality_metrics	—	推荐 ≥ 0.9	单元在其中放电的录制时段比例；< 0.9 = 漂移单元

最佳实践

1. 加载前始终检查可用的数据流：不同的采集系统会将 AP 数据、LFP 数据和辅助通道保存为单独的数据流。加载错误的数据流会悄无声息地得到看似有效但实际不正确的数据。

   streams = si.get_neo_streams("spikeglx", data_dir)
   print(streams)  # e.g. ['imec0.ap', 'imec0.lf', 'nidq']
   recording = si.read_spikeglx(data_dir, stream_name="imec0.ap")
   

2. 按惰性方式串联预处理；不要过早加载到内存：预处理对象采用惰性执行，会在读取时应用转换。在预处理前对原始记录调用 get_traces()，会不必要地将未过滤的数据加载到 RAM。请先构建完整的处理链，再访问数据。

3. 创建 SortingAnalyzer 时使用 sparse=True：对于通道数较多的探针（64–384 个通道），密集波形提取所需的 RAM 和磁盘空间会比稀疏提取高出 10–50 倍。稀疏模式只会在每个单元附近的通道上提取波形。

4. 运行容器化排序器以避免依赖冲突：Kilosort2/3（MATLAB）、IronClust 及其他排序器的依赖关系较为复杂。使用 docker_image=True 运行 run_sorter()，即可拉取官方容器并在隔离环境中运行排序器：

   sorting = ss.run_sorter("kilosort2_5", recording_clean,
                            output_folder="./ks25_out",
                            docker_image=True)
   

5. 按依赖顺序计算指标扩展：扩展之间存在依赖关系。标准顺序为：random_spikes → waveforms → templates → noise_levels → spike_amplitudes → principal_components。跳过较早的步骤会导致请求后续步骤时出现 MissingExtensionError。

6. 将大型记录的 SortingAnalyzer 保存到磁盘：进程退出后，内存中的分析器（format="memory"）将会丢失。对于时长超过 30 分钟或包含许多单元的记录，请务必指定 folder 路径，以便重新加载分析器：

   analyzer = si.load_sorting_analyzer("./analyzer_sc2")
   

7. 不要使用不匹配的预处理结果比较排序器：对排序器进行基准测试时，请让所有排序器都使用同一个经过预处理的 recording_clean 对象。在不同预处理链上运行排序器会使比较结果无效。

常见操作示例

示例：加载并检查多数据流记录

适用场景：快速检查不熟悉的记录中有哪些数据流，并在开始完整排序前确认通道数和时长。

import spikeinterface.full as si

data_dir = "/data/recording_session"

# Try SpikeGLX first; if it fails, try OpenEphys
try:
    streams = si.get_neo_streams("spikeglx", data_dir)
    fmt = "spikeglx"
except Exception:
    streams = si.get_neo_streams("openephys", data_dir)
    fmt = "openephys"

print(f"Format: {fmt}")
print(f"Streams: {streams}")

for stream in streams:
    try:
        rec = si.read_spikeglx(data_dir, stream_name=stream) if fmt == "spikeglx" \
              else si.read_openephys(data_dir, stream_name=stream)
        print(f"  {stream}: {rec.get_num_channels()} ch, "
              f"{rec.get_sampling_frequency()} Hz, "
              f"{rec.get_total_duration():.1f} s")
    except Exception as e:
        print(f"  {stream}: could not load ({e})")

示例：将质量指标报告导出为 CSV

适用场景：计算质量指标后，保存一个整洁的 CSV 文件，总结所有单元的指标及通过/未通过列，以便进行后续分析或与协作者分享。

import spikeinterface.qualitymetrics as sqm
import pandas as pd

metrics = sqm.compute_quality_metrics(
    analyzer,
    metric_names=["snr", "firing_rate", "isi_violation",
                  "presence_ratio", "amplitude_cutoff"],
)

# Add pass/fail column based on standard thresholds
metrics["pass_qc"] = (
    (metrics["snr"] >= 5) &
    (metrics["isi_violations_ratio"] <= 0.1) &
    (metrics["firing_rate"] >= 0.1) &
    (metrics["presence_ratio"] >= 0.9) &
    (metrics["amplitude_cutoff"] <= 0.1)
)

metrics.to_csv("unit_quality_metrics.csv")
n_pass = metrics["pass_qc"].sum()
print(f"QC report saved: {len(metrics)} total units, {n_pass} pass ({100*n_pass/len(metrics):.0f}%)")
print(metrics[metrics["pass_qc"]].describe())

示例：探针几何结构可视化

适用场景：排序前确认探针通道映射已正确加载。通道映射错误会悄无声息地降低高密度探针的排序质量。

import spikeinterface.full as si
import matplotlib.pyplot as plt
import probeinterface.plotting as pp

recording = si.read_spikeglx("/data/session_001", stream_name="imec0.ap")
probe = recording.get_probe()
print(f"Probe name: {probe.name}")
print(f"N contacts: {probe.get_contact_count()}")
print(f"Contact positions (first 5):\n{probe.contact_positions[:5]}")

fig, ax = plt.subplots(figsize=(3, 10))
pp.plot_probe(probe, ax=ax, with_channel_index=True)
ax.set_title(f"{probe.name} — channel map")
plt.tight_layout()
plt.savefig("probe_geometry.png", dpi=150)
print("Saved probe_geometry.png")

故障排查

问题	原因	解决方案
ValueError: stream_name not found	录制包含多个流，但未指定任何流	运行 si.get_neo_streams(format, path) 列出可用流；将正确的流传递给读取器
分选器输出的单元数为 0	检测阈值过高，或预处理移除了所有信号	确认 recording_clean.get_traces() 返回非零数据；降低检测阈值（例如，Kilosort4 使用 Th_learned=6）
MissingExtensionError	分析器扩展依赖尚未计算的前置项	按照规范的计算顺序执行：random_spikes → waveforms → templates → noise_levels → spike_amplitudes
Docker 分选器在启动时卡住	Docker 守护进程未运行，或尚未拉取镜像	运行 docker ps 确认 Docker 正在运行；使用 docker pull spikeinterface/kilosort4-compiled-base 手动拉取镜像
提取波形时出现 MemoryError	在高通道数探针上进行密集提取	在 create_sorting_analyzer 中使用 sparse=True；减小 max_spikes_per_unit（默认值为 500）
坏通道检测移除了过多通道	阈值过于激进，或录制时长太短	设置 method="std" 使用更简单的阈值；增大 bad_threshold 参数
单元比较显示分选器之间的一致率为 0%	时间差窗口过窄，或匹配分数要求过严	增大 delta_time（默认值为 0.4 ms），并降低 match_score（可尝试 0.3）
NWB 导出时，单元属性引发 TypeError	分选结果包含来自分选器的不可序列化属性	导出前移除有问题的属性：sorting.remove_unit_property("property_name")
read_spikeglx 无法读取 LF 流	LFP 流使用了不同的文件后缀（.lf.bin）	显式指定 stream_name="imec0.lf"；使用 ls data_dir/*.lf.bin 确认文件存在

相关技能

• neuropixels-analysis — 面向啮齿动物和灵长类动物实验的 Neuropixels 专用工作流，使用 SpikeInterface + Kilosort4，并提供 PSTH、调谐曲线和群体解码
• neurokit2 — 用于生物信号处理（ECG、EEG、EDA、EMG、PPG），而非尖峰分选；数据并非细胞外电生理数据时使用

参考资料

• SpikeInterface 文档 — 完整的 API 参考、教程和分选器专用指南
• SpikeInterface GitHub — 源代码、变更日志和问题跟踪器
• Buccino 等人（2020），eLife — SpikeInterface 论文 — 统一框架设计以及不同分选器的基准测试
• ProbeInterface 文档 — 探针几何信息处理和通道映射格式
• SpikeInterface 分选器列表 — 支持的分选器、要求和容器镜像
• NWB 文档 — 用于共享神经生理数据的 Neurodata Without Borders 格式