Skills › Scientific Agent Skills › neuropixels-analysis › explore_recording.py
explore_recording.py
175 lines · skills/neuropixels-analysis/scripts/explore_recording.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Quick exploration of Neuropixels recording.
5Usage:
6 python explore_recording.py /path/to/spikeglx/data
7"""
9import argparse
10import spikeinterface.full as si
11import matplotlib.pyplot as plt
12import numpy as np
13from _common import validate_recording
16def explore_recording(data_path: str, stream_name: str = 'imec0.ap'):
17 """Explore a Neuropixels recording."""
19 print(f"Loading: {data_path}")
20 recording = si.read_spikeglx(data_path, stream_name=stream_name)
22 validate_recording(recording)
24 # Basic info
25 print("\n" + "="*50)
26 print("RECORDING INFO")
27 print("="*50)
28 print(f"Channels: {recording.get_num_channels()}")
29 print(f"Duration: {recording.get_total_duration():.2f} s ({recording.get_total_duration()/60:.2f} min)")
30 print(f"Sampling rate: {recording.get_sampling_frequency()} Hz")
31 print(f"Total samples: {recording.get_num_samples()}")
33 # Probe info
34 probe = recording.get_probe()
35 print(f"\nProbe: {probe.manufacturer} {probe.model_name if hasattr(probe, 'model_name') else ''}")
36 print(f"Probe shape: {probe.ndim}D")
38 # Channel groups
39 if recording.get_channel_groups() is not None:
40 groups = np.unique(recording.get_channel_groups())
41 print(f"Channel groups (shanks): {len(groups)}")
43 # Check for bad channels
44 print("\n" + "="*50)
45 print("BAD CHANNEL DETECTION")
46 print("="*50)
47 bad_ids, labels = si.detect_bad_channels(recording)
48 if len(bad_ids) > 0:
49 print(f"Bad channels found: {len(bad_ids)}")
50 for ch, label in zip(recording.channel_ids, labels):
51 if ch not in bad_ids:
52 continue
53 print(f" Channel {ch}: {label}")
54 else:
55 print("No bad channels detected")
57 # Sample traces
58 print("\n" + "="*50)
59 print("SIGNAL STATISTICS")
60 print("="*50)
62 # Get 1 second of data
63 n_samples = min(int(recording.get_sampling_frequency()), recording.get_num_samples())
64 traces = recording.get_traces(start_frame=0, end_frame=n_samples)
66 print(f"Sample mean: {np.mean(traces):.2f}")
67 print(f"Sample std: {np.std(traces):.2f}")
68 print(f"Sample min: {np.min(traces):.2f}")
69 print(f"Sample max: {np.max(traces):.2f}")
71 return recording
74def plot_probe(recording, output_path=None):
75 """Plot probe layout."""
76 fig, ax = plt.subplots(figsize=(4, 12))
77 si.plot_probe_map(recording, ax=ax, with_channel_ids=False)
78 ax.set_title('Probe Layout')
80 if output_path:
81 plt.savefig(output_path, dpi=150, bbox_inches='tight')
82 print(f"Saved: {output_path}")
83 else:
84 plt.show()
87def plot_traces(recording, duration=1.0, output_path=None):
88 """Plot raw traces."""
89 n_samples = min(int(duration * recording.get_sampling_frequency()), recording.get_num_samples())
90 if n_samples < 1:
91 raise ValueError("Plot duration must contain at least one sample.")
92 traces = recording.get_traces(start_frame=0, end_frame=n_samples)
94 fig, ax = plt.subplots(figsize=(12, 8))
96 # Plot subset of channels
97 n_channels = min(20, recording.get_num_channels())
98 channel_idx = np.linspace(0, recording.get_num_channels()-1, n_channels, dtype=int)
100 time = np.arange(n_samples) / recording.get_sampling_frequency()
102 for i, ch in enumerate(channel_idx):
103 offset = i * 200 # Offset for visibility
104 ax.plot(time, traces[:, ch] + offset, 'k', linewidth=0.5)
106 ax.set_xlabel('Time (s)')
107 ax.set_ylabel('Stored trace units (offset)')
108 ax.set_title(f'Raw Traces ({n_channels} channels)')
110 if output_path:
111 plt.savefig(output_path, dpi=150, bbox_inches='tight')
112 print(f"Saved: {output_path}")
113 else:
114 plt.show()
117def plot_power_spectrum(recording, output_path=None):
118 """Plot power spectrum."""
119 from scipy import signal
121 # Get data from middle channel
122 mid_ch = recording.get_num_channels() // 2
123 n_samples = min(int(10 * recording.get_sampling_frequency()), recording.get_num_samples())
125 traces = recording.get_traces(
126 start_frame=0,
127 end_frame=n_samples,
128 channel_ids=[recording.channel_ids[mid_ch]]
129 ).flatten()
131 fs = recording.get_sampling_frequency()
133 # Compute power spectrum
134 freqs, psd = signal.welch(traces, fs, nperseg=4096)
136 fig, ax = plt.subplots(figsize=(10, 5))
137 ax.semilogy(freqs, psd)
138 ax.set_xlabel('Frequency (Hz)')
139 ax.set_ylabel('Power Spectral Density')
140 ax.set_title(f'Power Spectrum (Channel {mid_ch})')
141 ax.set_xlim(0, 5000)
142 ax.axvline(300, color='r', linestyle='--', alpha=0.5, label='300 Hz')
143 ax.axvline(6000, color='r', linestyle='--', alpha=0.5, label='6000 Hz')
144 ax.legend()
145 ax.grid(True, alpha=0.3)
147 if output_path:
148 plt.savefig(output_path, dpi=150, bbox_inches='tight')
149 print(f"Saved: {output_path}")
150 else:
151 plt.show()
154if __name__ == '__main__':
155 parser = argparse.ArgumentParser(description='Explore Neuropixels recording')
156 parser.add_argument('data_path', help='Path to SpikeGLX recording')
157 parser.add_argument('--stream', default='imec0.ap', help='Stream name (e.g. imec0.ap)')
158 parser.add_argument('--plot', action='store_true', help='Generate plots')
159 parser.add_argument('--output', default=None, help='Output directory for plots')
161 args = parser.parse_args()
163 recording = explore_recording(args.data_path, args.stream)
165 if args.plot:
166 import os
167 if args.output:
168 os.makedirs(args.output, exist_ok=True)
169 plot_probe(recording, f"{args.output}/probe_map.png")
170 plot_traces(recording, output_path=f"{args.output}/raw_traces.png")
171 plot_power_spectrum(recording, f"{args.output}/power_spectrum.png")
172 else:
173 plot_probe(recording)
174 plot_traces(recording)
175 plot_power_spectrum(recording)