Merge pull request #41 from Julie-Fabre/bleeding_edge

Ephys properties and cell type classification

bc_qualityMetrics_pipeline.m

@@ -105,6 +105,10 @@
 writetable(label_table,[savePath filesep 'templates._bc_unit_labels.tsv'],'FileType', 'text','Delimiter','\t');  
 
 
+%% optional: additionally compute ephys properties for each unit 
+rerunEP = 0;
+region = ''; % options include 'Striatum' and 'Cortex'
+[ephysProperties, unitClassif] = bc_ephysPropertiesPipeline(ephysKilosortPath, savePath, rerunEP, region);
 
 
 

ephysProperties/bc_classifyCells.asv

@@ -0,0 +1,17 @@
+function unitClassif = bc_classifyCells(ephysProperties, paramEP)
+% classify striatal, GPe and cortical cells
+unitClassif = cell(size(ephysProperties,1),1);
+
+if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum, classification as in Peters et al., Nature, 2021
+    unitClassif() = 'MSN';
+    unitClassif() = 'FSI';
+    unitClassif() = 'TAN';
+    unitClassif() = 'UIN';
+
+    figure();
+    scatter3(ephysProperties.waveformDuration_peakTrough_us, ephysProperties.postSpikeSuppression_ms, ephysProperties.proplongISI);
+elseif ismember(region, {'Ctx', 'Cortical'}) % cortex, classification as in Peters et al., Cell Reports, 2022
+elseif ismember(region, {'GPe', 'Globus Pallidus external'}) %striatum
+end
+
+end

ephysProperties/bc_classifyCells.m

@@ -0,0 +1,33 @@
+function unitClassif = bc_classifyCells(ephysProperties, paramEP)
+% classify striatal, GPe and cortical cells
+unitClassif = cell(size(ephysProperties,1),1);
+
+
+if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum, classification as in Peters et al., Nature, 2021
+    unitClassif(ephysProperties.waveformDuration_peakTrough_us > paramEP.templateDuration_CP_threshold &...
+        ephysProperties.postSpikeSuppression_ms < paramEP.postSpikeSup_CP_threshold) = {'MSN'};
+
+    unitClassif(ephysProperties.waveformDuration_peakTrough_us <= paramEP.templateDuration_CP_threshold &...
+        ephysProperties.proplongISI <= paramEP.propISI_CP_threshold) = {'FSI'};
+
+    unitClassif(ephysProperties.waveformDuration_peakTrough_us > paramEP.templateDuration_CP_threshold &...
+        ephysProperties.postSpikeSuppression_ms >= paramEP.postSpikeSup_CP_threshold) = {'TAN'};
+
+    unitClassif(ephysProperties.waveformDuration_peakTrough_us <= paramEP.templateDuration_CP_threshold &...
+        ephysProperties.proplongISI > paramEP.propISI_CP_threshold) = {'UIN'};
+
+    figure();
+    scatter3(ephysProperties.waveformDuration_peakTrough_us, ephysProperties.postSpikeSuppression_ms, ephysProperties.proplongISI, 4, 'filled'); hold on;
+    set(gca, 'YDir', 'reverse' );
+
+elseif ismember(region, {'Ctx', 'Cortex', 'Cortical'}) % cortex, classification as in Peters et al., Cell Reports, 2022
+     unitClassif(ephysProperties.waveformDuration_peakTrough_us > paramEP.templateDuration_Ctx_threshold) = {'Wide-spiking'};
+
+     unitClassif(ephysProperties.waveformDuration_peakTrough_us <= paramEP.templateDuration_Ctx_threshold) = {'Narrow-spiking'};
+
+
+% elseif ismember(region, {'GPe', 'Globus Pallidus external'}) % GPe - work
+% in progress  
+end
+
+end

ephysProperties/bc_computeACGprop.m

