Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

extracted peaks #116

Open
hsw28 opened this issue Dec 30, 2022 · 12 comments
Open

extracted peaks #116

hsw28 opened this issue Dec 30, 2022 · 12 comments

Comments

@hsw28
Copy link

hsw28 commented Dec 30, 2022

Hi,
I remember reading something in a readme a while ago that there is a more exact way to get peak times than using extractedPeaks, but for the life of me I cannot find it anymore. Is it deprecated?

Thank you!

@bahanonu
Copy link
Owner

Yep! Alternative peak-finding using ciapkg.signal_processing.computeSignalPeaks, see https://github.com/bahanonu/ciatah/blob/master/+ciapkg/+signal_processing/computeSignalPeaks.m.

Example usage on CNMF-e output loaded into Matlab, e.g.

  • Load file ending in _cnmfeAnalysis.mat and can use either cnmfeAnalysisOutput.extractedSignals or cnmfeAnalysisOutput.extractedSignalsEst.
  • Can vary numStdsForThresh and several other parameters to get ones that identify transients according to your criteria.
% signalPeaks: [nSignals frame] matrix. Binary matrix with 1 = peaks, 0 = non-peaks.
% signalPeaksArray: {1 nSignals} cell array. Each cell contains [1 nPeaks] vector that stores the frame locations of each peak.
% numStdsForThresh: number of standard deviations above the threshold to count as spike
% minTimeBtEvents: minimum number of time units (e.g. frames) between events.
% detectMethod:  detect on differential/1st derivative ('diff') or raw ('raw') trace.
% medianFilterLength: number of frames to calculate rolling median filter, correct for baseline drift over time.

[signalPeaks, signalPeaksArray] = ciapkg.signal_processing.computeSignalPeaks(cnmfeAnalysisOutput.extractedSignals,...
'numStdsForThresh',2,'minTimeBtEvents',8,'detectMethod','raw','medianFilterLength',100);

If you want to view the found peaks overlayed on each cell's trace quickly in a GUI, can set the makePlots to 1, e.g.

[signalPeaks, signalPeaksArray] = ciapkg.signal_processing.computeSignalPeaks(cnmfeAnalysisOutput.extractedSignals,...
'makePlots',1,...
'numStdsForThresh',2,'minTimeBtEvents',8,'detectMethod','raw','medianFilterLength',100);

image

Alternatively can load the cell extraction images/signals directly from MAT or NWB format into variables (for any cell extraction from CIAtah) using ciapkg.io.loadSignalExtraction, e.g.

cellExtractPath = 'PATH\TO\_cnmfeAnalysis.mat';

[cImages,cSignals,~,~,cSignals2] = ciapkg.io.loadSignalExtraction(cellExtractPath);

[signalPeaks, signalPeaksArray] = ciapkg.signal_processing.computeSignalPeaks(cSignals,...
'numStdsForThresh',2,'minTimeBtEvents',8,'detectMethod','raw','medianFilterLength',100);

I'll add a separate page in the CIAtah docs (https://git.io/ciatah_docs) on this.

@hmfeng
Copy link

hmfeng commented May 22, 2023

Hi @bahanonu , can you please explain more about this method? Why this computeSignalPeaks is more exact than the extracedPeaks generated by CNMFE? It looks that CNMFE uses deconvolution algorithm to compute the spikes from raw calcium trace. Did computeSignalPeaks also employ this deconvolution method? In my case, I want to see the spatial tuning of my neurons but it seems the computeSignalPeaks gives me much less spikes which brings difficulty in calculating the spatial encoding properties of the neurons.

Thank you!

@bahanonu
Copy link
Owner

bahanonu commented May 22, 2023

@hmfeng Will get back with a description in a moment, can you provide example traces where you ran ciapkg.signal_processing.computeSignalPeaks along with the exact calls/code you used? As the input parameters can have an impact on peaks identified.

@bahanonu
Copy link
Owner

bahanonu commented May 22, 2023

@hmfeng We describe the general method for that function's peak finding in the supplemental of the Corder/Ahanonu, 2019 paper, see https://www.science.org/action/downloadSupplement?doi=10.1126%2Fscience.aap8586&file=aap8586_corder_sm.pdf#page=9. Can also run on the raw trace instead of dx/dt.

See below, I'll also double check the current function as I might have updated it since then (e.g. flag to adjust peak locations to max [esp. if using 1st derivative], etc.).

image

@hmfeng
Copy link

hmfeng commented May 22, 2023

@hmfeng Will get back with a description in a moment, can you provide example traces where you ran ciapkg.signal_processing.computeSignalPeaks along with the exact calls/code you used? As the input parameters can have an impact on peaks identified.

I use the following command to compute peaks with ciatah function:

