<jupyter_start><jupyter_text>SpikeInterface pipeline for Movshon Lab - Blackrock<jupyter_code>from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from pprint import pprint

import spikeextractors as se
import spiketoolkit as st
import spikesorters as ss
import spikecomparison as sc
import spikewidgets as sw<jupyter_output><empty_output><jupyter_text>1) Load recordings, compute LFP, and inspect signals<jupyter_code># Data files directory
dir_path = Path('/home/luiz/storage/taufferconsulting/client_ben/project_movshon/data_blackrock')

# Test if file exists
file_path = dir_path / 'XX_LE_textures_20191128_002.ns6'
print(f'File exists: {file_path.is_file()}')

# Spikeinterface directory
dir_spikeinterface = dir_path / "spikeinterface"
dir_spikeinterface.mkdir(parents=True, exist_ok=True)
print(dir_spikeinterface)
# Choose Blackrock extension (ns*) to be used
nsx_to_load = 6

# Create recording extractor
recording = se.BlackrockRecordingExtractor(
    filename=str(file_path),
    nsx_to_load=nsx_to_load
)

# Set individual groups
#recording.set_channel_groups(range(recording.get_num_channels()))<jupyter_output><empty_output><jupyter_text>Stub recording for fast testing; set to False for running processing pipeline on entire data<jupyter_code>print(f"Num channels: {recording.get_num_channels()}")
print(f"Sampling rate: {recording.get_sampling_frequency()}")
print(f"Duration (s): {recording.get_num_frames() / recording.get_sampling_frequency()}")<jupyter_output>Num channels: 192
Sampling rate: 30000.0
Duration (s): 1942.1709333333333<jupyter_text>Compute LFP<jupyter_code>freq_min_lfp = 1
freq_max_lfp = 300
freq_resample_lfp = 1000.

# Apply bandpass filter
recording_lfp = st.preprocessing.bandpass_filter(recording, freq_min=freq_min_lfp, freq_max=freq_max_lfp)

# Resample lfp
recording_lfp = st.preprocessing.resample(recording_lfp, freq_resample_lfp)

print(f"Sampling frequency Raw: {recording.get_sampling_frequency()}")
print(f"Sampling frequency LF: {recording_lfp.get_sampling_frequency()}")<jupyter_output>Sampling frequency Raw: 30000.0
Sampling frequency LF: 1000.0<jupyter_text>Inspect signals<jupyter_code>w_ts_raw = sw.plot_timeseries(recording, trange=[0, 5], channel_ids=[1, 2, 3])
# w_ts_lf = sw.plot_timeseries(recording_lfp, trange=[10, 15])<jupyter_output><empty_output><jupyter_text>2) Pre-processing<jupyter_code>apply_filter = True
apply_cmr = True
freq_min_hp = 300
freq_max_hp = 3000
if apply_filter:
    recording_processed = st.preprocessing.bandpass_filter(recording, freq_min=freq_min_hp, freq_max=freq_max_hp)
else:
    recording_processed = recording

if apply_cmr:
    recording_processed = st.preprocessing.common_reference(recording_processed)
    
# Stub recording for fast testing; set to False for running processing pipeline on entire data
stub_test = True
nsec_stub = 20
subr_ids = [i + 1 for i in range(10)]
if stub_test:
    recording_processed = se.SubRecordingExtractor(
        parent_recording=recording_processed, 
        channel_ids=subr_ids,
        end_frame=int(nsec_stub*recording_processed.get_sampling_frequency())
    )
    recording_lfp = se.SubRecordingExtractor(recording_lfp, end_frame=int(nsec_stub*recording_lfp.get_sampling_frequency()))
    
print(f"Original signal length: {recording.get_num_frames()}")
print(f"Processed signal length: {recording_processed.get_num_frames()}")
w_ts_processed = sw.plot_timeseries(recording_processed, trange=[0, 5], channel_ids=[1, 2, 3])<jupyter_output><empty_output><jupyter_text>3) Run spike sorters<jupyter_code>ss.installed_sorters()
sorter_list = ['klusta']

# Inspect sorter-specific parameters and defaults
for sorter in sorter_list:
    print(f"{sorter} params description:")
    pprint(ss.get_params_description(sorter))
    print("Default params:")
    pprint(ss.get_default_params(sorter))
# Choose which recording to use for sorting
rec_to_sort = recording_processed

# spike sort by group
sorting = ss.run_klusta(
    recording=rec_to_sort, 
    grouping_property='group',
    output_folder=dir_spikeinterface / "si_output",
)<jupyter_output>WARNING! No channel location given. Add dummy location.
RUNNING SHELL SCRIPT: /home/luiz/storage/taufferconsulting/client_ben/project_movshon/data_blackrock/spikeinterface/si_output/0/run_klusta.sh<jupyter_text>4) Post-processing: extract waveforms, templates, quality metrics, extracellular features  Set postprocessing parameters<jupyter_code># Post-processing params
postprocessing_params = st.postprocessing.get_common_params()
pprint(postprocessing_params)
# (optional) change parameters
postprocessing_params['max_spikes_per_unit'] = 1000  # with None, all waveforms are extracted<jupyter_output><empty_output><jupyter_text>Set quality metric list<jupyter_code># Quality metrics
qc_list = st.validation.get_quality_metrics_list()
print(f"Available quality metrics: {qc_list}")
# (optional) define subset of qc
qc_list = ["snr", "isi_violation", "firing_rate"]<jupyter_output><empty_output><jupyter_text>Set extracellular features<jupyter_code># Extracellular features
ec_list = st.postprocessing.get_template_features_list()
print(f"Available EC features: {ec_list}")
# (optional) define subset of ec
ec_list = ["peak_to_valley", "halfwidth"]<jupyter_output><empty_output><jupyter_text>Postprocess all sorting outputs<jupyter_code># st.validation.compute_quality_metrics?
tmp_folder = dir_spikeinterface / 'tmp' / 'klusta'
tmp_folder.mkdir(parents=True, exist_ok=True)

# set local tmp folder
sorting.set_tmp_folder(tmp_folder)

# compute waveforms
waveforms = st.postprocessing.get_unit_waveforms(rec_to_sort, sorting, **postprocessing_params)

# compute templates
templates = st.postprocessing.get_unit_templates(rec_to_sort, sorting, **postprocessing_params)

# comput EC features
ec = st.postprocessing.compute_unit_template_features(rec_to_sort, sorting,
                                                      feature_names=ec_list, as_dataframe=True)
# compute QCs
qc = st.validation.compute_quality_metrics(sorting, recording=rec_to_sort, 
                                           metric_names=qc_list, as_dataframe=True)

# export to phy
phy_folder = dir_spikeinterface / 'phy' / 'klusta'
phy_folder.mkdir(parents=True, exist_ok=True)
st.postprocessing.export_to_phy(rec_to_sort, sorting, phy_folder)
# Visualize spike template waverforms
unit_id = 0

plt.figure()
spk = np.squeeze(templates[unit_id])
plt.plot(spk.T)
plt.ylabel('Avg spike trace per channel', fontsize=12)<jupyter_output><empty_output><jupyter_text>5) Compare with baseline - Extracting sorted spikes from Blackrock .nev<jupyter_code>spikes_file = str(dir_path / 'XX_LE_textures_20191128_002.nev')

br_spike_extractor = se.BlackrockSortingExtractor(
    filename=spikes_file, 
    nsx_to_load=6
)
cmp_KL_BR = sc.compare_two_sorters(
    sorting1=sorting, sorting2=br_spike_extractor,
    sorting1_name='klusta', sorting2_name='blackrock'
)
# sw.plot_agreement_matrix(cmp_KL_BR)
cmp_KL_BR.agreement_scores<jupyter_output><empty_output><jupyter_text>6) Automatic curation<jupyter_code># define curators and thresholds
isi_violation_threshold = 0.5
snr_threshold = 5
firing_rate_threshold = 0.1
num_frames = rec_to_sort.get_num_frames()

# firing rate threshold
sorting_curated = st.curation.threshold_firing_rates(
    sorting,
    duration_in_frames=num_frames,
    threshold=firing_rate_threshold, 
    threshold_sign='less'
)

# isi violation threshold
sorting_curated = st.curation.threshold_isi_violations(
    sorting,
    duration_in_frames=num_frames,
    threshold=isi_violation_threshold, 
    threshold_sign='greater'
)

# isi violation threshold
sorting_curated = st.curation.threshold_snrs(
    sorting,
    recording=recording_processed,
    threshold=snr_threshold, 
    threshold_sign='less'
)<jupyter_output><empty_output><jupyter_text>7) Quick save to NWBTo complete the full conversion with other types of data and extended options, use a NWB Converter class<jupyter_code>from nwb_conversion_tools.utils.spike_interface import write_recording, write_sorting

output_nwb = 'blackrock_si_results.nwb'

write_recording(
    recording=rec_to_sort,
    save_path=output_nwb,
    overwrite=True
)

write_sorting(
    sorting=sorting,
    save_path=output_nwb,
    overwrite=False
)
# Check NWB file with widgets
from pynwb import NWBFile, NWBHDF5IO
from nwbwidgets import nwb2widget

io = NWBHDF5IO(output_nwb, 'r')
nwbfile = io.read()
nwb2widget(nwbfile)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09058740268931352, 'fraction_numerical': 0.014626091059212079, 'mean_word_length': 6.843663274745606, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '354633', 'license': 'MIT', 'n_tokens_mistral': 3045, 'n_tokens_neox': 2887, 'n_words': 545, 'path': 'tutorials/blackrock_spikeinterface_pipeline.ipynb', 'repo_name': 'catalystneuro/movshon-lab-to-nwb', 'size': '8478'}