From 6388ec015331a99136b7e83f982dff3607ded5de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?C=C3=A9lian=20Bimbard?= <54884591+cbimbo@users.noreply.github.com> Date: Thu, 23 Nov 2023 16:37:38 +0000 Subject: [PATCH] Add split good vs. mua for non-somatic units --- qualityMetrics/bc_checkParameterFields.m | 2 + qualityMetrics/bc_getQualityUnitType.asv | 100 ------------------ qualityMetrics/bc_getQualityUnitType.m | 22 ++-- .../bc_qualityParamValuesForUnitMatch.m | 18 +++- 4 files changed, 33 insertions(+), 109 deletions(-) delete mode 100644 qualityMetrics/bc_getQualityUnitType.asv diff --git a/qualityMetrics/bc_checkParameterFields.m b/qualityMetrics/bc_checkParameterFields.m index 824fafa1..dd4352f1 100644 --- a/qualityMetrics/bc_checkParameterFields.m +++ b/qualityMetrics/bc_checkParameterFields.m @@ -28,6 +28,8 @@ defaultValues.saveAsTSV = 0; defaultValues.unitType_for_phy = 0; +% separate good from mua in non-somatic +defaultValues.splitGoodAndMua_NonSomatic = 0; %% Check for missing fields and add them with default value [param_complete, missingFields] = bc_addMissingFieldsWithDefault(param, defaultValues); diff --git a/qualityMetrics/bc_getQualityUnitType.asv b/qualityMetrics/bc_getQualityUnitType.asv deleted file mode 100644 index 512d15c5..00000000 --- a/qualityMetrics/bc_getQualityUnitType.asv +++ /dev/null @@ -1,100 +0,0 @@ -function [unitType, unitType_string] = bc_getQualityUnitType(param, qMetric, savePath) -% JF, Classify units into good/mua/noise/non-somatic -% ------ -% Inputs -% ------ -% param: struct containing paramaters for classification -% qMetric: struct containing fields of size nUnits x 1 -% ------ -% Outputs -% ------ -% unitType - nUnits x 1 double array indicating the type of each unit: -% unitType==0 defines all noise units -% unitType==1 defines all good units -% unitType==2 defines all MUA units -% unitType==3 defines all non-somatic units -% unitType_string - nUnits x 1 string array indicating the type of each unit (good, mua, noise, non-somatic). - - -%% Sanitize and check inputs -% sanitize parameter input -param = bc_checkParameterFields(param); - -% check whether to save this classification for automated loading by phy -if nargin < 3 && param.unitType_for_phy == 1 - savePath = pwd; - warning('no save path specified. using current working directory') -end - -% Process qMetric if it's a structure and required field is not computed -if isstruct(qMetric) % if saving failed, qMetric is a structure and the fractionRPVs_estimatedTauR field we need below is not computed yet - if ~isfield('fractionRPVs_estimatedTauR', qMetric) - qMetric.fractionRPVs_estimatedTauR = arrayfun(@(x) qMetric.fractionRPVs(x, qMetric.RPV_tauR_estimate(x)), 1:size(qMetric.fractionRPVs, 1)); - qMetric = rmfield(qMetric, 'fractionRPVs'); - end -end - -% Initialize unitType array -unitType = nan(length(qMetric.percentageSpikesMissing_gaussian), 1); - -% Classify noise units -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 - -% Classify non-somatic units -unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - -% 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 - -if param.computeDistanceMetrics && ~isnan(param.isoDmin) - unitType(qMetric.isoD <= param.isoDmin & ... - qMetric.Lratio > param.lratioMax & isnan(unitType)) = 2; -end - -if param.extractRaw - unitType(qMetric.rawAmplitude <= param.minAmplitude & qMetric.signalToNoiseRatio <= param.minSNR &... - isnan(unitType)) = 2; -end - -if param.computeDrift - unitType( qMetric.maxDriftEstimate > param.maxDrift & isnan(unitType)) = 2; -end - -% Classify good units -unitType(isnan(unitType)) = 1; % SINGLE SEXY UNIT - -% Get unit type string -unitType_string = cell(size(unitType, 1), 1); -unitType_string(unitType == 0) = {'NOISE'}; -unitType_string(unitType == 3) = {'NON-SOMA'}; -unitType_string(unitType == 1) = {'GOOD'}; -unitType_string(unitType == 2) = {'MUA'}; - - -% save unitType for phy if param.unitType_for_phy is equal to 1 -try - if isfield(param, 'unitType_for_phy') % ensure back-compatibility if users have a previous version of param - if param.unitType_for_phy == 1 - if isfield(param, 'saveAsTSV') % ensure back-compatibility if users have a previous version of param - if param.saveAsTSV == 1 - cluster_id_vector = qMetric.clusterID - 1; % from bombcell to phy nomenclature - if isfield(param, 'ephysKilosortPath') - saveTSV_path = param.ephysKilosortPath; - else - saveTSV_path = savePath; - end - - cluster_table = table(cluster_id_vector, unitType_string, 'VariableNames', {'cluster_id', 'bc_unitType'}); - writetable(cluster_table, [saveTSV_path, filesep, 'cluster_bc_unitType.tsv'], 'FileType', 'text', 'Delimiter', '\t'); - end - end - end - end -catch - warning('unable to save tsv file of unitTypes') -end -end \ No newline at end of file diff --git a/qualityMetrics/bc_getQualityUnitType.m b/qualityMetrics/bc_getQualityUnitType.m index 61ab8814..6989b77f 100644 --- a/qualityMetrics/bc_getQualityUnitType.m +++ b/qualityMetrics/bc_getQualityUnitType.m @@ -17,10 +17,10 @@ %% Sanitize and check inputs -% sanitize parameter input +% Sanitize parameter input param = bc_checkParameterFields(param); -% check whether to save this classification for automated loading by phy +% Check whether to save this classification for automated loading by phy if (nargin < 3 || isempty(savePath)) && param.unitType_for_phy == 1 savePath = pwd; warning('no save path specified. using current working directory') @@ -43,9 +43,6 @@ qMetric.spatialDecaySlope > param.minSpatialDecaySlope | qMetric.waveformDuration_peakTrough < param.minWvDuration | ... qMetric.waveformDuration_peakTrough > param.maxWvDuration | qMetric.waveformBaselineFlatness > param.maxWvBaselineFraction) = 0; % NOISE -% Classify non-somatic units -unitType(qMetric.isSomatic ~= param.somatic & isnan(unitType)) = 3; % NON-SOMATIC - % Classify mua units unitType((qMetric.percentageSpikesMissing_gaussian > param.maxPercSpikesMissing | qMetric.nSpikes < param.minNumSpikes & ... qMetric.fractionRPVs_estimatedTauR > param.maxRPVviolations | ... @@ -68,12 +65,25 @@ % Classify good units unitType(isnan(unitType)) = 1; % SINGLE SEXY UNIT +% Classify non-somatic units +if param.splitGoodAndMua_NonSomatic + unitType(qMetric.isSomatic ~= param.somatic & unitType == 1) = 3; % GOOD NON-SOMATIC + unitType(qMetric.isSomatic ~= param.somatic & unitType == 2) = 4; % MUA NON-SOMATIC +else + unitType(qMetric.isSomatic ~= param.somatic & ismember(unitType,[1 2])) = 3; % NON-SOMATIC +end + % Get unit type string unitType_string = cell(size(unitType, 1), 1); unitType_string(unitType == 0) = {'NOISE'}; -unitType_string(unitType == 3) = {'NON-SOMA'}; unitType_string(unitType == 1) = {'GOOD'}; unitType_string(unitType == 2) = {'MUA'}; +if param.splitGoodAndMua_NonSomatic + unitType_string(unitType == 3) = {'NON-SOMA GOOD'}; + unitType_string(unitType == 4) = {'NON-SOMA MUA'}; +else + unitType_string(unitType == 3) = {'NON-SOMA'}; +end %% Save classification for phy % save unitType for phy if param.unitType_for_phy is equal to 1 diff --git a/qualityMetrics/bc_qualityParamValuesForUnitMatch.m b/qualityMetrics/bc_qualityParamValuesForUnitMatch.m index af0d900b..90c86067 100644 --- a/qualityMetrics/bc_qualityParamValuesForUnitMatch.m +++ b/qualityMetrics/bc_qualityParamValuesForUnitMatch.m @@ -77,14 +77,22 @@ % recorded in the raw data. This is usually 384 or 385 for neuropixels % recordings paramBC.nSyncChannels = 1; -if ~isempty(ephysMetaDir) +if exist('ephysMetaDir','var') && ~isempty(ephysMetaDir) paramBC.ephysMetaFile = [ephysMetaDir.folder, filesep, ephysMetaDir.name]; paramBC.gain_to_uV = NaN; else paramBC.ephysMetaFile = 'NaN'; - paramBC.gain_to_uV = gain_to_uV; + if exist('gain_to_uV','var') + paramBC.gain_to_uV = gain_to_uV; + else + paramBC.gain_to_uV = NaN; + end +end +if exist('rawFile','var') + paramBC.rawFile = rawFile; +else + paramBC.rawFile = 'NaN'; end -paramBC.rawFile = rawFile; % distance metric parameters paramBC.computeDistanceMetrics = 0; % whether to compute distance metrics - this can be time consuming @@ -115,6 +123,10 @@ paramBC.isoDmin = 20; paramBC.lratioMax = 0.1; paramBC.ssMin = NaN; + +% split good and mua non-somatic +paramBC.splitGoodAndMua_NonSomatic = 1; + end % %bc_qualityParamValues % BCparam = struct;