Skills › Scientific Agent Skills › neuropixels-analysis › preprocess_recording.py
preprocess_recording.py
135 lines · skills/neuropixels-analysis/scripts/preprocess_recording.py at 5f2bc03
1#!/usr/bin/env python
2"""
3Preprocess Neuropixels recording.
5Usage:
6 python preprocess_recording.py /path/to/data --output preprocessed/ --format spikeglx
7"""
9import argparse
10from pathlib import Path
12import spikeinterface.full as si
13from _common import validate_recording, apply_phase_correction, reference_by_shank
16def preprocess_recording(
17 input_path: str,
18 output_dir: str,
19 format: str = 'auto',
20 stream_name: str = None,
21 freq_min: float = 300,
22 freq_max: float = 6000,
23 phase_shift: bool = True,
24 common_ref: bool = True,
25 detect_bad: bool = True,
26 n_jobs: int = -1,
27):
28 """Preprocess a Neuropixels recording."""
30 print(f"Loading recording from: {input_path}")
32 # Identify acquisition format explicitly; do not hide reader errors as fallback loads.
33 source = Path(input_path)
34 if format == 'auto':
35 if source.suffix.lower() == '.nwb':
36 format = 'nwb'
37 elif any(source.rglob('*.ap.meta')):
38 format = 'spikeglx'
39 elif any(source.rglob('*.oebin')):
40 format = 'openephys'
41 else:
42 format = 'saved'
43 if format == 'spikeglx':
44 recording = si.read_spikeglx(input_path, stream_name=stream_name or 'imec0.ap')
45 elif format == 'openephys':
46 recording = si.read_openephys(input_path, stream_name=stream_name)
47 elif format == 'nwb':
48 recording = si.read_nwb(input_path)
49 elif format == 'saved':
50 recording = si.load(input_path)
51 else:
52 raise ValueError(f"Unknown format: {format}")
54 validate_recording(recording)
55 if not 0 < freq_min < freq_max < recording.get_sampling_frequency() / 2:
56 raise ValueError("Require 0 < freq_min < freq_max < Nyquist frequency.")
58 print(f"Recording: {recording.get_num_channels()} channels, {recording.get_total_duration():.1f}s")
60 # Preprocessing chain
61 rec = recording
63 # Bandpass filter
64 print(f"Applying bandpass filter ({freq_min}-{freq_max} Hz)...")
65 rec = si.bandpass_filter(rec, freq_min=freq_min, freq_max=freq_max)
67 # Phase shift correction (for Neuropixels ADC)
68 if phase_shift:
69 print("Applying phase shift correction...")
70 rec = apply_phase_correction(rec)
72 # Bad channel detection
73 if detect_bad:
74 print("Detecting bad channels...")
75 bad_channel_ids, bad_labels = si.detect_bad_channels(rec)
76 if len(bad_channel_ids) > 0:
77 print(f" Removing {len(bad_channel_ids)} bad channels: {bad_channel_ids[:10]}...")
78 rec = rec.remove_channels(bad_channel_ids)
80 # Common median reference
81 if common_ref:
82 print("Applying common median reference...")
83 rec = reference_by_shank(rec)
85 # Save preprocessed
86 output_path = Path(output_dir)
87 output_path.mkdir(parents=True, exist_ok=True)
89 print(f"Saving preprocessed recording to: {output_path}")
90 rec = rec.save(folder=output_path / 'preprocessed', n_jobs=n_jobs)
92 # Save probe info
93 probe = rec.get_probe()
94 if probe is not None:
95 from probeinterface import write_probeinterface
96 write_probeinterface(output_path / 'probe.json', probe)
98 print("Done!")
99 print(f" Output channels: {rec.get_num_channels()}")
100 print(f" Output duration: {rec.get_total_duration():.1f}s")
102 return rec
105def main():
106 parser = argparse.ArgumentParser(description='Preprocess Neuropixels recording')
107 parser.add_argument('input', help='Path to input recording')
108 parser.add_argument('--output', '-o', default='preprocessed/', help='Output directory')
109 parser.add_argument('--format', '-f', default='auto', choices=['auto', 'spikeglx', 'openephys', 'nwb'])
110 parser.add_argument('--stream-name', default=None, help='Stream name for multi-probe recordings (e.g. imec0.ap)')
111 parser.add_argument('--freq-min', type=float, default=300, help='Highpass cutoff (Hz)')
112 parser.add_argument('--freq-max', type=float, default=6000, help='Lowpass cutoff (Hz)')
113 parser.add_argument('--no-phase-shift', action='store_true', help='Skip phase shift correction')
114 parser.add_argument('--no-cmr', action='store_true', help='Skip common median reference')
115 parser.add_argument('--no-bad-channel', action='store_true', help='Skip bad channel detection')
116 parser.add_argument('--n-jobs', type=int, default=-1, help='Number of parallel jobs')
118 args = parser.parse_args()
120 preprocess_recording(
121 args.input,
122 args.output,
123 format=args.format,
124 stream_name=args.stream_name,
125 freq_min=args.freq_min,
126 freq_max=args.freq_max,
127 phase_shift=not args.no_phase_shift,
128 common_ref=not args.no_cmr,
129 detect_bad=not args.no_bad_channel,
130 n_jobs=args.n_jobs,
131 )
134if __name__ == '__main__':
135 main()