From 6912812e353ea03e05061df3e221dbccff7fe6f0 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 18 Sep 2023 12:57:36 +0100 Subject: [PATCH] update function header - decompression --- bc_qualityMetrics_pipeline.asv | 94 ----------------- decompressData/bc_extractCbinData.m | 7 +- qualityMetrics/bc_getQualityUnitType.asv | 122 ----------------------- 3 files changed, 3 insertions(+), 220 deletions(-) delete mode 100644 bc_qualityMetrics_pipeline.asv delete mode 100644 qualityMetrics/bc_getQualityUnitType.asv diff --git a/bc_qualityMetrics_pipeline.asv b/bc_qualityMetrics_pipeline.asv deleted file mode 100644 index b3dd0f8..0000000 --- a/bc_qualityMetrics_pipeline.asv +++ /dev/null @@ -1,94 +0,0 @@ -%% ~~ Example bombcell pipeline ~~ -% Adjust the paths in the 'set paths' section and the parameters in bc_qualityParamValues -% This pipeline will: -% (1) load your ephys data, -% (2) decompress your raw data if it is in .cbin format -% (3) run bombcell on your data and save the output and -% (4) bring up summary plots and a GUI to flip through classified cells. -% The first time, this pipeline will be significantly slower (10-20' more) -% than after because it extracts raw waveforms. Subsequent times these -% pre-extracted waveforms are simply loaded in. -% We recommend running this pipeline on a few datasets and deciding on -% quality metric thresholds depending on the summary plots (histograms -% of the distributions of quality metrics for each unit) and GUI. - - -%% set paths - EDIT THESE -ephysKilosortPath = '/home/netshare/zinu/CB016/2021-10-07/ephys/CB016_2021-10-07_NatImages_g0/CB016_2021-10-07_NatImages_g0/pyKS/output/';% path to your kilosort output files -ephysRawDir = dir('/home/netshare/zinu/CB016/2021-10-07/ephys/CB016_2021-10-07_NatImages_g0/CB016_2021-10-07_NatImages_g0/*.*bin'); % path to yourraw .bin or .dat data -ephysMetaDir = dir('/home/netshare/zinu/CB016/2021-10-07/ephys/CB016_2021-10-07_NatImages_g0/CB016_2021-10-07_NatImages_g0/*.*meta'); % path to your .meta or .oebin meta file -saveLocation = '/media/julie/ExtraHD/CB016'; % where you want to save the quality metrics -savePath = fullfile(saveLocation, 'qMetrics'); -decompressDataLocal = '/media/julie/ExtraHD/decompressedData'; % where to save raw decompressed ephys data - -%% load data -[spikeTimes_samples, spikeTemplates, templateWaveforms, templateAmplitudes, pcFeatures, ... - pcFeatureIdx, channelPositions] = bc_loadEphysData(ephysKilosortPath); - -%% detect whether data is compressed, decompress locally if necessary -rawFile = bc_manageDataCompression(ephysRawDir, decompressDataLocal); - -%% which quality metric parameters to extract and thresholds -param = bc_qualityParamValues(ephysMetaDir, rawFile, ephysKilosortPath); -% param = bc_qualityParamValuesForUnitMatch(ephysMetaDir, rawFile) % Run this if you want to use UnitMatch after - -%% compute quality metrics -rerun = 0; -qMetricsExist = ~isempty(dir(fullfile(savePath, 'qMetric*.mat'))) || ~isempty(dir(fullfile(savePath, 'templates._bc_qMetrics.parquet'))); - -if qMetricsExist == 0 || rerun - [qMetric, unitType] = bc_runAllQualityMetrics(param, spikeTimes_samples, spikeTemplates, ... - templateWaveforms, templateAmplitudes,pcFeatures,pcFeatureIdx,channelPositions, savePath); -else - [param, qMetric] = bc_loadSavedMetrics(savePath); - unitType = bc_getQualityUnitType(param, qMetric); -end - -%% view units + quality metrics in GUI -% load data for GUI -loadRawTraces = 0; % default: don't load in raw data (this makes the GUI significantly faster) -bc_loadMetricsForGUI; - -% GUI guide: -% left/right arrow: toggle between units -% g : go to next good unit -% m : go to next multi-unit -% n : go to next noise unit -% up/down arrow: toggle between time chunks in the raw data -% u: brings up a input dialog to enter the unit you want to go to -unitQualityGuiHandle = bc_unitQualityGUI(memMapData, ephysData, qMetric, forGUI, rawWaveforms, ... - param, probeLocation, unitType, loadRawTraces); - - -%% example: get the quality metrics for one unit -% this is an example to get the quality metric for the unit with the -% original kilosort and phy label of xx (0-indexed), which corresponds to -% the unit with qMetric.clusterID == xx + 1, and to -% qMetric.phy_clusterID == xx . This is *NOT NECESSARILY* the -% (xx + 1)th row of the structure qMetric - some of the clusters that kilosort -% outputs are empty, because they were dropped in the last stages of the -% algorithm. These empty clusters are not included in the qMetric structure -% there are two ways to do this: -% 1: -original_id_we_want_to_load = 0; -id_we_want_to_load_1_indexed = original_id_we_want_to_load + 1; -number_of_spikes_for_this_cluster = qMetric.nSpikes(qMetric.clusterID == id_we_want_to_load_1_indexed); -% or 2: -original_id_we_want_to_load = 0; -number_of_spikes_for_this_cluster = qMetric.nSpikes(qMetric.phy_clusterID == original_id_we_want_to_load); - - -%% example: get unit labels -% the output of `uunitType = bc_getQualityUnitType(param, qMetric);` gives -% the unitType in a number format. 1 inidicates good units, 2 inidicates mua units, 3 -% indicates non-somatic units and 0 indicates noise units (see below) -% -goodUnits = unitType == 1; -muaUnits = unitType == 2; -noiseUnits = unitType == 0; -nonSomaticUnits = unitType == 3; - -% example: get all good units number of spikes -all_good_units_number_of_spikes = qMetric.nSpikes(goodUnits); - -% (for use with another language: output a .tsv file of labels. You can then simply load t) diff --git a/decompressData/bc_extractCbinData.m b/decompressData/bc_extractCbinData.m index e83ad52..3567e4d 100644 --- a/decompressData/bc_extractCbinData.m +++ b/decompressData/bc_extractCbinData.m @@ -1,5 +1,4 @@ function decompDataFile = bc_extractCbinData(fileName, sStartEnd, allChannelIndices, doParfor, saveFileFolder, onlySaveSyncChannel, verbose) -% adapted from a script by Micheal Krumin % % requires the zmat package: % https://github.com/fangq/zmat @@ -16,14 +15,14 @@ % using parfor regularly crashed on my computer, so I have % disabled by default. % saveFileFolder - where to save your data -% onlySaveSyncChannels - if true, only save the sync channel -% +% onlySaveSyncChannel - if true, only save the sync channel +% verbose - if true, dispply information about progress % ------ % Outputs % ------ % decompDataFile - nSamples x nChannels array of decompressed data % - +% adapted from a script by Micheal Krumin % JF: added sanity checks, more options including parfor, % save output as matrix % added size, byte, method info (zmat version used previously diff --git a/qualityMetrics/bc_getQualityUnitType.asv b/qualityMetrics/bc_getQualityUnitType.asv deleted file mode 100644 index fb0142d..0000000 --- a/qualityMetrics/bc_getQualityUnitType.asv +++ /dev/null @@ -1,122 +0,0 @@ -function unitType = bc_getQualityUnitType(param, qMetric) -% JF, Classify units into good/mua/noise/non-somatic -% ------ -% Inputs -% ------ -% -% ------ -% Outputs -% ------ -if param.computeDistanceMetrics && ~isnan(param.isoDmin) && param.computeDrift && param.extractRaw - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.signalToNoiseRatio >= param.minSNR &... - qMetric.presenceRatio >= param.minPresenceRatio & qMetric.maxDriftEstimate <= param.maxDrift & ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.isoD >= param.isoDmin &... - qMetric.Lratio <= param.lratioMax & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.computeDistanceMetrics && ~isnan(param.isoDmin) && param.computeDrift - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian' <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & ... - qMetric.presenceRatio >= param.minPresenceRatio & qMetric.maxDriftEstimate <= param.maxDrift & ... - qMetric.isoD >= param.isoDmin &... - qMetric.Lratio <= param.lratioMax & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.computeDrift && param.extractRaw - - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.signalToNoiseRatio >= param.minSNR &... - qMetric.presenceRatio >= param.minPresenceRatio & qMetric.maxDriftEstimate <= param.maxDrift & ... - qMetric.rawAmplitude > param.minAmplitude & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.computeDistanceMetrics && ~isnan(param.isoDmin) && param.extractRaw - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.signalToNoiseRatio >= param.minSNR &... - qMetric.presenceRatio >= param.minPresenceRatio & isnan(unitType)& ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.isoD >= param.isoDmin &... - qMetric.Lratio <= param.lratioMax & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.computeDrift - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian' <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & ... - qMetric.presenceRatio >= param.minPresenceRatio & qMetric.maxDriftEstimate <= param.maxDrift & ... - isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.computeDistanceMetrics && ~isnan(param.isoDmin) - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & ... - qMetric.presenceRatio >= param.minPresenceRatio & ... - qMetric.isoD >= param.isoDmin &... - qMetric.Lratio <= param.lratioMax & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -elseif param.extractRaw - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.rawAmplitude > param.minAmplitude & qMetric.signalToNoiseRatio >= param.minSNR &... - qMetric.presenceRatio >= param.minPresenceRatio & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -else - unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - unitType(qMetric.nPeaks > param.maxNPeaks | qMetric.nTroughs > param.maxNTroughs | ... - qMetric.spatialDecaySlope >= param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration |... - qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness >= param.maxWvBaselineFraction) = 0; % NOISE - unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - unitType(qMetric.percentageSpikesMissing_gaussian <= param.maxPercSpikesMissing & qMetric.nSpikes > param.minNumSpikes & ... - qMetric.fractionRPVs_estimatedTauR <= param.maxRPVviolations & ... - qMetric.presenceRatio >= param.minPresenceRatio & isnan(unitType)) = 1; % SINGLE SEXY UNIT - unitType(isnan(unitType)) = 2; % MULTI UNIT - -end