julie-fabre/bombcell

Automated quality control, curation and neuron classification of spike-sorted electrophysiology data

Hi there 👋

I am PhD student in Kenneth Harris and Matteo Carandini's lab at University College London.

Here are a few of my open source ressources:

• For Neuropixels probes:
  • automated quality control bombcell
  • a curated list of awesome_neuropixels resources
  • YAHT: Yet Another Histology Tool
• For science/academia:
  • a curated list of awesome_science resources
  • a set of code to prettify MATLAB plots and scripts
  • a LaTeX PhD thesis template for UCL

bombcell's Issues

how to make a .oebin file for multiple day recordings and after cutting bad channels

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.

error with bc_loadMetricsForGUI

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

bug in 'bc_getQualityUnitType'?

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

fontname function introduced after R2022b

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

fontname in prettify_plot.m requires one more input

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.

Error in plotting quality metrics of data sorted with kilosort4

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

question about SNR calculation

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

issue with computation of timeChunks

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

Kilosort templates can be all NaNs

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

function RFV - imaginary number

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

unitType in Phy

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?

Link to prettify-matlab repo is broken

Is param.saveAsTSV already implemented?

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

bc_classifyCells needs `region` as third input argument

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

classification of MUA

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 :)

Depth relative to probe tip?

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!

Issue extracting quality metrics for data spike sorted using Kilosort4

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

another small bug

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

bug report

Hi Julie,

Some small bugs:

• function: 'isMATLABReleaseOlderThan' is missing from repo
• in 'bc_removeDuplicateSpikes' line 71 you'd need to e.g. add || ~isnan(pcFeatures) to the if statement (at least for me pcFeatures=NaN)

Cheers

Laurenz

tauC value

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.

Unrecognized function or variable 'Vrange' error while running pipeline

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.

Error in bc_waveformShape

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

bug in bc_waveformShape

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

"makepretty" is nonexistent

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.

max drift.tsv has NaN values

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?

Bug Report

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

mismatch between qcMetrics rows and number of units

classification in somatic and non-somatic units

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

Question regarding Unit IDs

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.