Skill Tracer
Skills › Scientific Agent Skills › neuropixels-analysis › ANALYSIS.md

ANALYSIS.md

398 lines · skills/neuropixels-analysis/references/ANALYSIS.md at 5f2bc03

Open on GitHub
1# Post-Processing & Analysis Reference
2
3Targets SI 0.105.0 (reviewed 2026-10-01). Recording-dependent examples are illustrative.
4Use [QUALITY_METRICS.md](QUALITY_METRICS.md) for metric dependencies and caveats;
5result columns below are not necessarily computation names.
6

Post-Processing & Analysis Reference

Targets SI 0.105.0 (reviewed 2026-10-01). Recording-dependent examples are illustrative. Use QUALITY_METRICS.md for metric dependencies and caveats; result columns below are not necessarily computation names.

Raw lines 1 to 6
7## Sorting Analyzer
8
9The `SortingAnalyzer` is the central object for post-processing.

Sorting Analyzer

The SortingAnalyzer is the central object for post-processing.

Raw lines 7 to 10
11### Create Analyzer
12```python
13import spikeinterface.full as si
15# Create analyzer
16analyzer = si.create_sorting_analyzer(
17 sorting,
18 recording,
19 sparse=True, # Use sparse representation
20 format='binary_folder', # Storage format
21 folder='analyzer_output' # Save location
23```

Create Analyzer

import spikeinterface.full as si

# Create analyzer
analyzer = si.create_sorting_analyzer(
    sorting,
    recording,
    sparse=True,                    # Use sparse representation
    format='binary_folder',         # Storage format
    folder='analyzer_output'        # Save location
)
Raw lines 11 to 24
25### Compute Extensions
26```python
27# Compute all standard extensions
28analyzer.compute('random_spikes') # Random spike selection
29analyzer.compute('waveforms') # Extract waveforms
30analyzer.compute('templates') # Compute templates
31analyzer.compute('noise_levels') # Noise estimation
32analyzer.compute('principal_components') # PCA
33analyzer.compute('spike_amplitudes')
34analyzer.compute('amplitude_scalings') # Amplitude per spike
35analyzer.compute('correlograms') # Auto/cross correlograms
36analyzer.compute('unit_locations') # Unit locations
37analyzer.compute('spike_locations') # Per-spike locations
38analyzer.compute('template_similarity') # Template similarity matrix
39analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate']) # Quality metrics
41# Or compute multiple at once
42analyzer.compute([
43 'random_spikes', 'waveforms', 'templates', 'noise_levels',
44 'principal_components', 'spike_amplitudes', 'amplitude_scalings', 'correlograms',
45 'unit_locations', 'quality_metrics'
46])
47```

Compute Extensions

# Compute all standard extensions
analyzer.compute('random_spikes')       # Random spike selection
analyzer.compute('waveforms')           # Extract waveforms
analyzer.compute('templates')           # Compute templates
analyzer.compute('noise_levels')        # Noise estimation
analyzer.compute('principal_components')  # PCA
analyzer.compute('spike_amplitudes')
analyzer.compute('amplitude_scalings')    # Amplitude per spike
analyzer.compute('correlograms')        # Auto/cross correlograms
analyzer.compute('unit_locations')      # Unit locations
analyzer.compute('spike_locations')     # Per-spike locations
analyzer.compute('template_similarity') # Template similarity matrix
analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate'])     # Quality metrics

# Or compute multiple at once
analyzer.compute([
    'random_spikes', 'waveforms', 'templates', 'noise_levels',
    'principal_components', 'spike_amplitudes', 'amplitude_scalings', 'correlograms',
    'unit_locations', 'quality_metrics'
])
Raw lines 25 to 48
49### Save and Load
50```python
51# Save
52analyzer.save_as(folder='analyzer_saved', format='binary_folder')
54# Load
55analyzer = si.load_sorting_analyzer('analyzer_saved')
56```

Save and Load

# Save
analyzer.save_as(folder='analyzer_saved', format='binary_folder')

# Load
analyzer = si.load_sorting_analyzer('analyzer_saved')
Raw lines 49 to 57
58## Quality Metrics

Quality Metrics

Raw lines 58 to 59
60### Compute Metrics
61```python
62analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate'])
63qm = analyzer.get_extension('quality_metrics').get_data()
64print(qm)
65```

Compute Metrics

analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate'])
qm = analyzer.get_extension('quality_metrics').get_data()
print(qm)
Raw lines 60 to 66
67### Available Metrics
69| Metric | Description | Good Values |
70|--------|-------------|-------------|
71| `snr` | Signal-to-noise ratio | > 5 |
72| `isi_violations_ratio` | Rate-normalized ISI ratio | < 0.01 (not 1%) |
73| `isi_violations_count` | ISI violation count | Low |
74| `presence_ratio` | Fraction of recording with spikes | > 0.9 |
75| `firing_rate` | Spikes per second | 0.1-50 Hz |
76| `amplitude_cutoff` | Estimated missed spikes | < 0.1 |
77| `amplitude_median` | Median spike amplitude | - |
78| `amplitude_cv` | Coefficient of variation | < 0.5 |
79| `drift_ptp` | Peak-to-peak drift (um) | < 40 |
80| `drift_std` | Standard deviation of drift | < 10 |
81| `drift_mad` | Median absolute deviation | < 10 |
82| `sliding_rp_violation` | Sliding refractory period | < 0.05 |
83| `sync_spike_2` | Synchrony with other units | < 0.5 |
84| `isolation_distance` | Mahalanobis distance | > 20 |
85| `l_ratio` | L-ratio (isolation) | < 0.1 |
86| `d_prime` | Discriminability | > 5 |
87| `nn_hit_rate` | Nearest neighbor hit rate | > 0.9 |
88| `nn_miss_rate` | Nearest neighbor miss rate | < 0.1 |
89| `silhouette_score` | Cluster silhouette | > 0.5 |

Available Metrics

MetricDescriptionGood Values
snrSignal-to-noise ratio> 5
isi_violations_ratioRate-normalized ISI ratio< 0.01 (not 1%)
isi_violations_countISI violation countLow
presence_ratioFraction of recording with spikes> 0.9
firing_rateSpikes per second0.1-50 Hz
amplitude_cutoffEstimated missed spikes< 0.1
amplitude_medianMedian spike amplitude-
amplitude_cvCoefficient of variation< 0.5
drift_ptpPeak-to-peak drift (um)< 40
drift_stdStandard deviation of drift< 10
drift_madMedian absolute deviation< 10
sliding_rp_violationSliding refractory period< 0.05
sync_spike_2Synchrony with other units< 0.5
isolation_distanceMahalanobis distance> 20
l_ratioL-ratio (isolation)< 0.1
d_primeDiscriminability> 5
nn_hit_rateNearest neighbor hit rate> 0.9
nn_miss_rateNearest neighbor miss rate< 0.1
silhouette_scoreCluster silhouette> 0.5
Raw lines 67 to 90
91### Compute Specific Metrics
92```python
93analyzer.compute(
94 'quality_metrics',
95 metric_names=['snr', 'isi_violation', 'presence_ratio', 'firing_rate']
97```

Compute Specific Metrics

analyzer.compute(
    'quality_metrics',
    metric_names=['snr', 'isi_violation', 'presence_ratio', 'firing_rate']
)
Raw lines 91 to 98
99### Custom Quality Thresholds
100```python
101qm = analyzer.get_extension('quality_metrics').get_data()
103# Define quality criteria
104quality_criteria = {
105 'snr': ('>', 5),
106 'isi_violations_ratio': ('<', 0.01),
107 'presence_ratio': ('>', 0.9),
108 'firing_rate': ('>', 0.1),
109 'amplitude_cutoff': ('<', 0.1),
112# Filter good units
113good_units = qm.query(
114 "(snr > 5) & (isi_violations_ratio < 0.01) & (presence_ratio > 0.9)"
115).index.tolist()
117print(f"Good units: {len(good_units)}/{len(qm)}")
118```

Custom Quality Thresholds

qm = analyzer.get_extension('quality_metrics').get_data()

# Define quality criteria
quality_criteria = {
    'snr': ('>', 5),
    'isi_violations_ratio': ('<', 0.01),
    'presence_ratio': ('>', 0.9),
    'firing_rate': ('>', 0.1),
    'amplitude_cutoff': ('<', 0.1),
}

# Filter good units
good_units = qm.query(
    "(snr > 5) & (isi_violations_ratio < 0.01) & (presence_ratio > 0.9)"
).index.tolist()

