julie-fabre / bombcell Goto Github PK
View Code? Open in Web Editor NEWAutomated quality control, curation and neuron classification of spike-sorted electrophysiology data
License: GNU General Public License v3.0
Automated quality control, curation and neuron classification of spike-sorted electrophysiology data
License: GNU General Public License v3.0
Hi! Thanks for developing this toolbox, it looks very promising. I'm running into an issue while using it. I'm getting the following error while running the bc_qualityMetrics_pipeline
:
Unrecognized function or variable 'Vrange'.
Error in bc_readSpikeGLXMetaFile (line 54)
scalingFactor = Vrange / (2 ^ bits_encoding) / gain;
Error in bc_getRawAmplitude (line 21)
scalingFactor = bc_readSpikeGLXMetaFile(metaFile);
Error in bc_runAllQualityMetrics (line 198)
qMetric.rawAmplitude(iUnit) = bc_getRawAmplitude(rawWaveformsFull(iUnit,rawWaveformsPeakChan(iUnit),:), ...
Error in bc_qualityMetrics_pipeline (line 40)
[qMetric, unitType] = bc_runAllQualityMetrics(param, spikeTimes_samples, spikeTemplates, ...
I tracked it down to the problem that my metadata file does not contain imDatPrb_type
or imProbeOpt
so it can't determine the probe type. Here is my metadata content:
acqMnMaXaDw=0,0,1,1
appVersion=20230815
fileCreateTime=2023-10-12T12:58:23
fileName=D:/NeuropixelData/450409/20231012/raw_ephys_data/_spikeGLX_ephysData_g2/_spikeGLX_ephysData_g2_t0.nidq.bin
fileSHA1=5A87AC275C00486880E99F81E94876AB03263229
fileSizeBytes=321170604
fileTimeSecs=2676.1540578326762
firstSample=265462178
gateMode=Immediate
imErrFlags_IS_CT_SR_LK_PP_SY=1 2 0 0 0 0
nDataDirs=1
nSavedChans=2
niAiRangeMax=5
niAiRangeMin=-5
niAiTermination=Default
niClockLine1=Internal
niClockSource=PXI1Slot2_2_1ch_Int : 30003.000300
niDev1=PXI1Slot2_2
niDev1ProductName=PXIe-6341
niMAChans1=
niMAGain=1
niMNChans1=
niMNGain=200
niMaxInt=32768
niMuxFactor=32
niSampRate=30003.0003
niStartEnable=false
niStartLine=PXI1Slot2_2/port0/line0
niXAChans1=0
niXDBytes1=1
niXDChans1=0:7
snsMnMaXaDw=0,0,1,1
snsSaveChanSubset=all
syncNiChan=3
syncNiChanType=0
syncNiThresh=1.1
syncSourceIdx=3
syncSourcePeriod=1
trigMode=Immediate
typeImEnabled=2
typeNiEnabled=1
typeObEnabled=0
typeThis=nidq
userNotes=
~snsChanMap=(0,0,32,1,1)(XA0;0:0)(XD0;1:1)
~snsShankMap=(1,2,0)
I'm using 1.0 probes with SpikeGLX v20230815-phase30.
Hi Julie,
First of all thanks for generating this great tool!
I am running you pipeline for the first time on my data with the example script so it might be an error on my side. The computing of the quality metrics seems to work, however when I try to start the GUI with bc_loadMetricsForGUI I get the following error:
Operands to the logical AND (&&) and OR (||) operators must be convertible to logical scalar values. Use the ANY or ALL functions to reduce operands to logical scalar values.
Error in bc_removeDuplicateSpikes (line 72)
if ~isempty(pcFeatures) || ~isnan(pcFeatures)
Error in bc_loadMetricsForGUI (line 49)
bc_removeDuplicateSpikes(ephysData.spike_times_samples, ephysData.spike_templates, ephysData.template_amplitudes,...
Error in bc_qualityMetrics_pipeline (line 62)
bc_loadMetricsForGUI;
I hope you can figure out with this information what the problem is :)
Best,
Nora
Hi Julie,
Some small bugs:
Cheers
Laurenz
Hello Julie,
I discovered that Kilosort is capable of outputting all NaNs for a template. As a result, bombcell will run into an error in "bc_waveformShape" when it tries to index into "max_waveform_value" which is empty.
My working solution is to put a try-catch around where "bc_waveformShape" is called in "bc_runAllQualityMetrics" and continue past the unit with the all NaN template.
I think it would be helpful for bombcell to be able to handle all NaN templates.
Thank you,
Luke
Hi,
Bombcell outputs all NaN values for max drift. I used kilosort3 to extract units which outputs pc features as empty variable. Is there a way to get max drift with kilosort 3 output?
First, thanks for making Bombcell. It is a great toolbox! I ran the example pipeline on my data and have some questions: How do the unit IDs relate to the ones in Phy? Can I look at both GUIs simultaneously and compare the units? Looking at the unitType variable, I see that the indexes ==1 don't relate to the somatic unit numbers in the Bombcell GUI. For example, I search for unit 11 in the GUI, and the title of the unit is a different number. And another question. I may have missed it, but which output provides the spike times per unit? I want to make rasters of the units.
Hi,
Thanks for putting togther this very useful toolbox. I have a question about the signal to noise calc in bc_extractRawWaveformsFast line 159:
abs(nanmax(squeeze(rawWaveformsFull(X,rawWaveformsPeakChan(X),:))))...
Is this meant to get the amplitude of the average waveform? If so, shouldn't it be nanmin or the abs and squeeze calls have to be swapped around? Atm this will more likely get the amplitude of the AHP of the spike which will be positive rather than the negative peak of the AP. Apologies if I am misunderstanding the SNR calculation. Cheers
Hi Julie, bombcell is great! One question - I have been interpreting depth values from 4-shank NP2.0 probes and kilosort (v2.5) as probe tip = 0 (b/c when I record from the bottom of a 4-shank I get depth values from 0-~750 but when recording from an entire shank I get values up to ~2800). I believe the depth values bombcell are the same - so, when you plot the depth in the GUI - does 0 at the top of the plot represent the tip of the probe and depth (increasing downwards) represents further along the probe (but less depth in terms of brain)? Thanks!
Hi!
First, thanks for your tool!
I use R2021b to run Kilosort (for CUDA reasons...) and in that version fontname
does not exist yet, used in prettify_plot.m
.
If deemed useful, I could do a pull request specifying the MATLAB version used in order to use, or not, this function.
Best,
Axel
Hello, I'm trying to set up bombcell on data spike-sorted with Kilosort4. However, I run into the following error when trying to compute the quality metrics:
Index exceeds the number of array elements. Index must not exceed 0.
Error in bc_upSetPlot (line 27)
pType = sPPos(dPPos);
Error in bc_upSetPlot_wrapper (line 62)
bc_upSetPlot(UpSet_data_mua, UpSet_labels_mua, figHandle_mua);
Error in bc_plotGlobalQualityMetric (line 17)
bc_upSetPlot_wrapper(qMetric, param, unitType)
Error in bc_runAllQualityMetrics (line 200)
bc_plotGlobalQualityMetric(qMetric, param, unitType, uniqueTemplates, forGUI.tempWv);
Error in bc_qualityMetrics_pipeline (line 55)
[qMetric, unitType] = bc_runAllQualityMetrics(param, spikeTimes_samples, spikeTemplates, ...
I tried with both, the matlab package available via file exchange and the cloned repo from github but the error remains.
Any idea why this is the case? Thank you so much,
Johanna
Thanks for the amazing tool. I have been trying to use and validate it on our data and I found that the max channel between the templates and the average waveform turn out to be in different places for most deeper cells. Is there a mistake in using kilosort's output templates?
I am using Npx 1.0, and KS version 2.0.
plot(squeeze(rawWaveforms.average(iCluster, :, :))') reveals there are waveforms present in the channels:
However they are not at the chansToPlot : plot(squeeze(rawWaveforms.average(iCluster, chansToPlot, :))')
Possible causes: discrepancy between KS2 template max amplitude and average waveform? or is it at the bc_extractAverageWF?
We noticed that the censored period tauC parameter is set to a default value of 0.1/1000 seconds, which corresponds to 0.1 milliseconds. We were wondering why this parameter is set to a quite low default value.
Hello,
Thank you for sharing this tool I think that it will be incredibly useful and I've already started to apply it to my own data. However, I ran into a few problems.
The first problem is that the function "makepretty" appears to be nonexistent. It's called in multiple files, but I do not believe it is declared. I was able to produce figures by commenting all "makepretty" lines out. I would appreciate it if these files could be included.
2 other minor problems: (1) the GUI at the bottom of "bc_qualityMetric_pipeline" seems to expect data for raw waveforms, but by default param.extractRaw in "bc_qualityParamValues" is set to 0 so these data don't exist. Changing the variable and rerunning solves the problem. (2) in "bc_qualityMetric_pipeline" the comment says that ephysMetaDir should have a ".meta" file but OpenEphys only outputs a ".oebin" file. Using the ".oebin" file instead works, but I had some confusion hunting for the nonexistent ".meta" file.
Hi Julie,
was going over the logical (in bc_getQualityUnitType) that defines MUA clusters and was wondering why there is an & for one of the conditions. Should all condition statements be ORs?
% Classify mua units
unitType((qMetric.percentageSpikesMissing_gaussian > param.maxPercSpikesMissing | qMetric.nSpikes < param.minNumSpikes & ...
qMetric.fractionRPVs_estimatedTauR > param.maxRPVviolations | ...
qMetric.presenceRatio < param.minPresenceRatio) & isnan(unitType)) = 2; % MUA
should it be:
% Classify mua units
unitType((qMetric.percentageSpikesMissing_gaussian > param.maxPercSpikesMissing | qMetric.nSpikes < param.minNumSpikes | ...
qMetric.fractionRPVs_estimatedTauR > param.maxRPVviolations | ...
qMetric.presenceRatio < param.minPresenceRatio) & isnan(unitType)) = 2; % MUA
Thanks :)
Hi! Thanks for this great package that you have written.
I just had a question about bc_getDistanceMetrics.
Near the beginning of the code (line 48) you check for
currentID = uniqueIDs(iID);
% Skip if current ID matches the unit of interest
if currentID == thisUnit
continue;
end
where 'thisUnit' is the ID of the unit (cluster number) and currentID
is the first column from the pc_features.npy
input. I'm kind of confused what you are comparing here. Isn't the first column from the pc_features.npy
a channel number? So you would be comparing channel number with the ID of the cluster here? Maybe I am misunderstanding something here. Any help with the interpretation would be great. Thanks a lot!
I recently started to use this very good matlab package, but I have some doubts regarding SNR and the way it is calculated. I saw that the minimum "good" value is of 0.1, that is quite low at least considering those I have always seen that are usually in the order of 3,4,5 of above, so I have the doubt that this parameter is calculated in a different way than what is classically seen, is it true? And in that case, do you have any suggestions on ranges of values that can be considered to be acceptable (other than "higher than 0.1")?
Hi Julie,
I think I found another bug. In line 153 of bc_waveformShape (spatialDecayFit = polyfit(channelPositions_relative(sortexChanPosIdx), spatialDecayPoints_norm,1)), the dimensions of the 2 inputs to polyfit don't agree. 1 is a column the other a row vector - I am using Matlab 2019 and it crashes unless I transpose one of them. I hope that helps.
Cheers
Laurenz
Hi Julie!
I think the function bc_classifyCells needs region
as third input argument.
Accordingly, where the function is called, here, would need to be changed.
Alessandro
Hi Julie,
The exporting of the quality metrics to phy works great but it's missing the classification of whether it's a good unit (unitType). This is the one I'm actually most interested in to see in phy. Could this be added?
Hi Julie,
In 'bc_getQualityUnitType' you do qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations in multiples places. Shouldn't it be qMetric.fractionRPVs <= param.maxRPVviolations. At least for me 'qMetric.fractionRPVs_estimatedTauR' doesn't get computed anymore. At least I can't see it. Thanks for your help.
Cheers
Laurenz
Hi Julie
I think that the functionality to save the tsv file to use phy after bombcell is not implemented yet, is it?
https://github.com/Julie-Fabre/bombcell/wiki/Compatibility-with-phy
Hello,
Usually there is correspondence between Kilosort's ContamPct and Bombcell's frac_RPV metrics, but occasionally there will be units where the values are completely different (see cluster IDs 10, 55, 71, and 101 in the attached screenshot). As a result, I'm not sure whether these units actually meet the criteria for inclusion. Have you noticed a discrepancy before?
Hello
This line in prettify-matlab/prettify_plot/prettify_plot.m
throws me an error because it requires one more input (I have Matlab R2022b)
Replacing fontname(options.Font)
with fontname(currFig, options.Font)
fixed it.
Hi Julie,
Sorry for all the messages. There is an issue with the way you compute the good time chunks for the quality metrics. First it seems you cannot have discontinuous time chunks. But even aside from that there is a bug. If you, e.g., have 4 time chunks and the first one is deemed bad, then currently your code will also drop the last time chunk.
Cheers
Laurenz
Hello,
I'm running into the following error when trying to extract quality metrics using the pipeline:
Error using :
Colon operands must be real scalars.
Error in bc_percSpikesMissing (line 61)
surrogate_bins = [fliplr(bins(maxAmpli_val_smooth)-bin_steps:-bin_steps:bins(maxAmpli_val_smooth) -
bin_steps*floor(size(surrogate_amplitudes,2)/2)),...
Error in bc_runAllQualityMetrics (line 118)
bc_percSpikesMissing(theseAmplis, theseSpikeTimes, timeChunks, param.plotDetails);
Error in bc_qualityMetrics_pipeline (line 51)
[qMetric, unitType] = bc_runAllQualityMetrics(param, spikeTimes_samples, spikeTemplates, ...
The extracting raw waveforms step runs fine, but this error is thrown in the subsequent step. The pipeline runs perfectly on data sorted with kilosort2.5 and all quality parameters are set to their default value. Please let me know if you have any idea what's causing this issue, or if you require any further information. Thank you!
-Cam
Hi Julie,
When using Hill's equation to calculate refractory period violations, I got sometimes imaginary numbers. I noticed that you somehow solved this issue.
I was wondering how did you get to this equation: nRPVs / (2 * (tauR(iTauR_value) - tauC) * (N_chunk - nRPVs)) and how you usually treat those units.
(I pasted the code below)
if ~isreal(fractionRPVs(iTimeChunk,iTauR_value)) % function returns imaginary number if r is too high: overestimate number.
overestimateBool(iTimeChunk,iTauR_value) = 1;
if nRPVs < N_chunk %to not get a negative wierd number or a 0 denominator
fractionRPVs(iTimeChunk,iTauR_value) = nRPVs / (2 * (tauR(iTauR_value) - tauC) * (N_chunk - nRPVs));
else
fractionRPVs(iTimeChunk,iTauR_value) = 1;
end
end
Thanks,
Sere
Hi
I am Sungmin Kang who is a post doc in Cardiff University, UK.
First of all, thanks for sharing a wonderful code.
I have recorded some ephys data using a neuropixel 1.0 probe for multiple days.
To keep same cell for days, we merged all data into one big file and run kilosort.
When I had a look for bc_qualityMetrics_pipeline.m, in line 20, I need to set the path for .oebin meta file.
However, unfortunately, as the data is merged, I feel like I need to make another .oebin file for merged data set
Do you have any idea or way to make a .oebin file for multiple recordings? or this program can only run in single data recording.
Another question is that sometimes we cut some bad channels from recorded data.
In this case, do I need to change some setting?
Thanks for again sharing your code.
Hi Julie,
I wanted to express my gratitude for the fantastic tool you've provided! It's been incredibly helpful for checking my ephys data.
While using the pipeline, I came across a potential issue in the code. In the file "bombcell/loading/bc_loadMetricsForGUI.m," specifically between lines 46 and 52, it seems that after removing duplicated spikes, the "spike_times" in the "ephysData" structure retains its original length (without removing duplicate spikes). Meanwhile, "spike_times_samples" and "spike_templates" correctly reflect the new length after the removal of duplicate spikes. This inconsistency leads to an incorrect plot of ACG and ISI in "bombcell/visualizationTools/bc_unitQualityGUI.m" (line 427). It appears that the "theseSpikeTimes" variable is calculated from "spike_times" with original spikes and "spike_templates" with removed spikes, causing the discrepancy.
Could you please verify if my observation is accurate?
Thanks!
Hui
Hi Julie,
Thank you very much for this great tool! I love it!
I am running your pipeline on my data with the example script and I get the following error:
Operands to the logical AND (&&) and OR (||) operators must be convertible to logical scalar values. Use the ANY or ALL functions to reduce operands to logical scalar values.
Error in bc_waveformShape (line 125)
if peakLoc > troughLoc && max(TRS) > max(PKS)
Error in bc_runAllQualityMetrics (line 159)
forGUI.tempWv(iUnit, :)] = bc_waveformShape(templateWaveforms, thisUnit, qMetric.maxChannels(thisUnit), ...
Error in Pipeline_code (line 49)
[qMetric, unitType] = bc_runAllQualityMetrics(param, spikeTimes_samples, spikeTemplates, ...
For additional info: peakLoc and troughLoc are empty variables while TRS and PKS are nan values. This happens for a specific unit only, but works fine on other units.
I think it might be related to issue #53 but I might be wrong. It might also just be me.
Hope it's an easy fix and thanks for your help,
Valerio
Hi Julie,
I encountered another small bug in your pipeline. In bc_saveQMetric, line 59, _fieldsRename = {'%spikes_missing', 'presence_ratio', 'max_drift', 'n_peaks', 'n_troughs', ..., '%_spikes_missing' is an invalid variable name (the '%' isn't allowed I presume).
Cheers
Laurenz
Hi Julie,
I think the way you characterise non-somatic units might be a bit too strict. See below an example where, imo, the way you have implemented it in your pipeline misclassifies the unit, I see fair amount of these in my data. Maybe this works better for your data but I thought I'd let you know.
Cheers
Laurenz
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.