Skills › Scientific Agent Skills › neuropixels-analysis › analysis_template.py
analysis_template.py
271 lines · skills/neuropixels-analysis/assets/analysis_template.py at 5f2bc03
1#!/usr/bin/env python
2"""
3Neuropixels Analysis Template
5Complete analysis workflow from raw data to curated units.
6Illustrative template: copy, customize, and validate on your own recording and sorter.
8Usage:
9 1. Copy this file to your analysis directory
10 2. Copy the skill scripts/ alongside it as neuropixels_scripts/
11 3. Update the PARAMETERS section
12 4. Run: python analysis_template.py
13"""
15# =============================================================================
16# PARAMETERS - Customize these for your analysis
17# =============================================================================
19# Input/Output paths
20DATA_PATH = '/path/to/your/spikeglx/data/'
21OUTPUT_DIR = 'analysis_output/'
22DATA_FORMAT = 'spikeglx' # 'spikeglx', 'openephys', or 'nwb'
23STREAM_NAME = 'imec0.ap' # Select the exact AP stream discovered in acquisition metadata
25# Preprocessing parameters
26FREQ_MIN = 300 # Highpass filter (Hz)
27FREQ_MAX = 6000 # Lowpass filter (Hz)
28APPLY_PHASE_SHIFT = True
29APPLY_CMR = True
30DETECT_BAD_CHANNELS = True
32# Motion correction
33CORRECT_MOTION = True
34MOTION_PRESET = 'nonrigid_accurate' # 'kilosort_like', 'nonrigid_fast_and_accurate'
36# Spike sorting
37SORTER = 'kilosort4' # 'kilosort4', 'spykingcircus2', 'mountainsort5'
38SORTER_PARAMS = {
39 'batch_size': 30000,
40 'nblocks': 1, # Internal Kilosort correction only when CORRECT_MOTION=False
41}
43# Quality metrics and curation
44CURATION_METHOD = 'allen' # 'allen', 'ibl', 'strict'
46# Processing
47N_JOBS = -1 # -1 = all cores
49# =============================================================================
50# ANALYSIS PIPELINE - Usually no need to modify below
51# =============================================================================
53from pathlib import Path
54import json
55import sys
57# When used inside this skill, use ../scripts; when copied, use the bundled copy.
58script_dir = Path(__file__).resolve().parents[1] / "scripts"
59if not (script_dir / "compute_metrics.py").is_file():
60 script_dir = Path(__file__).resolve().parent / "neuropixels_scripts"
61sys.path.insert(0, str(script_dir))
62from compute_metrics import curate_units, METRIC_NAMES
63from _common import validate_recording, apply_phase_correction, reference_by_shank
65import spikeinterface.full as si
66from spikeinterface.exporters import export_to_phy
69def main():
70 """Run the full analysis pipeline."""
72 output_path = Path(OUTPUT_DIR)
73 output_path.mkdir(parents=True, exist_ok=True)
75 # =========================================================================
76 # 1. LOAD DATA
77 # =========================================================================
78 print("=" * 60)
79 print("1. LOADING DATA")
80 print("=" * 60)
82 if DATA_FORMAT == 'spikeglx':
83 recording = si.read_spikeglx(DATA_PATH, stream_name=STREAM_NAME)
84 elif DATA_FORMAT == 'openephys':
85 recording = si.read_openephys(DATA_PATH, stream_name=STREAM_NAME)
86 elif DATA_FORMAT == 'nwb':
87 recording = si.read_nwb(DATA_PATH)
88 else:
89 raise ValueError(f"Unknown format: {DATA_FORMAT}")
91 validate_recording(recording)
92 si.set_global_job_kwargs(n_jobs=N_JOBS, chunk_duration="1s")
94 print(f"Recording: {recording.get_num_channels()} channels")
95 print(f"Duration: {recording.get_total_duration():.1f} seconds")
96 print(f"Sampling rate: {recording.get_sampling_frequency()} Hz")
98 # =========================================================================
99 # 2. PREPROCESSING
100 # =========================================================================
101 print("\n" + "=" * 60)
102 print("2. PREPROCESSING")
103 print("=" * 60)
105 rec = recording
107 # Bandpass filter
108 print(f"Applying bandpass filter ({FREQ_MIN}-{FREQ_MAX} Hz)...")
109 rec = si.bandpass_filter(rec, freq_min=FREQ_MIN, freq_max=FREQ_MAX)
111 # Phase shift correction
112 if APPLY_PHASE_SHIFT:
113 print("Applying phase shift correction...")
114 rec = apply_phase_correction(rec)
116 # Bad channel detection
117 if DETECT_BAD_CHANNELS:
118 print("Detecting bad channels...")
119 bad_ids, _ = si.detect_bad_channels(rec)
120 if len(bad_ids) > 0:
121 print(f" Removing {len(bad_ids)} bad channels")
122 rec = rec.remove_channels(bad_ids)
124 # Common median reference
125 if APPLY_CMR:
126 print("Applying common median reference...")
127 rec = reference_by_shank(rec)
129 # Save preprocessed
130 print("Saving preprocessed recording...")
131 rec = rec.save(folder=output_path / 'preprocessed', n_jobs=N_JOBS)
133 # =========================================================================
134 # 3. MOTION CORRECTION
135 # =========================================================================
136 if CORRECT_MOTION:
137 print("\n" + "=" * 60)
138 print("3. MOTION CORRECTION")
139 print("=" * 60)
141 print(f"Estimating and correcting motion (preset: {MOTION_PRESET})...")
142 rec = si.correct_motion(
143 rec,
144 preset=MOTION_PRESET,
145 folder=output_path / 'motion',
146 )
148 # =========================================================================
149 # 4. SPIKE SORTING
150 # =========================================================================
151 print("\n" + "=" * 60)
152 print("4. SPIKE SORTING")
153 print("=" * 60)
155 print(f"Running {SORTER}...")
156 sorter_params = SORTER_PARAMS.copy()
157 if CORRECT_MOTION:
158 # Disable the sorter's own motion stage after external interpolation.
159 if SORTER in {"kilosort2_5", "kilosort3", "kilosort4"}:
160 sorter_params["do_correction"] = False
161 elif SORTER == "spykingcircus2":
162 sorter_params["apply_motion_correction"] = False
163 sorting = si.run_sorter(
164 SORTER,
165 rec,
166 folder=output_path / f'{SORTER}_output',
167 verbose=True,
168 **sorter_params,
169 )
171 print(f"Found {len(sorting.unit_ids)} units")
173 # =========================================================================
174 # 5. POSTPROCESSING
175 # =========================================================================
176 print("\n" + "=" * 60)
177 print("5. POSTPROCESSING")
178 print("=" * 60)
180 print("Creating SortingAnalyzer...")
181 analyzer = si.create_sorting_analyzer(
182 sorting,
183 rec,
184 format='binary_folder',
185 folder=output_path / 'analyzer',
186 sparse=True,
187 )
189 print("Computing extensions...")
190 analyzer.compute('random_spikes', max_spikes_per_unit=500)
191 analyzer.compute('waveforms', ms_before=1.0, ms_after=2.0)
192 analyzer.compute('templates', operators=['average', 'std'])
193 analyzer.compute('noise_levels')
194 analyzer.compute('spike_amplitudes')
195 analyzer.compute('amplitude_scalings')
196 analyzer.compute('correlograms', window_ms=50.0, bin_ms=1.0)
197 analyzer.compute('unit_locations', method='monopolar_triangulation')
199 # =========================================================================
200 # 6. QUALITY METRICS
201 # =========================================================================
202 print("\n" + "=" * 60)
203 print("6. QUALITY METRICS")
204 print("=" * 60)
206 print("Computing quality metrics...")
207 metrics = si.compute_quality_metrics(
208 analyzer,
209 metric_names=METRIC_NAMES,
210 n_jobs=N_JOBS,
211 )
213 metrics.to_csv(output_path / 'quality_metrics.csv')
214 print(f"Saved metrics to: {output_path / 'quality_metrics.csv'}")
216 # Print summary
217 print("\nMetrics summary:")
218 for col in ['snr', 'isi_violations_ratio', 'presence_ratio', 'firing_rate']:
219 if col in metrics.columns:
220 print(f" {col}: {metrics[col].median():.4f} (median)")
222 # =========================================================================
223 # 7. CURATION
224 # =========================================================================
225 print("\n" + "=" * 60)
226 print("7. CURATION")
227 print("=" * 60)
229 labels = curate_units(metrics, method=CURATION_METHOD)
231 # Save labels
232 with open(output_path / 'curation_labels.json', 'w') as f:
233 json.dump({str(uid): label for uid, label in labels.items()}, f, indent=2)
235 # Count
236 good_count = sum(1 for v in labels.values() if v == 'good')
237 mua_count = sum(1 for v in labels.values() if v == 'mua')
238 noise_count = sum(1 for v in labels.values() if v == 'noise')
240 print(f"\nCuration results:")
241 print(f" Good: {good_count}")
242 print(f" MUA: {mua_count}")
243 print(f" Noise: {noise_count}")
244 print(f" Total: {len(labels)}")
246 # =========================================================================
247 # 8. EXPORT
248 # =========================================================================
249 print("\n" + "=" * 60)
250 print("8. EXPORT")
251 print("=" * 60)
253 analyzer.sorting.set_property('curation_label', [labels[uid] for uid in analyzer.sorting.unit_ids])
254 good_ids = [uid for uid, label in labels.items() if label == 'good']
255 sorting.select_units(good_ids).save(folder=output_path / 'sorting_curated')
256 print("Exporting all units with screening labels to Phy...")
257 export_to_phy(
258 analyzer,
259 output_folder=output_path / 'phy_export',
260 copy_binary=True,
261 additional_properties=['curation_label'],
262 )
264 print(f"\nAnalysis complete!")
265 print(f"Results saved to: {output_path}")
266 print(f"\nTo open in Phy:")
267 print(f" phy template-gui {output_path / 'phy_export' / 'params.py'}")
270if __name__ == '__main__':
271 main()