print(f"Good units: {len(good_units)}/{len(qm)}")
Raw lines 99 to 119
120## Waveforms & Templates

Waveforms & Templates

Raw lines 120 to 121
122### Extract Waveforms
123```python
124analyzer.compute('random_spikes', max_spikes_per_unit=500, seed=42)
125analyzer.compute('waveforms', ms_before=1.5, ms_after=2.5)
127# Get waveforms for a unit
128waveforms = analyzer.get_extension('waveforms').get_waveforms_one_unit(unit_id=0)
129print(f"Shape: {waveforms.shape}") # (n_spikes, n_samples, n_channels)
130```

Extract Waveforms

analyzer.compute('random_spikes', max_spikes_per_unit=500, seed=42)
analyzer.compute('waveforms', ms_before=1.5, ms_after=2.5)

# Get waveforms for a unit
waveforms = analyzer.get_extension('waveforms').get_waveforms_one_unit(unit_id=0)
print(f"Shape: {waveforms.shape}")  # (n_spikes, n_samples, n_channels)
Raw lines 122 to 131
132### Compute Templates
133```python
134analyzer.compute('templates', operators=['average', 'std', 'median'])
136# Get template
137templates_ext = analyzer.get_extension('templates')
138template = templates_ext.get_unit_template(unit_id=0, operator='average')
139```

Compute Templates

analyzer.compute('templates', operators=['average', 'std', 'median'])

# Get template
templates_ext = analyzer.get_extension('templates')
template = templates_ext.get_unit_template(unit_id=0, operator='average')
Raw lines 132 to 140
141### Template Similarity
142```python
143analyzer.compute('template_similarity')
144sim = analyzer.get_extension('template_similarity').get_data()
145# Matrix of cosine similarities between templates
146```

Template Similarity

analyzer.compute('template_similarity')
sim = analyzer.get_extension('template_similarity').get_data()
# Matrix of cosine similarities between templates
Raw lines 141 to 147
148## Unit Locations

Unit Locations

Raw lines 148 to 149
150### Compute Locations
151```python
152analyzer.compute('unit_locations', method='monopolar_triangulation')
153locations = analyzer.get_extension('unit_locations').get_data()
154print(locations) # x, y coordinates per unit
155```

Compute Locations

analyzer.compute('unit_locations', method='monopolar_triangulation')
locations = analyzer.get_extension('unit_locations').get_data()
print(locations)  # x, y coordinates per unit
Raw lines 150 to 156
157### Spike Locations
158```python
159analyzer.compute('spike_locations', method='center_of_mass')
160spike_locs = analyzer.get_extension('spike_locations').get_data()
161```

Spike Locations

analyzer.compute('spike_locations', method='center_of_mass')
spike_locs = analyzer.get_extension('spike_locations').get_data()
Raw lines 157 to 162
163### Location Methods
164- `'center_of_mass'` - Fast, less accurate
165- `'monopolar_triangulation'` - More accurate, slower
166- `'grid_convolution'` - Good balance

Location Methods

  • 'center_of_mass' - Fast, less accurate
  • 'monopolar_triangulation' - More accurate, slower
  • 'grid_convolution' - Good balance
Raw lines 163 to 167
168## Correlograms

Correlograms

Raw lines 168 to 169
170### Auto-correlograms
171```python
172analyzer.compute('correlograms', window_ms=50, bin_ms=1)
173correlograms, bins = analyzer.get_extension('correlograms').get_data()
175# correlograms shape: (n_units, n_units, n_bins)
176# Auto-correlogram for unit i: correlograms[i, i, :]
177# Cross-correlogram units i,j: correlograms[i, j, :]
178```

Auto-correlograms

analyzer.compute('correlograms', window_ms=50, bin_ms=1)
correlograms, bins = analyzer.get_extension('correlograms').get_data()

# correlograms shape: (n_units, n_units, n_bins)
# Auto-correlogram for unit i: correlograms[i, i, :]
# Cross-correlogram units i,j: correlograms[i, j, :]
Raw lines 170 to 179
180## Visualization

Visualization

Raw lines 180 to 181
182### Probe Map
183```python
184si.plot_probe_map(recording, with_channel_ids=True)
185```

Probe Map

si.plot_probe_map(recording, with_channel_ids=True)
Raw lines 182 to 186
187### Unit Templates
188```python
189# All units
190si.plot_unit_templates(analyzer)
192# Specific units
193si.plot_unit_templates(analyzer, unit_ids=[0, 1, 2])
194```