[signalPeaks, signalPeaksArray] = ciapkg.signal_processing.computeSignalPeaks(cnmfeAnalysisOutput.extractedSignals(:,1:4500),'numStdsForThresh',2.5,'minTimeBtEvents',0,'detectMethod','raw','medianFilterLength',100); % I set minTimebtEvents to 0 here just for more spikes

I compared the spikes generated by this computeSignalPeaks function and the extractedPeaks from cnmfe output dataset (only keep events with amplitude larger than 2.5*std of whole trace of this signal), as shown in below picture:

computed vs extracted peaks

@hmfeng
Copy link

hmfeng commented May 22, 2023

@hmfeng We describe the general method for that function's peak finding in the supplemental of the Corder/Ahanonu, 2019 paper, see https://www.science.org/action/downloadSupplement?doi=10.1126%2Fscience.aap8586&file=aap8586_corder_sm.pdf#page=9. Can also run on the raw trace instead of dx/dt.

See below, I'll also double check the current function as I might have updated it since then (e.g. flag to adjust peak locations to max [esp. if using 1st derivative], etc.).

image

Thank you! Actually I did refer the detailed description you mentioned in this paper as my starting point of peak extraction.

@bahanonu
Copy link
Owner

@hmfeng Will get back with a description in a moment, can you provide example traces where you ran ciapkg.signal_processing.computeSignalPeaks along with the exact calls/code you used? As the input parameters can have an impact on peaks identified.

I use the following command to compute peaks with ciatah function:

[signalPeaks, signalPeaksArray] = ciapkg.signal_processing.computeSignalPeaks(cnmfeAnalysisOutput.extractedSignals(:,1:4500),'numStdsForThresh',2.5,'minTimeBtEvents',0,'detectMethod','raw','medianFilterLength',100); % I set minTimebtEvents to 0 here just for more spikes

I compared the spikes generated by this computeSignalPeaks function and the extractedPeaks from cnmfe output dataset (only keep events with amplitude larger than 2.5*std of whole trace of this signal), as shown in below picture:

computed vs extracted peaks

Thanks, could you send over the MAT-file with the traces and cnmfe event times? Want to get a sense for whether there are actually smaller events in some of the cases where it looks like multi-events are annotated by cnmfe during the transient rise (see below).

  • Can lower medianFilterLength (the cnmfe traces you're using already have a flat baseline).
  • Given the trace is relatively clean, can reduce numStdsForThresh.

image

@hmfeng
Copy link

hmfeng commented May 27, 2023

@bahanonu sorry for the late reply. Here is the share link to the zipped data file: https://drive.google.com/file/d/1Eftew6uN-Sm8szWTbBB0g4PboxmfDNyL/view?usp=share_link
The cnmfeAnalysisoutput is the raw output and the ManualSort is the sorted tag. The trace I show in this post is from cell index 4.
Thank you!

@bahanonu
Copy link
Owner

bahanonu commented Jun 12, 2023

@hmfeng Thanks for sending the data.

Took a look at the CNMF vs. computeSignalPeaks identified transients (I change following params from your values: numStdsForThresh=1.5 and medianFilterLength=5). Looks like the CNMF events are potentially detecting noise or are also indicating "events" at each point along the rise of a given transient.

image

For example, see zoom in on the below transients where computeSignalPeaks correctly identifies the peaks but CNMF events include those along with "events" during the rise of the transient. This depends on whether you believe you can make claims about those events on the rise (e.g. are part of the burst producing the rise/transient) or if you only want to make claims based on the transient peak and amplitude.
image

image

Also, given how densely packed your cells appear to be (what region is this? CA1?), it is possible the sub-0.01 transients are actually cross-talk. I would go to cell #4 using CIAtah's computeManualSortSignals module (see https://bahanonu.github.io/ciatah/pipeline_detailed_signal_sorting_manual/) with the movie loaded so that you can check whether those smaller transients are real (you can adjust computeSignalPeaks threshold to detect those as well).

image

If you send the movie, I can take a quick look at this in terms of false positives/negatives.

@hmfeng
Copy link

hmfeng commented Oct 19, 2023

@bahanonu Thank you so much!
I found the extractedPeaks generated from cnmfe is not the real 'peaks' we are talking about. They are the 'inferred spikes' come from an algorithm OASIS (https://github.com/zhoupc/CNMF_E/tree/master/OASIS_matlab). I found many recent papers just use this inferred spikes as the quantitative measurement of the neuronal activity, especially when they are doing the place cell analysis (e.g. https://www.nature.com/articles/s41467-022-34114-x).

@bahanonu
Copy link
Owner

bahanonu commented Oct 19, 2023

@hmfeng Sure thing! Let me know if have any follow-up questions.

And that's correct, if I'm using CNMF-e I often will use the extractedSignalsEst instead of extractedSignals (within cnmfeAnalysisOutput) as depending on the situation a closer estimate of fluorescence can be more useful (even if the denoised/deconvolved OASIS traces look "nicer").

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants