SpikeInterface

# Basic installation
pixi add spikeinterface

# With specific sorters
pixi add kilosort  # For Kilosort3/4
pixi add mountainsort5

# Or with pip
pip install spikeinterface[full]

import spikeinterface as si
import spikeinterface.extractors as se

# Load from file
recording = se.read_intan('data.rhd')  # Intan
# or
recording = se.read_spikeglx('data_folder')  # Neuropixels

# Check properties
print(f"Duration: {recording.get_total_duration()}s")
print(f"Sampling rate: {recording.get_sampling_frequency()}Hz")
print(f"Channels: {recording.get_num_channels()}")

import spikeinterface.preprocessing as spre

# Bandpass filter
recording_f = spre.bandpass_filter(recording, freq_min=300, freq_max=6000)

# Common average reference
recording_cmr = spre.common_reference(recording_f, reference='global')

# Remove artifacts
recording_clean = spre.remove_artifacts(
    recording_cmr,
    triggers=artifact_times,
    ms_before=0.5,
    ms_after=2.0
)

import spikeinterface.sorters as ss

# Run Kilosort3
sorting = ss.run_sorter(
    'kilosort3',
    recording_clean,
    output_folder='sorted_output',
    verbose=True,
    docker_image=True  # Use Docker for isolation
)

# Check results
print(f"Found {len(sorting.get_unit_ids())} units")

import spikeinterface.postprocessing as spost

# Extract waveforms
waveforms = spost.extract_waveforms(
    recording_clean,
    sorting,
    folder='waveforms',
    max_spikes_per_unit=500
)

# Compute templates
templates = waveforms.get_all_templates()

# Compute quality metrics
import spikeinterface.qualitymetrics as sqm

metrics = sqm.compute_quality_metrics(
    waveforms,
    metric_names=['snr', 'isi_violation', 'firing_rate']
)
print(metrics)

import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.preprocessing as spre
import spikeinterface.sorters as ss
import spikeinterface.postprocessing as spost
import spikeinterface.qualitymetrics as sqm
import spikeinterface.curation as sc

# 1. Load data
recording = se.read_spikeglx('neuropixels_data/')

# 2. Preprocessing pipeline
recording = spre.phase_shift(recording)  # Correct for ADC shifts
recording = spre.bandpass_filter(recording, freq_min=300, freq_max=6000)
recording = spre.common_reference(recording, reference='global')

# 3. Spike sorting
sorting = ss.run_sorter(
    'kilosort3',
    recording,
    output_folder='sorting_output',
    detect_threshold=6,
    projection_threshold=[10, 4],
    preclust_threshold=8
)

# 4. Extract waveforms
waveforms = spost.extract_waveforms(
    recording,
    sorting,
    folder='waveforms',
    max_spikes_per_unit=500,
    ms_before=1.0,
    ms_after=2.0
)

# 5. Compute metrics
metrics = sqm.compute_quality_metrics(
    waveforms,
    metric_names=[
        'snr',
        'isi_violation',
        'firing_rate',
        'presence_ratio',
        'amplitude_cutoff'
    ]
)

# 6. Automatic curation
sorting_curated = sc.auto_curation(
    sorting,
    waveforms,
    min_snr=5,
    max_isi_violation=0.5,
    min_firing_rate=0.1
)

print(f"Original units: {len(sorting.get_unit_ids())}")
print(f"Curated units: {len(sorting_curated.get_unit_ids())}")

# 7. Export
sorting_curated.save(folder='final_sorting')

import spikeinterface.comparison as sc

# Run multiple sorters
sorting_ks3 = ss.run_sorter('kilosort3', recording, output_folder='ks3')
sorting_ms5 = ss.run_sorter('mountainsort5', recording, output_folder='ms5')

# Compare results
comparison = sc.compare_two_sorters(
    sorting_ks3,
    sorting_ms5,
    sorting1_name='Kilosort3',
    sorting2_name='MountainSort5'
)

# Get agreement matrix
agreement = comparison.get_agreement_fraction()
print(f"Agreement: {agreement:.2%}")

# Find consensus units
matches = comparison.get_matching()
print(f"Matched units: {len(matches)}")

# Visualize
import matplotlib.pyplot as plt
sc.plot_agreement_matrix(comparison)
plt.show()

metrics = sqm.compute_quality_metrics(
    waveforms,
    metric_names=[
        'snr',                    # Signal-to-noise ratio
        'isi_violation',          # ISI violations (%)
        'firing_rate',            # Mean firing rate
        'presence_ratio',         # Fraction of time active
        'amplitude_cutoff',       # Est. false negative rate
        'nn_isolation',           # Nearest-neighbor isolation
        'nn_noise_overlap'        # Overlap with noise
    ]
)

import pandas as pd

# Set thresholds
high_quality = metrics[
    (metrics['snr'] > 5) &
    (metrics['isi_violation'] < 0.5) &
    (metrics['firing_rate'] > 0.1) &
    (metrics['presence_ratio'] > 0.9)
]

print(f"High-quality units: {len(high_quality)}")

import spikeinterface.widgets as sw

# Plot unit waveforms
sw.plot_unit_waveforms(waveforms, unit_ids=[0, 1, 2])

# Plot unit templates
sw.plot_unit_templates(waveforms, unit_ids=[0, 1, 2])

# Plot amplitudes over time
sw.plot_amplitudes(waveforms, unit_ids=[0, 1, 2])

# Plot ISI distribution
sw.plot_isi_distribution(sorting, unit_ids=[0, 1, 2])

# Plot autocorrelograms
sw.plot_autocorrelograms(sorting, unit_ids=[0, 1, 2])

# Plot crosscorrelograms
sw.plot_crosscorrelograms(sorting, unit_ids=[0, 1])

# Export for manual curation in Phy
spost.export_to_phy(
    waveforms,
    output_folder='phy_output',
    compute_pc_features=True,
    compute_amplitudes=True
)

# After manual curation in Phy, load back:
sorting_phy = se.read_phy('phy_output')

# Get spike times for a unit
spike_times = sorting.get_unit_spike_train(unit_id=0)
print(f"Unit 0: {len(spike_times)} spikes")

# Get all spike trains
for unit_id in sorting.get_unit_ids():
    spike_train = sorting.get_unit_spike_train(unit_id)
    rate = len(spike_train) / recording.get_total_duration()
    print(f"Unit {unit_id}: {rate:.2f} Hz")

# Convert to Neo format
import neo

spike_train_neo = neo.SpikeTrain(
    spike_times / recording.get_sampling_frequency(),
    t_stop=recording.get_total_duration(),
    units='s'
)

# Highpass filter
recording_hp = spre.highpass_filter(recording, freq_min=300)

# Bandpass filter
recording_bp = spre.bandpass_filter(recording, freq_min=300, freq_max=6000)

# Notch filter (remove 50/60 Hz)
recording_notch = spre.notch_filter(recording, freq=60, q=30)

# Common average reference
recording_car = spre.common_reference(recording, reference='global')

# Common median reference
recording_cmr = spre.common_reference(recording, reference='median')

# Specific reference channels
recording_ref = spre.common_reference(recording, ref_channel_ids=[0, 1, 2])

# Remove time periods
recording_clean = spre.remove_artifacts(
    recording,
    list_triggers=artifact_frames,
    ms_before=1.0,
    ms_after=3.0
)

# Interpolate artifacts
recording_interp = spre.remove_artifacts(
    recording,
    list_triggers=artifact_frames,
    mode='linear'  # or 'cubic'
)