Unit Templates

# All units
si.plot_unit_templates(analyzer)

# Specific units
si.plot_unit_templates(analyzer, unit_ids=[0, 1, 2])
Raw lines 187 to 195
196### Waveforms
197```python
198# Plot waveforms with template
199si.plot_unit_waveforms(analyzer, unit_ids=[0])
201# Waveform density
202si.plot_unit_waveforms_density_map(analyzer, unit_ids=[0])
203```

Waveforms

# Plot waveforms with template
si.plot_unit_waveforms(analyzer, unit_ids=[0])

# Waveform density
si.plot_unit_waveforms_density_map(analyzer, unit_ids=[0])
Raw lines 196 to 204
205### Raster Plot
206```python
207si.plot_rasters(sorting, time_range=(0, 10)) # First 10 seconds
208```

Raster Plot

si.plot_rasters(sorting, time_range=(0, 10))  # First 10 seconds
Raw lines 205 to 209
210### Amplitudes
211```python
212analyzer.compute('spike_amplitudes')
213analyzer.compute('amplitude_scalings')
214si.plot_amplitudes(analyzer)
216# Distribution
217si.plot_all_amplitudes_distributions(analyzer)
218```

Amplitudes

analyzer.compute('spike_amplitudes')
analyzer.compute('amplitude_scalings')
si.plot_amplitudes(analyzer)

# Distribution
si.plot_all_amplitudes_distributions(analyzer)
Raw lines 210 to 219
220### Correlograms
221```python
222# Auto-correlograms
223si.plot_autocorrelograms(analyzer, unit_ids=[0, 1, 2])
225# Cross-correlograms
226si.plot_crosscorrelograms(analyzer, unit_ids=[0, 1])
227```

Correlograms

# Auto-correlograms
si.plot_autocorrelograms(analyzer, unit_ids=[0, 1, 2])

# Cross-correlograms
si.plot_crosscorrelograms(analyzer, unit_ids=[0, 1])
Raw lines 220 to 228
229### Quality Metrics
230```python
231# Summary plot
232si.plot_quality_metrics(analyzer)
234# Specific metric distribution
235import matplotlib.pyplot as plt
236qm = analyzer.get_extension('quality_metrics').get_data()
237plt.hist(qm['snr'], bins=50)
238plt.xlabel('SNR')
239plt.ylabel('Count')
240```

Quality Metrics

# Summary plot
si.plot_quality_metrics(analyzer)

# Specific metric distribution
import matplotlib.pyplot as plt
qm = analyzer.get_extension('quality_metrics').get_data()
plt.hist(qm['snr'], bins=50)
plt.xlabel('SNR')
plt.ylabel('Count')
Raw lines 229 to 241
242### Unit Locations on Probe
243```python
244si.plot_unit_locations(analyzer)
245```

Unit Locations on Probe

si.plot_unit_locations(analyzer)
Raw lines 242 to 246
247### Drift Map
248```python
249analyzer.compute("spike_locations", method="center_of_mass")
250si.plot_drift_raster_map(sorting_analyzer=analyzer)
251```

Drift Map

analyzer.compute("spike_locations", method="center_of_mass")
si.plot_drift_raster_map(sorting_analyzer=analyzer)
Raw lines 247 to 252
253### Summary Plot
254```python
255# Comprehensive unit summary
256si.plot_unit_summary(analyzer, unit_id=0)
257```

Summary Plot

# Comprehensive unit summary
si.plot_unit_summary(analyzer, unit_id=0)
Raw lines 253 to 258
259## LFP Analysis

LFP Analysis

Raw lines 259 to 260
261### Load LFP Data
262```python
263lfp = si.read_spikeglx('/path/to/data', stream_name='imec0.lf')
264print(f"LFP: {lfp.get_sampling_frequency()} Hz")
265```

Load LFP Data

lfp = si.read_spikeglx('/path/to/data', stream_name='imec0.lf')
print(f"LFP: {lfp.get_sampling_frequency()} Hz")
Raw lines 261 to 266
267### Basic LFP Processing
268```python
269# Downsample if needed
270lfp_ds = si.resample(lfp, resample_rate=1000)
272# Common average reference
273lfp_car = si.common_reference(lfp_ds, reference='global', operator='median')
274```

Basic LFP Processing

# Downsample if needed
lfp_ds = si.resample(lfp, resample_rate=1000)

# Common average reference
lfp_car = si.common_reference(lfp_ds, reference='global', operator='median')
Raw lines 267 to 275
276### Extract LFP Traces
277```python
278import numpy as np
280# Get traces (samples x channels), in calibrated microvolts
281traces = lfp.get_traces(start_frame=0, end_frame=min(30000, lfp.get_num_samples()), return_in_uV=True)
283# Specific channels
284traces = lfp.get_traces(channel_ids=[0, 1, 2])
285```

Extract LFP Traces

import numpy as np

# Get traces (samples x channels), in calibrated microvolts
traces = lfp.get_traces(start_frame=0, end_frame=min(30000, lfp.get_num_samples()), return_in_uV=True)

# Specific channels
traces = lfp.get_traces(channel_ids=[0, 1, 2])
Raw lines 276 to 286
287### Spectral Analysis
288```python
289from scipy import signal
290import matplotlib.pyplot as plt
292# Get single channel
293trace = lfp.get_traces(channel_ids=[lfp.channel_ids[0]],
294 end_frame=min(lfp.get_num_samples(), int(60 * lfp.get_sampling_frequency())),
295 return_in_uV=True).flatten()
296fs = lfp.get_sampling_frequency()
298# Power spectrum
299freqs, psd = signal.welch(trace, fs, nperseg=4096)
300plt.semilogy(freqs, psd)
301plt.xlabel('Frequency (Hz)')
302plt.ylabel('PSD (uV squared / Hz)')
303plt.xlim(0, 100)
304```

Spectral Analysis

from scipy import signal
import matplotlib.pyplot as plt

# Get single channel
trace = lfp.get_traces(channel_ids=[lfp.channel_ids[0]],
                       end_frame=min(lfp.get_num_samples(), int(60 * lfp.get_sampling_frequency())),
                       return_in_uV=True).flatten()
fs = lfp.get_sampling_frequency()

# Power spectrum
freqs, psd = signal.welch(trace, fs, nperseg=4096)
plt.semilogy(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (uV squared / Hz)')
plt.xlim(0, 100)
Raw lines 287 to 305
306### Spectrogram
307```python
308f, t, Sxx = signal.spectrogram(trace, fs, nperseg=2048, noverlap=1024)
309plt.pcolormesh(t, f, 10*np.log10(Sxx), shading='gouraud')
310plt.ylabel('Frequency (Hz)')
311plt.xlabel('Time (s)')
312plt.ylim(0, 100)
313plt.colorbar(label='Power (dB)')
314```

Spectrogram

f, t, Sxx = signal.spectrogram(trace, fs, nperseg=2048, noverlap=1024)
plt.pcolormesh(t, f, 10*np.log10(Sxx), shading='gouraud')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.ylim(0, 100)
plt.colorbar(label='Power (dB)')
Raw lines 306 to 315
316## Export Formats

Export Formats

Raw lines 316 to 317
318### Export to Phy
319```python
320si.export_to_phy(
321 analyzer,
322 output_folder='phy_export',
323 compute_pc_features=True,
324 compute_amplitudes=True,
325 copy_binary=True
327# Then: phy template-gui phy_export/params.py
328```

Export to Phy

si.export_to_phy(
    analyzer,
    output_folder='phy_export',
    compute_pc_features=True,
    compute_amplitudes=True,
    copy_binary=True
)
# Then: phy template-gui phy_export/params.py
Raw lines 318 to 329
330### Export to NWB
332There is no `spikeinterface.exporters.export_to_nwb` in SI 0.105.0.
333Use [NeuroConv's NWBConverter](https://neuroconv.readthedocs.io/en/main/user_guide/nwbconverter.html)
334with the source recording and sorting interfaces, then validate the NWB output.
335Supply session identity, timezone-aware start time, electrode/probe metadata,
336calibration, and a shared clock for ElectricalSeries and units. A units-only
337conversion does not preserve raw voltage/electrode metadata automatically. This
338optional workflow is documentation-verified, not executed by the bundled scripts.

Export to NWB

There is no spikeinterface.exporters.export_to_nwb in SI 0.105.0. Use NeuroConv's NWBConverter with the source recording and sorting interfaces, then validate the NWB output. Supply session identity, timezone-aware start time, electrode/probe metadata, calibration, and a shared clock for ElectricalSeries and units. A units-only conversion does not preserve raw voltage/electrode metadata automatically. This optional workflow is documentation-verified, not executed by the bundled scripts.

Raw lines 330 to 339
340### Export Report
341```python
342si.export_report(
343 analyzer,
344 output_folder='report',
345 remove_if_exists=True,
346 format='png'
348```

Export Report

si.export_report(
    analyzer,
    output_folder='report',
    remove_if_exists=True,
    format='png'
)
Raw lines 340 to 349
350## Complete Analysis Pipeline
352```python
353import spikeinterface.full as si
355def analyze_sorting(recording, sorting, output_dir):
356 """Complete post-processing pipeline."""
358 # Create analyzer
359 analyzer = si.create_sorting_analyzer(
360 sorting, recording,
361 sparse=True, format='binary_folder',
362 folder=f'{output_dir}/analyzer'
363 )
365 # Compute all extensions
366 print("Computing extensions...")
367 analyzer.compute(['random_spikes', 'waveforms', 'templates', 'noise_levels'])
368 analyzer.compute(['principal_components', 'spike_amplitudes', 'amplitude_scalings'])
369 analyzer.compute(['correlograms', 'unit_locations', 'template_similarity'])
370 analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate'])
372 # Get quality metrics
373 qm = analyzer.get_extension('quality_metrics').get_data()
375 # Filter good units
376 good_units = qm.query(
377 "(snr > 5) & (isi_violations_ratio < 0.01) & (presence_ratio > 0.9)"
378 ).index.tolist()
380 print(f"Quality filtering: {len(good_units)}/{len(qm)} units passed")
382 # Export
383 si.export_to_phy(analyzer, f'{output_dir}/phy')
384 si.export_report(analyzer, f'{output_dir}/report')
386 # Save metrics
387 qm.to_csv(f'{output_dir}/quality_metrics.csv')
389 return analyzer, qm, good_units
391# Usage
392analyzer, qm, good_units = analyze_sorting(recording, sorting, 'output/')
393```
395SpikeInterface 0.105.0 has an observed `read_phy` bug for exported nonnumeric
396unit IDs (`np.isnan` TypeError). Phy export preserves `cluster_si_unit_ids.tsv`;
397keep that mapping and use a validated importer/fixed release for string-ID
398readback. Numeric-ID Phy export/reload was tested on synthetic data.

Complete Analysis Pipeline

import spikeinterface.full as si

def analyze_sorting(recording, sorting, output_dir):
    """Complete post-processing pipeline."""

    # Create analyzer
    analyzer = si.create_sorting_analyzer(
        sorting, recording,
        sparse=True, format='binary_folder',
        folder=f'{output_dir}/analyzer'
    )

    # Compute all extensions
    print("Computing extensions...")
    analyzer.compute(['random_spikes', 'waveforms', 'templates', 'noise_levels'])
    analyzer.compute(['principal_components', 'spike_amplitudes', 'amplitude_scalings'])
    analyzer.compute(['correlograms', 'unit_locations', 'template_similarity'])
    analyzer.compute('quality_metrics', metric_names=['snr', 'isi_violation', 'presence_ratio', 'amplitude_cutoff', 'firing_rate'])

    # Get quality metrics
    qm = analyzer.get_extension('quality_metrics').get_data()

    # Filter good units
    good_units = qm.query(
        "(snr > 5) & (isi_violations_ratio < 0.01) & (presence_ratio > 0.9)"
    ).index.tolist()

    print(f"Quality filtering: {len(good_units)}/{len(qm)} units passed")

    # Export
    si.export_to_phy(analyzer, f'{output_dir}/phy')
    si.export_report(analyzer, f'{output_dir}/report')

    # Save metrics
    qm.to_csv(f'{output_dir}/quality_metrics.csv')

    return analyzer, qm, good_units

# Usage
analyzer, qm, good_units = analyze_sorting(recording, sorting, 'output/')

SpikeInterface 0.105.0 has an observed read_phy bug for exported nonnumeric unit IDs (np.isnan TypeError). Phy export preserves cluster_si_unit_ids.tsv; keep that mapping and use a validated importer/fixed release for string-ID readback. Numeric-ID Phy export/reload was tested on synthetic data.

Raw lines 350 to 398