detect_peaks and extracting waveforms without sorting

[SusanL82]

Good evening (and apologies for yet another question),

I am trying to extract spikes from axona recordings and to save the spike times and waveforms to some filetype that I can then read into matlab. There, I hope to use the Neuralynx utilities to write the spikes to an .ntt file that we can then use for manual clustering in SpikeSort3D. We want to do so that we can directly compare our manual to automatic scores, but also because the manual clustering is very fast.

I have two question about the first step:

1. When I run detect_peaks (see below), I get the following warning:

   C:\Users\susan\anaconda3\envs\SIsort\lib\site-packages\neo\rawio\axonarawio.py:316: RuntimeWarning: overflow encountered in scalar multiply
   offset = i_start // 3 * (bin_dict['bytes_packet'] // 2)

   What does this mean, and should I worry about it?

2. I would like to get waveforms on all 4 tetrode channels for the detected peaks. The waveform extractor seems to require a sorted recording. Should I try to collect small sections of my traces with get_traces instead (e.g. centered on the detected peak)?

Thanks for your help, Susan

import spikeinterface as si
import spikeinterface.extractors as se
from scipy.io import savemat
import numpy as np
from spikeinterface.sortingcomponents.peak_detection import detect_peaks
from spikeinterface.preprocessing import bandpass_filter

InFolder = "D:/SEPSIS/"
OutFolder = "D:/SEPSIS/ntt_test/"
BaseName ="HD190-1511_"

numsess = 2 #number of sessions/ .bin files
Tetnums = [2,4,5,6,7,8,9,11] #tetrodes to process
spike_thresh = 5 # detection threshold for spike detection is spike_thresh* signal SD
spike_sign = 'neg' #detect negative peaks only. can be: 'neg', 'pos', 'both'

# !!! need to manually edit the .set files so that all channels get read. 

##############################################################################

allsess = range(numsess) #count sessions zero-based

# make file list
files = list()
for f in allsess:
   thissess = allsess[f]+1 #1-based session number
   InName = InFolder+BaseName+"%02d" % thissess+".bin"
   list.append(files, InName)

# make list of SI extractor recordings
recording_list = list()
for f in allsess:
   list.append(recording_list, se.AxonaRecordingExtractor(files[f]))

# merge recording objects 
MergeRec = si.concatenate_recordings(recording_list)

# add highpass filter
MergeRec_f = bandpass_filter(MergeRec, freq_min=300.0, freq_max=6000.0, margin_ms=5.0, dtype=None)

# detect peaks per tetrode and get trace segments for waveforms
for tet in range(np.size(Tetnums)):

        thistet = Tetnums[tet]
        gettet = thistet-1
        a = 4;
        chanIDs1 = gettet*a
        chanIDs2 = chanIDs1 + 4
        chanlist = range(chanIDs1,chanIDs2)
        chanlist = ["{:01d}".format(x) for x in chanlist]

        TheseSigs = MergeRec_f.channel_slice(channel_ids = chanlist)
        TetPeaks = detect_peaks(TheseSigs, method='by_channel', pipeline_nodes = None, gather_mode='memory', folder=None, names=None, 
                                peak_sign=spike_sign, detect_threshold= spike_thresh, exclude_sweep_ms = 0.1)

[SusanL82]

I've ended up getting waveforms via get_traces, which seems to work fine (see below). I'm still curious if something similar could be done with the waveform extractor (and profit from some of the extra features that the waveform extractor has).

  allWFs = np.zeros([32,4,len(TetPeaks['sample_index'])], dtype = 'int16') 
        for p in tqdm(range(len(TetPeaks['sample_index'])), desc="collecting waveforms"):
            sf = TetPeaks['sample_index'][p] - prepeak
            ef = TetPeaks['sample_index'][p] + postpeak

            thisWF = TheseSigs.get_traces(segment_index = None, start_frame = sf, end_frame = ef)

            #write only complete spike waveforms (might skip last spike if too short)
            if np.size(thisWF,0) == 32:
                allWFs[:,:,p] = thisWF           

[samuelgarcia]

Hi Susan, your need is a corner case but this is should be doable.

Question 1: you could open an issue in https://github.com/NeuralEnsemble/python-neo Apparently we have an issue with overflow : internally some counter that could int32 should int64 or somethign similar.

Question 2 : you sedonc soultion is totally valid. We have similar machinery to do this in parralel and maybe faster in spikeinterface for this you should have a look in https://github.com/SpikeInterface/spikeinterface/blob/main/src/spikeinterface/core/waveform_tools.py

importantly : detect_peaks(TheseSigs, method='by_channel', ) you are detecting peaks by channels. If you are using tetrode or silicon probe a spike can be seen on sevral channel. To avoid multiple detection of the same spike but diffrents channel you should use detect_peaks(TheseSigs, method='locally_exclusive', radius_um=XXX )