@@ -0,0 +1,55 @@
+function [postSpikeSup_ms, tauRise_ms, tauDecay_ms, refractoryPeriod_ms] = bc_computeACGprop(thisACG, ACGbinSize, ACGduration)
+
+%% post spike suppression
+postSpikeSup = find(thisACG(500:1000) >= nanmean(thisACG(600:900))); % nanmean(ephysParams.ACG(iUnit, 900:1000)) also works.
+if ~isempty(postSpikeSup)
+    postSpikeSup = postSpikeSup(1);
+else
+    postSpikeSup = NaN;
+end
+postSpikeSup_ms = postSpikeSup  / ACGbinSize / 1000;
+
+%% tau rise and decay
+% Assuming an exponential rise, we fit the ACG to an exponential function
+% and estimate tau from the fit parameters
+ACGlags = -ACGduration/2:ACGbinSize:ACGduration/2;
+[tauRise, ~] = estimateTau(ACGlags, thisACG, true); % True for rising phase
+[tauDecay, ~] = estimateTau(ACGlags, thisACG, false); % False for decaying phase
+tauRise_ms  = tauRise   * 1000;
+tauDecay_ms = tauDecay  * 1000;
+
+%% refractory period
+% Assuming the refractory period is the time lag at which the ACG starts to rise
+ACGlags_from0 = ACGlags(ACGlags>0);
+refractoryPeriod = ACGlags_from0(find(thisACG(ACGlags>0) > (min(thisACG) + std(thisACG)), 1));
+if isempty(refractoryPeriod)
+    refractoryPeriod_ms = NaN;
+else
+    refractoryPeriod_ms = refractoryPeriod  *1000;
+end
+
+
+% figure();
+% plot(ACGlags, thisACG)
+
+%% internal functions 
+% Function to estimate tau rise/decay from ACG
+function [tau, fitResult] = estimateTau(lags, autoCorr, isRise)
+    % Fit the rising or decaying part of the ACG
+    if isRise
+        relevantPart = autoCorr(lags >= 0);
+        relevantLags = lags(lags >= 0);
+    else
+        relevantPart = autoCorr(lags <= 0);
+        relevantLags = lags(lags <= 0);
+    end
+
+    % Exponential fit
+    fitFunc = fittype('a*exp(-b*x)', 'independent', 'x');
+    [fitResult, ~] = fit(relevantLags', relevantPart', fitFunc, 'StartPoint', [max(relevantPart), 1]);
+
+    % Tau is the inverse of the b parameter in the exponential
+    tau = 1 / fitResult.b ;
+end
+
+end

ephysProperties/bc_computeAllEphysProperties.m

@@ -1,84 +1,91 @@
-function ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms_whitened, winv, paramEP, savePath)
+function [ephysProperties, unitClassif] = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms,...
+     templateAmplitudes, pcFeatures, channelPositions, paramEP, savePath)
 
 ephysProperties = struct;
-uniqueTemplates = unique(spikeTemplates);
+
+% get unit max channels
+maxChannels = bc_getWaveformMaxChannel(templateWaveforms);
+
+% extract and save or load in raw waveforms
+[rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio] = bc_extractRawWaveformsFast(paramEP, ...
+    spikeTimes_samples, spikeTemplates, paramEP.reextractRaw, savePath, paramEP.verbose); % takes ~10' for
+% an average dataset, the first time it is run, <1min after that
+
+% remove any duplicate spikes
+[uniqueTemplates, ~, spikeTimes_samples, spikeTemplates, ~, ~, ~, ~, ~, ...
+    ephysProperties.maxChannels] = ...
+    bc_removeDuplicateSpikes(spikeTimes_samples, spikeTemplates, templateAmplitudes, ...
+    pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ...
+    maxChannels, paramEP.removeDuplicateSpikes, paramEP.duplicateSpikeWindow_s, ...
+    paramEP.ephys_sample_rate, paramEP.saveSpikes_withoutDuplicates, savePath, paramEP.recomputeDuplicateSpikes);
+
 spikeTimes = spikeTimes_samples ./ paramEP.ephys_sample_rate; %convert to seconds after using sample indices to extract raw waveforms
-%timeChunks = min(spikeTimes):param.deltaTimeChunk:max(spikeTimes);
-[maxChannels, templateWaveforms] = bc_getWaveformMaxChannelEP(templateWaveforms_whitened, winv);
-%% loop through units and get ephys properties
-% QQ divide in time chunks , add plotThis 
+
+% Work in progress - divide recording into time chunks like in quality  metrics
+% spikeTimes_seconds = spikeTimes_samples ./ param.ephys_sample_rate; %convert to seconds after using sample indices to extract raw waveforms
+% if param.computeTimeChunks
+%     timeChunks = [min(spikeTimes_seconds):param.deltaTimeChunk:max(spikeTimes_seconds), max(spikeTimes_seconds)];
+% else
+%     timeChunks = [min(spikeTimes_seconds), max(spikeTimes_seconds)];
+% end
 
 fprintf('\n Extracting ephys properties ... ')
 
 for iUnit = 1:length(uniqueTemplates)
     clearvars thisUnit theseSpikeTimes theseAmplis
