SpikeInterface / spikeinterface

A Python-based module for creating flexible and robust spike sorting pipelines.
https://spikeinterface.readthedocs.io
MIT License
531 stars 188 forks source link

problem of amplitude on localy exclusive and by_channel peak detection? #2319

Closed Kayv-cmb closed 11 months ago

Kayv-cmb commented 11 months ago

Hi Spikeinterface,

I have huge 'spike' in my recording more than 1.5mV which for me is higly unprobable considering electrophysiology recording. (Recorded with open_ephys and neuropixel 2.0). I wanted to have a look at those outside of the traditionnal sorter (Kilosort and stuff). However i have issue with the detect_peak function. I wrote a code (see code below) that will plot those peak by channel. However what is plot does not always have the amplitude mentioned in the detect peak (I tried both locally exclusive and by_channel) function. I filtered for <2000 uV peak but what is plotted doesnt have the amplitude expected see figure below


import numpy as np
import pandas as pd
import warnings
import tqdm 
def fxn():
    warnings.warn("deprecated", DeprecationWarning)

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    fxn()

from spikeinterface.preprocessing import common_reference,bandpass_filter
import os
import spikeinterface.full as si
import spikeinterface.extractors as se 
import spikeinterface.sorters as ss
import spikeinterface.widgets as sw
import matplotlib.pyplot as plt
# Load the preprocessed data
directory_pre = '/nemo/project/proj-schaefera-ephys/sorted_batch/NP_231120/2023-11-23_12-58-34/recording1/experiment1/'
pre_rec = si.load_extractor(directory_pre+'preprocessing_slice')
from spikeinterface.sortingcomponents.peak_detection import detect_peaks

job_kwargs = dict(chunk_duration='1s', n_jobs=40, progress_bar=True)

peaks_re = detect_peaks(
    recording=pre_rec,
    method='by_channel',
    peak_sign='neg',
    detect_threshold=5,
    exclude_sweep_ms=0.2,
    #radius_um=75,
    noise_levels=None,
    random_chunk_kwargs={},
    **job_kwargs,
)

filtered_peaks = peaks_re[peaks_re['amplitude'] < -2000]
channel_location = pre_rec.get_channel_locations()
list_channel = [] 
for i in tqdm.trange(channel_location.shape[0]):
    channel_peaks = filtered_peaks[filtered_peaks['channel_index']==i]
    plt.figure()
    for p in channel_peaks:
        channel = 'CH'+str(i)
        sample_peak = p['sample_index']
        s = int(sample_peak - 0.002*30000)
        e = int(sample_peak + 0.003*30000)
        if e>pre_rec.get_total_samples():
            e = pre_rec.get_total_samples()
        trace = pre_rec.get_traces(segment_index=0,channel_ids=[channel],start_frame=s,end_frame=e)
        plt.plot(trace)
    plt.show()

image image

but some are more or less what I expect.

image

Can I get some help to figure out what is happening, I don't understand why the plot doesnt have the amplitude I expect (more than 2000uV) ? Also many thanks for always helping me with my question !!!!

alejoe91 commented 11 months ago

HI @Kayv-cmb

To plot in uV, you need get_traces(..., return_scaled=True). Yuo are currently plotting the "raw" int16 values

Kayv-cmb commented 11 months ago

I corrected that but now all of them are below the 2000uV threshold (which is what I expect for a spike), I had but why the amplitude computed by the detect_peak still found something around 2000uV or 3000uV for some of them

alejoe91 commented 11 months ago

The peaks["amplitude"] is unscaled ;)

Kayv-cmb commented 11 months ago

I agree with that. That peaks['amplitude'] is unscaled but why when I plot with the unscaled trace the amplitude of the trace is not 2000, which is the amplitude found in the detect_peaks?

alejoe91 commented 11 months ago

I see what you mean..can you check the filtered_peaks["amplitude"] distribution?

Kayv-cmb commented 11 months ago

That's the distribution image which all of them got an unscaled amplitude below 2000

alejoe91 commented 11 months ago

Mmm the distribution is fine and I don't see anything apparently wrong in the code. Maybe @samuelgarcia can take a look?

Kayv-cmb commented 11 months ago

Also I just realized that I based one of my preprocessing method on the get_traces function but hadn't notice that it was collecting the unscale trace, so maybe the documentation can clarify this point because there is only a variable saying scale or not.

alejoe91 commented 11 months ago

Can you try to substitute this line:

channel = 'CH'+str(i)

with this:

channel = pre_rec.channel_ids[i]

Did you apply bad channel removal? This migh explain the mismatch!

Kayv-cmb commented 11 months ago

With the substitution it's working correctly I only have the filtered one why is that?

alejoe91 commented 11 months ago

