spikeinterface-electrophysiology

SpikeInterface — 細胞外電気生理学の統合フレームワーク

概要

SpikeInterface は、20 種類以上のファイル形式から細胞外記録を読み込み、生の電圧トレースを前処理し、10 種類以上のスパイクソーターを実行し、ソート済みユニットを後処理して品質管理を行い、結果をエクスポートするための共通 Python API を提供します。形式ごとのコードは不要です。モジュール設計により、パイプラインを書き直すことなく、ソーター、形式、前処理ステップを切り替えられます。SpikeInterface は遅延評価に対応した連結可能なオブジェクトを中心に構成されています。Recording は生データを保持し、Sorting はスパイク時刻を保持し、SortingAnalyzer は波形とメトリクスの計算のために両者を結び付けます。

使用する場面

• 形式ごとのパーサーではなく、統一 API を使って複数の記録システム（SpikeGLX、OpenEphys、Intan、NWB、Maxwell MEA、Blackrock）から記録を読み込む場合
• 異なるハードウェアで記録された実験に対して、同じ前処理とソーティングのパイプラインを実行する場合
• 同じ記録で 2 種類以上のスパイクソーターを比較し、一致度を評価して最適な出力を選ぶ場合
• ローカルにインストールせず、Docker または Singularity を介してコンテナ化されたソーター（Kilosort、IronClust、MountainSort5）を実行する場合
• 標準的な品質メトリクス（SNR、ISI 違反、発火率、存在率、振幅カットオフ）を計算し、しきい値に基づくキュレーションを適用する場合
• 合成またはハイブリッドのグラウンドトゥルース記録を使って、スパイクソーティングの精度を検証する場合
• データ共有用にソート済みの結果を NWB にエクスポートする場合、または手動キュレーション用に Phy にエクスポートする場合
• PSTH 計算、チューニング曲線、集団デコーディングを含む、Neuropixels に特化した完全な Kilosort4 ワークフローには、代わりに neuropixels-analysis を使用してください
• EEG、ECG、またはその他の生体信号の処理（スパイクソーティング以外）には、代わりに neurokit2 を使用してください

前提条件

• Python パッケージ：spikeinterface[full]>=0.101、probeinterface、numpy、matplotlib
• オプションのソーター依存関係：kilosort（pip）、またはコンテナ化されたソーター用の Docker/Singularity
• データ要件：生のバイナリー記録ファイルとプローブのジオメトリ（.prb、.json、または形式から自動検出）
• ハードウェア：Kilosort4 には GPU が必要です。それ以外のソーターはすべて 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：記録の入出力

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：スパイクソーティング

ss.run_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：品質メトリクス

品質メトリクスは、ユニットの分離品質を要約します。スパイク時刻のみを必要とするメトリクス（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：比較とエクスポート

ソーター同士、またはグラウンドトゥルースとソーターを比較し、結果を共有可能な形式でエクスポートします。

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 データで複数のソーターを比較する

目的：OpenEphys の記録を読み込み、前処理を行い、2 種類のソーターを実行して一致度を比較し、より多くのユニットが得られた出力をキュレーションして、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 しきい値	sqm.compute_quality_metrics	—	5–10 を推奨	振幅／ノイズ比です。5 を超える場合、ユニットが十分に分離されていることを示します
isi_violations_ratio	sqm.compute_quality_metrics	—	≤ 0.1 を推奨	不応期（1.5 ms）未満の ISI の割合です。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 などのソーターには複雑な依存関係があります。run_sorter() で docker_image=True を使用すると、公式コンテナを取得して、ソーターを分離環境で実行できます：

   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) を実行して利用可能なストリームを一覧表示し、正しいストリームをリーダーに渡します
ソーターの出力ユニット数がゼロです	検出しきい値が高すぎるか、前処理ですべての信号が除去されています	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 — PSTH、チューニング曲線、集団デコーディングを含む、げっ歯類および霊長類の実験向け Neuropixels 専用パイプライン。SpikeInterface と Kilosort4 を使用します
• neurokit2 — スパイクソーティングではなく、生体信号（ECG、EEG、EDA、EMG、PPG）の処理向けです。データが細胞外電気生理データではない場合に使用します

参考資料

• SpikeInterface のドキュメント — API リファレンス、チュートリアル、ソーター別ガイド一式
• SpikeInterface の GitHub — ソースコード、変更履歴、課題トラッカー
• Buccino et al. (2020), eLife — SpikeInterface の論文 — ソーター間で統一されたフレームワークの設計とベンチマーク
• ProbeInterface のドキュメント — プローブのジオメトリ処理とチャンネルマップ形式
• SpikeInterface のソーター一覧 — 対応ソーター、要件、コンテナイメージ
• NWB のドキュメント — 神経生理データ共有のための Neurodata Without Borders 形式