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）の記録をスパイクソーティングして、単一ユニット活動を抽出する
• 自動品質管理指標（ISI 違反、SNR、発火率）を適用して、ソート済みユニットをキュレーションする
• 実験イベントに合わせた刺激周辺時間ヒストグラム（PSTH）を計算する
• 集団コーディングを解析する：発火率から刺激変数または行動変数をデコードする
• ソート済みデータを共有用に NWB（Neurodata Without Borders）形式に変換する
• 手法を検証するため、同じデータセットで複数のスパイクソーターを比較する
• ユニットの波形、自己相関ヒストグラム、プローブチャネル上の空間分布を可視化する
• 共通 API と比較ツールを備え、10 種類以上のスパイクソーターに対応する統一フレームワークが必要な場合は、代わりに 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（VRAM 4 GB 以上）が必要です。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 は、単一ユニットが分離されていることを示します
• 存在率：ユニットが少なくとも 1 回スパイクを発火した記録エポックの割合です。<0.9 は、ユニットが消失または失われた可能性を示します
• 振幅カットオフ：検出しきい値によって検出されなかったスパイクの割合です。<0.1 であれば、ほぼ完全な検出が保証されます

よく使うレシピ

レシピ：同じ記録で 2 つのソーターを比較する

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）の設定を使用します
ユニットが 0 個見つかる	しきい値が高すぎる、または前処理ですべての信号が除去された	Th_learned を 6 に下げます。recording_clean のチャンネルデータがゼロでないことを確認します
ISI 違反が過剰（>50%）	マルチユニット活動が単一ユニットとして統合された	Phy で手動キュレーションを行います。より保守的になるよう Th_learned を上げます
波形抽出で OOM が発生する	スパイク数 × チャンネル数 × サンプル数が多すぎる	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 et al. (2024), Nature Methods — ドリフト補正とテンプレート学習
• SpikeInterface の論文：Buccino et al. (2020), eLife — 細胞外電気生理学の統合フレームワーク
• ProbeInterface のドキュメント — プローブのジオメトリとチャンネルマップの処理
• Neuropixels のドキュメント — ハードウェア仕様、imec リファレンス設計
• NWB のドキュメント — アーカイブと共有のための神経生理学データ標準