+
     thisUnit = uniqueTemplates(iUnit);
-    ephysProperties.clusterID(iUnit) = thisUnit;
+    ephysProperties.phy_clusterID(iUnit) = thisUnit - 1; % this is the cluster ID as it appears in phy
+    ephysProperties.clusterID(iUnit) = thisUnit; % this is the cluster ID as it appears in phy, 1-indexed (adding 1)
     theseSpikeTimes = spikeTimes(spikeTemplates == thisUnit);
 
-    %% compute ACG
+    %% ACG-based properties  
     ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis);
-
-    %% compute post spike suppression
-    ephysProperties.postSpikeSuppression(iUnit) = bc_computePSS(ephysProperties.acg(iUnit, :));
 
-    %% compute template duration
-    ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),...
-        paramEP.ephys_sample_rate);
+    [ephysProperties.postSpikeSuppression_ms(iUnit), ephysProperties.tauRise_ms(iUnit), ephysProperties.tauDecay_ms(iUnit),...
+        ephysProperties.refractoryPeriod_ms(iUnit)] = bc_computeACGprop(ephysProperties.acg(iUnit, :), paramEP.ACGbinSize, paramEP.ACGduration);
 
-    %% compute firing rate
-    ephysProperties.spike_rateSimple(iUnit) = bc_computeFR(theseSpikeTimes);
-
-    %% compute proportion long ISIs
-    ephysProperties.propLongISI(iUnit) = bc_computePropLongISI(theseSpikeTimes, paramEP.longISI);
+    %% ISI-based properties
+    ISIs = diff(spikeTimes);
 
-    %% cv, cv2
+    [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),...
+         ephysProperties.coefficient_variation2(iUnit),  ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes);
 
-    %% Fano factor
+    %% Waveform-based properties
+    % Work in progress: add option to use mean raw waveform
+    [ephysProperties.waveformDuration_peakTrough_us(iUnit), ephysProperties.halfWidth_ms(iUnit), ...
+        ephysProperties.peakTroughRatio(iUnit), ephysProperties.firstPeakTroughRatio(iUnit),...
+        ephysProperties.nPeaks(iUnit), ephysProperties.nTroughs(iUnit), ephysProperties.isSomatic(iUnit)] =...
+        bc_computeWaveformProp(templateWaveforms,thisUnit, ephysProperties.maxChannels(thisUnit),...
+        paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs);
 
-    %% skewISI
+    %% Burstiness properties
+    % Work in progress
 
-    %% max firing rate
+    %% Spike properties
+    [ephysProperties.mean_firingRate(iUnit), ephysProperties.fanoFactor(iUnit),...
+        ephysProperties.max_FiringRate(iUnit), ephysProperties.min_FiringRate(iUnit)] = bc_computeSpikeProp(theseSpikeTimes);
+
 
-    %% bursting things
+    %% Progress info
     if ((mod(iUnit, 100) == 0) || iUnit == length(uniqueTemplates)) && paramEP.verbose
        fprintf(['\n   Finished ', num2str(iUnit), ' / ', num2str(length(uniqueTemplates)), ' units.']);
     end
 end
 
 %% save ephys properties
+ephysProperties.maxChannels = ephysProperties.maxChannels(uniqueTemplates)';
+
 fprintf('\n Finished extracting ephys properties')
 try
     bc_saveEphysProperties(paramEP, ephysProperties, savePath);
     fprintf('\n Saved ephys properties to %s \n', savePath)
-    %% get some summary plots
-
+
 catch
     warning('\n Warning, ephys properties not saved! \n')
 end
-%% plot
-paramEP.plotThis=0;
-if paramEP.plotThis
-    % QQ plot histograms of each metric with the cutoffs set in params
-    figure();
-    subplot(311)
-    scatter(abs(ephysProperties.templateDuration), ephysProperties.postSpikeSuppression);
-    xlabel('waveform duration (us)')
-    ylabel('post spike suppression')
-    makepretty;
-
-    subplot(312)
-    scatter(ephysProperties.postSpikeSuppression, ephysProperties.propLongISI);
-    xlabel('post spike suppression')
-    ylabel('prop long ISI')
-    makepretty;
-
-    subplot(313)
-    scatter(abs(ephysProperties.templateDuration), ephysProperties.propLongISI);
-    xlabel('waveform duration (us)')
-    ylabel('prop long ISI')
-    makepretty;
-end
+
+%% get some summary plots - work in progress 
 
 
 end