From 1bb65ad39e5201baadb5ea8d5a7944a8dc01221b Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Wed, 5 Jun 2024 18:14:20 +0100 Subject: [PATCH] more robust way of computing spatial decay slope --- bc_qualityMetrics_pipeline.m | 12 ++--- qualityMetrics/bc_checkParameterFields.m | 4 +- qualityMetrics/bc_qualityParamValues.m | 5 +- qualityMetrics/bc_runAllQualityMetrics.m | 2 +- qualityMetrics/bc_waveformShape.m | 4 +- qualityMetrics/helpers/bc_getSpatialDecay.asv | 53 +++++++++++++++++++ qualityMetrics/helpers/bc_getSpatialDecay.m | 13 ++++- 7 files changed, 81 insertions(+), 12 deletions(-) create mode 100644 qualityMetrics/helpers/bc_getSpatialDecay.asv diff --git a/bc_qualityMetrics_pipeline.m b/bc_qualityMetrics_pipeline.m index 6224d250..ca09af2b 100644 --- a/bc_qualityMetrics_pipeline.m +++ b/bc_qualityMetrics_pipeline.m @@ -15,14 +15,14 @@ %% set paths - EDIT THESE % '/home/netshare/zaru/JF093/2023-03-06/ephys/kilosort2/site1 -ephysKilosortPath = '/home/netshare/zaru/JF093/2023-03-06/ephys/kilosort2/site1';% path to your kilosort output files -ephysRawDir = dir('/home/netshare/zaru/JF093/2023-03-06/ephys/site1/*ap*.*bin'); % path to your raw .bin or .dat data -ephysMetaDir = dir('/home/netshare/zaru/JF093/2023-03-06/ephys/site1/*ap*.*meta'); % path to your .meta or .oebin meta file -savePath = '/media/julie/ExtraHD/JF093/2023-03-06/ephys/site1/qMetrics'; % where you want to save the quality metrics +ephysKilosortPath = '/home/netshare/zaru/AV043/2024-02-01/ephys/AV043_2024-02-01_1_ActivePassive_g0/AV043_2024-02-01_1_ActivePassive_g0_imec1/kilosort4/';% path to your kilosort output files +ephysRawDir = dir('/home/netshare/zaru/AV043/2024-02-01/ephys/AV043_2024-02-01_1_ActivePassive_g0/AV043_2024-02-01_1_ActivePassive_g0_imec1/*ap*.*bin'); % path to your raw .bin or .dat data +ephysMetaDir = dir('/home/netshare/zaru/AV043/2024-02-01/ephys/AV043_2024-02-01_1_ActivePassive_g0/AV043_2024-02-01_1_ActivePassive_g0_imec1/*ap*.*meta'); % path to your .meta or .oebin meta file +savePath = '/home/netshare/zaru/AV043/2024-02-01/ephys/AV043_2024-02-01_1_ActivePassive_g0/AV043_2024-02-01_1_ActivePassive_g0_imec1/kilosort4/qMetrics_testJF'; % where you want to save the quality metrics decompressDataLocal = '/media/julie/ExtraHD/decompressedData'; % where to save raw decompressed ephys data -gain_to_uV = 0.195; % use this if you are not using spikeGLX or openEphys to record your data. You then must leave the ephysMetaDir +gain_to_uV = NaN;%0.195; % use this if you are not using spikeGLX or openEphys to record your data. You then must leave the ephysMetaDir % empty(e.g. ephysMetaDir = '') -kilosortVersion = 2;% if using kilosort4, you need to change this value. Otherwise it does not matter. +kilosortVersion = 4;% if using kilosort4, you need to change this value. Otherwise it does not matter. %% check MATLAB version if exist('isMATLABReleaseOlderThan', 'file') == 0 % function introduced in MATLAB 2020b. diff --git a/qualityMetrics/bc_checkParameterFields.m b/qualityMetrics/bc_checkParameterFields.m index ad7e54fd..9f3b262a 100644 --- a/qualityMetrics/bc_checkParameterFields.m +++ b/qualityMetrics/bc_checkParameterFields.m @@ -42,7 +42,9 @@ % it must be at least firstPeakRatio times larger than the peak after % the trough to qualify as a non-somatic unit. 0 means this value is % not used. - +defaultValues.normalizeSpDecay = 0;% whether to normalize spatial decay points relative to +% maximum - this makes the spatrial decay slop calculation more invariant to the +% spike-sorting algorithm used %% Check for missing fields and add them with default value [param_complete, missingFields] = bc_addMissingFieldsWithDefault(param, defaultValues); diff --git a/qualityMetrics/bc_qualityParamValues.m b/qualityMetrics/bc_qualityParamValues.m index 3ab2351b..aca886af 100644 --- a/qualityMetrics/bc_qualityParamValues.m +++ b/qualityMetrics/bc_qualityParamValues.m @@ -116,6 +116,9 @@ % matlab's findpeaks function. param.firstPeakRatio = 1.1; % if units have an initial peak before the trough, % it must be at least firstPeakRatio times larger than the peak after the trough to qualify as a non-somatic unit. +param.normalizeSpDecay = 1; % whether to normalize spatial decay points relative to +% maximum - this makes the spatrial decay slop calculation more invariant to the +% spike-sorting algorithm used % recording parameters param.ephys_sample_rate = 30000; % samples per second @@ -147,7 +150,7 @@ param.somatic = 1; % keep only somatic units, and reject non-somatic ones param.minWvDuration = 100; % in us param.maxWvDuration = 1000; % in us -param.minSpatialDecaySlope = -0.003; % in a.u./um +param.minSpatialDecaySlope = -0.005; % in a.u./um param.maxWvBaselineFraction = 0.3; % maximum absolute value in waveform baseline % should not exceed this fraction of the waveform's abolute peak value diff --git a/qualityMetrics/bc_runAllQualityMetrics.m b/qualityMetrics/bc_runAllQualityMetrics.m index c4d74820..baa63033 100644 --- a/qualityMetrics/bc_runAllQualityMetrics.m +++ b/qualityMetrics/bc_runAllQualityMetrics.m @@ -159,7 +159,7 @@ forGUI.spatialDecayPoints(iUnit, :), qMetric.spatialDecaySlope(iUnit), qMetric.waveformBaselineFlatness(iUnit), ... . forGUI.tempWv(iUnit, :)] = bc_waveformShape(templateWaveforms, thisUnit, qMetric.maxChannels(thisUnit), ... param.ephys_sample_rate, channelPositions, param.maxWvBaselineFraction, waveformBaselineWindow, ... - param.minThreshDetectPeaksTroughs, param.firstPeakRatio, param.plotDetails); %do we need tempWv ? + param.minThreshDetectPeaksTroughs, param.firstPeakRatio, param.normalizeSpDecay, param.plotDetails); %do we need tempWv ? %% amplitude if param.extractRaw diff --git a/qualityMetrics/bc_waveformShape.m b/qualityMetrics/bc_waveformShape.m index 9837751a..bc664277 100644 --- a/qualityMetrics/bc_waveformShape.m +++ b/qualityMetrics/bc_waveformShape.m @@ -1,7 +1,7 @@ function [nPeaks, nTroughs, isSomatic, peakLocs, troughLocs, waveformDuration_peakTrough, ... spatialDecayPoints, spatialDecaySlope, waveformBaseline, thisWaveform] = bc_waveformShape(templateWaveforms, ... thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... - waveformBaselineWindow, minThreshDetectPeaksTroughs, firstPeakRatio, plotThis) + waveformBaselineWindow, minThreshDetectPeaksTroughs, firstPeakRatio, normalizeSpDecay, plotThis) % JF % Get the number of troughs and peaks for each waveform, % determine whether waveform is likely axonal/dendritic (biggest peak before @@ -164,7 +164,7 @@ % (get waveform spatial decay accross channels) linearFit =1; [spatialDecaySlope, spatialDecayFit, spatialDecayPoints, spatialDecayPoints_loc, estimatedUnitXY] = ... - bc_getSpatialDecay(templateWaveforms, thisUnit, maxChannel, channelPositions, linearFit); + bc_getSpatialDecay(templateWaveforms, thisUnit, maxChannel, channelPositions, linearFit, normalizeSpDecay); % (get waveform baseline fraction) diff --git a/qualityMetrics/helpers/bc_getSpatialDecay.asv b/qualityMetrics/helpers/bc_getSpatialDecay.asv new file mode 100644 index 00000000..b6d5981f --- /dev/null +++ b/qualityMetrics/helpers/bc_getSpatialDecay.asv @@ -0,0 +1,53 @@ +function [spatialDecaySlope, spatialDecayFit, spatialDecayPoints, spatialDecayPoints_loc, estimatedUnitXY] = ... + bc_getSpatialDecay(templateWaveforms, thisUnit, maxChannel, channelPositions, linearFit) + + if linearFit % linear of fit of first 6 channels (at same X position). + % In real, good units, these points decrease linearly sharply (and, in further away channels they then decrease exponentially). + % In noise artefacts they are mostly flat. + channels_withSameX = find(channelPositions(:, 1) <= channelPositions(maxChannel, 1)+33 & ... + channelPositions(:, 1) >= channelPositions(maxChannel, 1)-33); % for 4 shank probes + if numel(channels_withSameX) >= 5 + if find(channels_withSameX == maxChannel) > 5 + channels_forSpatialDecayFit = channels_withSameX( ... + find(channels_withSameX == maxChannel):-1:find(channels_withSameX == maxChannel)-5); + else + channels_forSpatialDecayFit = channels_withSameX( ... + find(channels_withSameX == maxChannel):1:min(find(channels_withSameX == maxChannel)+5, size(channels_withSameX, 1))); + end + + % get maximum value %QQ could we do value at detected trough is peak + % waveform? + spatialDecayPoints = max(abs(squeeze(templateWaveforms(thisUnit, :, channels_forSpatialDecayFit)))); + + estimatedUnitXY = channelPositions(maxChannel, :); + relativePositionsXY = channelPositions(channels_forSpatialDecayFit, :) - estimatedUnitXY; + channelPositions_relative = sqrt(nansum(relativePositionsXY.^2, 2)); + + [~, sortexChanPosIdx] = sort(channelPositions_relative); + spatialDecayPoints_norm = spatialDecayPoints(sortexChanPosIdx); + spatialDecayPoints_loc = channelPositions_relative(sortexChanPosIdx); + + % normalize spatial decay points + spatialDecayPoints_norm = spatialDecayPoints_norm ./ max(spatialDecayPoints_norm); + % linear fit + spatialDecayFit = polyfit(spatialDecayPoints_loc, spatialDecayPoints_norm', 1); % fit first order polynomial to data. first output is slope of polynomial, second is a constant + spatialDecaySlope = spatialDecayFit(1); + if length(spatialDecayPoints) < 6 + if length(spatialDecayPoints) > 1 + spatialDecayPoints = [spatialDecayPoints_norm, nan(21-length(spatialDecayPoints_norm),1)]; + else + spatialDecayPoints = [spatialDecayPoints_norm; nan(21-length(spatialDecayPoints_norm),1)]; + end + end + else + warning('No other good channels with same x location') + spatialDecayFit = NaN; + spatialDecaySlope = NaN; + spatialDecayPoints = nan(1, 6); + + end + else % not yet implemented. exponential fit? + + + end +end \ No newline at end of file diff --git a/qualityMetrics/helpers/bc_getSpatialDecay.m b/qualityMetrics/helpers/bc_getSpatialDecay.m index 84bebad8..3a6437f5 100644 --- a/qualityMetrics/helpers/bc_getSpatialDecay.m +++ b/qualityMetrics/helpers/bc_getSpatialDecay.m @@ -1,5 +1,9 @@ function [spatialDecaySlope, spatialDecayFit, spatialDecayPoints, spatialDecayPoints_loc, estimatedUnitXY] = ... - bc_getSpatialDecay(templateWaveforms, thisUnit, maxChannel, channelPositions, linearFit) + bc_getSpatialDecay(templateWaveforms, thisUnit, maxChannel, channelPositions, linearFit, normalizePoints) + +if nargin < 6 || isempty(normalizePoints) + normalizePoints = 0; +end if linearFit % linear of fit of first 6 channels (at same X position). % In real, good units, these points decrease linearly sharply (and, in further away channels they then decrease exponentially). @@ -26,6 +30,13 @@ [~, sortexChanPosIdx] = sort(channelPositions_relative); spatialDecayPoints_norm = spatialDecayPoints(sortexChanPosIdx); spatialDecayPoints_loc = channelPositions_relative(sortexChanPosIdx); + + % normalize spatial decay points + if normalizePoints + spatialDecayPoints_norm = spatialDecayPoints_norm ./ max(spatialDecayPoints_norm); + end + + % linear fit spatialDecayFit = polyfit(spatialDecayPoints_loc, spatialDecayPoints_norm', 1); % fit first order polynomial to data. first output is slope of polynomial, second is a constant spatialDecaySlope = spatialDecayFit(1); if length(spatialDecayPoints) < 6