Merge pull request #94 from Julie-Fabre/bleeding_edge

more robust way of computing spatial decay slope

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.

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

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
 

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

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)

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

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