so your preprocessing is only filtering? Not sure, but the second approach is better because it slices the internal list of channels.

Can you print:

print(pre_rec.channel_ids)

?

Kayv-cmb commented 11 months ago

my preprocessing is removing everything from the recording that was higher than 1500uV on the positive side because I had huge artefact, + filtering cmr and phase shift. However I based my homemade removal on the get_traces so I removed a lot of my big units because I was applying what the threshold to what I thought was scaled.

But the trace from the CH+'channel_id' is way different than the. pre_rec.channel_ids

print(pre_rec.channel_ids) that's my print ['CH1' 'CH2' 'CH3' 'CH4' 'CH5' 'CH6' 'CH7' 'CH8' 'CH9' 'CH10' 'CH11' 'CH12' 'CH13' 'CH14' 'CH15' 'CH16' 'CH17' 'CH18' 'CH19' 'CH20' 'CH21' 'CH22' 'CH23' 'CH24' 'CH25' 'CH26' 'CH27' 'CH28' 'CH29' 'CH30' 'CH31' 'CH32' 'CH33' 'CH34' 'CH35' 'CH36' 'CH37' 'CH38' 'CH39' 'CH40' 'CH41' 'CH42' 'CH43' 'CH44' 'CH45' 'CH46' 'CH47' 'CH48' 'CH49' 'CH50' 'CH51' 'CH52' 'CH53' 'CH54' 'CH55' 'CH56' 'CH57' 'CH58' 'CH59' 'CH60' 'CH61' 'CH62' 'CH63' 'CH64' 'CH65' 'CH66' 'CH67' 'CH68' 'CH69' 'CH70' 'CH71' 'CH72' 'CH73' 'CH74' 'CH75' 'CH76' 'CH77' 'CH78' 'CH79' 'CH80' 'CH81' 'CH82' 'CH83' 'CH84' 'CH85' 'CH86' 'CH87' 'CH88' 'CH89' 'CH90' 'CH91' 'CH92' 'CH93' 'CH94' 'CH95' 'CH96' 'CH97' 'CH98' 'CH99' 'CH100' 'CH101' 'CH102' 'CH103' 'CH104' 'CH105' 'CH106' 'CH107' 'CH108' 'CH109' 'CH110' 'CH111' 'CH112' 'CH113' 'CH114' 'CH115' 'CH116' 'CH117' 'CH118' 'CH119' 'CH120' 'CH121' 'CH122' 'CH123' 'CH124' 'CH125' 'CH126' 'CH127' 'CH128' 'CH129' 'CH130' 'CH131' 'CH132' 'CH133' 'CH134' 'CH135' 'CH136' 'CH137' 'CH138' 'CH139' 'CH140' 'CH141' 'CH142' 'CH143' 'CH144' 'CH145' 'CH146' 'CH147' 'CH148' 'CH149' 'CH150' 'CH151' 'CH152' 'CH153' 'CH154' 'CH155' 'CH156' 'CH157' 'CH158' 'CH159' 'CH160' 'CH161' 'CH162' 'CH163' 'CH164' 'CH165' 'CH166' 'CH167' 'CH168' 'CH169' 'CH170' 'CH171' 'CH172' 'CH173' 'CH174' 'CH175' 'CH176' 'CH177' 'CH178' 'CH179' 'CH180' 'CH181' 'CH182' 'CH183' 'CH184' 'CH185' 'CH186' 'CH187' 'CH188' 'CH189' 'CH190' 'CH191' 'CH192' 'CH193' 'CH194' 'CH195' 'CH196' 'CH197' 'CH198' 'CH199' 'CH200' 'CH201' 'CH202' 'CH203' 'CH204' 'CH205' 'CH206' 'CH207' 'CH208' 'CH209' 'CH210' 'CH211' 'CH212' 'CH213' 'CH214' 'CH215' 'CH216' 'CH217' 'CH218' 'CH219' 'CH220' 'CH221' 'CH222' 'CH223' 'CH224' 'CH225' 'CH226' 'CH227' 'CH228' 'CH229' 'CH230' 'CH231' 'CH232' 'CH233' 'CH234' 'CH235' 'CH236' ... 'CH354' 'CH355' 'CH356' 'CH357' 'CH358' 'CH359' 'CH360' 'CH361' 'CH362' 'CH363' 'CH364' 'CH365' 'CH366' 'CH367' 'CH368' 'CH369' 'CH370' 'CH371' 'CH372' 'CH373' 'CH374' 'CH375' 'CH376' 'CH377' 'CH378' 'CH379' 'CH380' 'CH381' 'CH382' 'CH383' 'CH384']

alejoe91 commented 11 months ago

Ah you see it starts from 1, not 0!

Kayv-cmb commented 11 months ago

Oh this makes a lot of sense now. Thank you for the help I will close it is solved