From b4a70bac541cdc55475237796a18bceafc3fb331 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Fri, 24 Nov 2023 18:31:20 +0000 Subject: [PATCH 1/8] ephys properties --- bc_qualityMetrics_pipeline.m | 4 + ephysProperties/bc_computeACGprop.m | 51 +++++++++++++ .../bc_computeAllEphysProperties.asv | 70 +++++++++++++++--- .../bc_computeAllEphysProperties.m | 73 ++++++++++++++++--- ephysProperties/bc_computeISIprop.m | 52 +++++++++++++ ephysProperties/bc_computeWaveformProp.asv | 8 ++ ephysProperties/bc_computeWaveformProp.m | 18 +++++ ephysProperties/bc_ephysPropertiesPipeline.m | 14 ++-- ephysProperties/bc_getWaveformMaxChannelEP.m | 7 +- 9 files changed, 261 insertions(+), 36 deletions(-) create mode 100644 ephysProperties/bc_computeACGprop.m create mode 100644 ephysProperties/bc_computeISIprop.m create mode 100644 ephysProperties/bc_computeWaveformProp.asv create mode 100644 ephysProperties/bc_computeWaveformProp.m diff --git a/bc_qualityMetrics_pipeline.m b/bc_qualityMetrics_pipeline.m index 852c5cdb..03f00a3d 100644 --- a/bc_qualityMetrics_pipeline.m +++ b/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 = ''; +[ephysProperties, unitClassif] = bc_ephysPropertiesPipeline(ephysKilosortPath, savePath, rerunEP, region); diff --git a/ephysProperties/bc_computeACGprop.m b/ephysProperties/bc_computeACGprop.m new file mode 100644 index 00000000..9a4479f9 --- /dev/null +++ b/ephysProperties/bc_computeACGprop.m @@ -0,0 +1,51 @@ +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)); +refractoryPeriod_ms = refractoryPeriod *1000; + + +% 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 diff --git a/ephysProperties/bc_computeAllEphysProperties.asv b/ephysProperties/bc_computeAllEphysProperties.asv index d23d3c92..0181ccea 100644 --- a/ephysProperties/bc_computeAllEphysProperties.asv +++ b/ephysProperties/bc_computeAllEphysProperties.asv @@ -1,13 +1,13 @@ -function ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms_whitened, winv, paramEP, savePath) +function [ephysProperties, unitClassif] = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, paramEP, savePath) ephysProperties = struct; uniqueTemplates = unique(spikeTemplates); 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); +maxChannels = bc_getWaveformMaxChannelEP(templateWaveforms); %% loop through units and get ephys properties -% QQ didvide in time chunks , add plotThis -paramEP.plotThis = 0; +% QQ divide in time chunks , add plotThis + fprintf('\n Extracting ephys properties ... ') for iUnit = 1:length(uniqueTemplates) @@ -16,21 +16,67 @@ for iUnit = 1:length(uniqueTemplates) ephysProperties.clusterID(iUnit) = thisUnit; theseSpikeTimes = spikeTimes(spikeTemplates == thisUnit); - %% compute ACG + %% ACG-based metrics ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); + + [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 post spike suppression - ephysProperties.postSpikeSuppression(iUnit) = bc_computePSS(ephysProperties.acg(iUnit, :)); + %% ISI-based metrics + ISIs = diff(spikeTimes); - %% compute template duration - ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... + [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),... + ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); + + %% Waveform-based metrics + bc_computeWaveformProp() + ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... paramEP.ephys_sample_rate); - %% compute firing rate + %% Burstiness metrics + + + %% compute spike metrics + % firing rate ephysProperties.spike_rateSimple(iUnit) = bc_computeFR(theseSpikeTimes); - %% compute proportion long ISIs - ephysProperties.propLongISI(iUnit) = bc_computePropLongISI(theseSpikeTimes, paramEP.longISI); + % fano factor + spikeCounts = histcounts(spikeTimes, 'BinWidth', window); + FanoFactor = var(spikeCounts) / mean(spikeCounts); + + + % Assuming 'waveform' is your waveform data vector +% And 'time' is a corresponding time vector + +% Find the peak and trough +[peakAmplitude, peakIndex] = max(waveform); +[troughAmplitude, troughIndex] = min(waveform(peakIndex:end)); +troughIndex = troughIndex + peakIndex - 1; % Adjust index + +% Compute Peak-to-Trough Duration +peakToTroughDuration = time(troughIndex) - time(peakIndex); + +% Compute Half-Width +halfAmplitude = peakAmplitude / 2; +aboveHalfIndices = find(waveform >= halfAmplitude); +halfWidthStartIndex = aboveHalfIndices(find(aboveHalfIndices < peakIndex, 1, 'last')); +halfWidthEndIndex = aboveHalfIndices(find(aboveHalfIndices > peakIndex, 1)); +halfWidth = time(halfWidthEndIndex) - time(halfWidthStartIndex); + +% Compute Rise Time +riseTime = time(peakIndex) - time(halfWidthStartIndex); + +% Compute Decay Time +decayTime = time(halfWidthEndIndex) - time(peakIndex); + +% Compute Rise Slope (Max Slope during the Rising Phase) +riseSlope = max(diff(waveform(halfWidthStartIndex:peakIndex)) ./ diff(time(halfWidthStartIndex:peakIndex))); + +% Compute Decay Slope (Max Slope during the Falling Phase) +decaySlope = min(diff(waveform(peakIndex:halfWidthEndIndex)) ./ diff(time(peakIndex:halfWidthEndIndex))); + + + %% cv, cv2 diff --git a/ephysProperties/bc_computeAllEphysProperties.m b/ephysProperties/bc_computeAllEphysProperties.m index 1d8a0d6b..f79ab203 100644 --- a/ephysProperties/bc_computeAllEphysProperties.m +++ b/ephysProperties/bc_computeAllEphysProperties.m @@ -1,10 +1,10 @@ -function ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms_whitened, winv, paramEP, savePath) +function ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, paramEP, savePath) ephysProperties = struct; uniqueTemplates = unique(spikeTemplates); 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); +maxChannels = bc_getWaveformMaxChannelEP(templateWaveforms); %% loop through units and get ephys properties % QQ divide in time chunks , add plotThis @@ -16,21 +16,76 @@ ephysProperties.clusterID(iUnit) = thisUnit; theseSpikeTimes = spikeTimes(spikeTemplates == thisUnit); - %% compute ACG + %% ACG-based metrics ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); + %units? -|> ms convert + [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 post spike suppression - ephysProperties.postSpikeSuppression(iUnit) = bc_computePSS(ephysProperties.acg(iUnit, :)); + %% ISI-based metrics + ISIs = diff(spikeTimes); - %% compute template duration - ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... + % prop long isi + bc_computePropLongISI + + % Coefficient of Variation (CV) of ISI + ISI_CV = std(ISIs) / mean(ISIs); + + % Coefficient of Variation 2 (CV2) of ISI + ISI_CV2 = 2 * mean(abs(diff(ISIs))) / mean([ISIs(1:end-1); ISIs(2:end)]); + + % ISI Skewness + ISI_Skewness = skewness(ISIs); + +% Fano Factor for Spike Counts +% Assuming 'window' is the time window over which you want to compute Fano Factor +spikeCounts = histcounts(spikeTimes, 'BinWidth', window); +FanoFactor = var(spikeCounts) / mean(spikeCounts); + + + %% Waveform-based metrics + ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... paramEP.ephys_sample_rate); + %% Burstiness metrics + + %% compute firing rate ephysProperties.spike_rateSimple(iUnit) = bc_computeFR(theseSpikeTimes); - %% compute proportion long ISIs - ephysProperties.propLongISI(iUnit) = bc_computePropLongISI(theseSpikeTimes, paramEP.longISI); + + % Assuming 'waveform' is your waveform data vector +% And 'time' is a corresponding time vector + +% Find the peak and trough +[peakAmplitude, peakIndex] = max(waveform); +[troughAmplitude, troughIndex] = min(waveform(peakIndex:end)); +troughIndex = troughIndex + peakIndex - 1; % Adjust index + +% Compute Peak-to-Trough Duration +peakToTroughDuration = time(troughIndex) - time(peakIndex); + +% Compute Half-Width +halfAmplitude = peakAmplitude / 2; +aboveHalfIndices = find(waveform >= halfAmplitude); +halfWidthStartIndex = aboveHalfIndices(find(aboveHalfIndices < peakIndex, 1, 'last')); +halfWidthEndIndex = aboveHalfIndices(find(aboveHalfIndices > peakIndex, 1)); +halfWidth = time(halfWidthEndIndex) - time(halfWidthStartIndex); + +% Compute Rise Time +riseTime = time(peakIndex) - time(halfWidthStartIndex); + +% Compute Decay Time +decayTime = time(halfWidthEndIndex) - time(peakIndex); + +% Compute Rise Slope (Max Slope during the Rising Phase) +riseSlope = max(diff(waveform(halfWidthStartIndex:peakIndex)) ./ diff(time(halfWidthStartIndex:peakIndex))); + +% Compute Decay Slope (Max Slope during the Falling Phase) +decaySlope = min(diff(waveform(peakIndex:halfWidthEndIndex)) ./ diff(time(peakIndex:halfWidthEndIndex))); + + + %% cv, cv2 diff --git a/ephysProperties/bc_computeISIprop.m b/ephysProperties/bc_computeISIprop.m new file mode 100644 index 00000000..1275165b --- /dev/null +++ b/ephysProperties/bc_computeISIprop.m @@ -0,0 +1,52 @@ +function [propLongISI, coefficient_variation, coefficient_variation2, isi_skewness] = bc_computeISIprop(ISIs, theseSpikes) + +% refs: +% Yamin/Cohen 2013 +% Stalnaker/Schoenbaum 2016 + +% (for whole session right now, can modified in the future) +spiking_stat_window = max(theseSpikes) - min(theseSpikes); +spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; + +% Get firing rate across the session +bin_spikes = ... + histcounts(theseSpikes, ... + spiking_stat_bins); + +min_spikes = 10; +use_spiking_stat_bins = bsxfun(@ge, bin_spikes, prctile(bin_spikes, 80, 2)) & bin_spikes > min_spikes; + +longISI = 2; +long_isi_total = 0; +%isi_ratios = []; +for curr_bin = find(use_spiking_stat_bins) + curr_spike_times = theseSpikes( ... + theseSpikes > spiking_stat_bins(curr_bin) & ... + theseSpikes < spiking_stat_bins(curr_bin+1)); + curr_isi = diff(curr_spike_times); + + long_isi_total = long_isi_total + sum(curr_isi(curr_isi > longISI)); + + % just compute on the first bin if there are several + if curr_bin==1 + % Coefficient of Variation (CV) of ISI + coefficient_variation = std(curr_isi) / mean(curr_isi); + + % Coefficient of Variation 2 (CV2) of ISI + coefficient_variation2 = 2 * mean(abs(diff(curr_isi))) / mean([curr_isi(1:end-1); curr_isi(2:end)]); + + % ISI Skewness + isi_skewness = skewness(ISIs); + end + + + %isi_ratios = [isi_ratios; (2 * abs(curr_isi(2:end)-curr_isi(1:end-1))) ./ ... + % (curr_isi(2:end) + curr_isi(1:end-1))]; %WRONG, see Holt 1996 +end + + +propLongISI = long_isi_total / ... + (sum(use_spiking_stat_bins(:)) * spiking_stat_window); + + +end \ No newline at end of file diff --git a/ephysProperties/bc_computeWaveformProp.asv b/ephysProperties/bc_computeWaveformProp.asv new file mode 100644 index 00000000..8034beac --- /dev/null +++ b/ephysProperties/bc_computeWaveformProp.asv @@ -0,0 +1,8 @@ +function bc_computeWaveformProp()- + +[nPeaks, nTroughs, isSomatic, peakLocs, troughLocs, waveformDuration_peakTrough, ... + spatialDecayPoints, spatialDecaySlope, ... + waveformBaseline, thisWaveform] = bc_waveformShape(templateWaveforms, ... + thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... + waveformBaselineWindow, minThreshDetectPeaksTroughs, plotThis) +end \ No newline at end of file diff --git a/ephysProperties/bc_computeWaveformProp.m b/ephysProperties/bc_computeWaveformProp.m new file mode 100644 index 00000000..40028426 --- /dev/null +++ b/ephysProperties/bc_computeWaveformProp.m @@ -0,0 +1,18 @@ +function bc_computeWaveformProp(templateWaveforms, ... + thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... + waveformBaselineWindow, minThreshDetectPeaksTroughs) + +% get waveform, waveform duration and peak/trough locations +plotThis = false; +[nPeaks, nTroughs, isSomatic, peakLocs, troughLocs, waveformDuration_peakTrough, ... + ~, ~, ~, thisWaveform] = bc_waveformShape(templateWaveforms, ... + thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... + waveformBaselineWindow, minThreshDetectPeaksTroughs, plotThis); + + +% get full width half max + +% rise time + +% asymetry +end \ No newline at end of file diff --git a/ephysProperties/bc_ephysPropertiesPipeline.m b/ephysProperties/bc_ephysPropertiesPipeline.m index 6d6055ff..48551ffe 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.m +++ b/ephysProperties/bc_ephysPropertiesPipeline.m @@ -1,26 +1,22 @@ -function ephysProperties = bc_ephysPropertiesPipeline_JF(animal, day, site, recording, experiment, rerun, runEP, region) +function ephysProperties = bc_ephysPropertiesPipeline(ephysPath, savePath, rerunEP, region) %% load ephys data paramEP = bc_ephysPropValues; %% compute ephys properties -ephysDirPath = AP_cortexlab_filenameJF(animal, day, experiment, 'ephys_dir',site, recording); -savePath = fullfile(ephysDirPath, 'ephysProperties'); ephysPropertiesExist = dir(fullfile(savePath, 'templates._bc_ephysProperties.parquet')); -if (runEP && isempty(ephysPropertiesExist)) || rerun +if isempty(ephysPropertiesExist) || rerunEP - ephysPath = AP_cortexlab_filenameJF(animal,day,experiment,'ephys',site, recording); - [spikeTimes_samples, spikeTemplates, templateWaveforms_whitened,~, ~, ~, ~] = bc_loadEphysData(ephysPath); - winv = readNPY([ephysPath filesep 'whitening_mat_inv.npy']); - ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms_whitened, winv, paramEP, savePath); + [spikeTimes_samples, spikeTemplates, templateWaveforms,~, ~, ~, ~] = bc_loadEphysData(ephysPath); + ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, paramEP, savePath); elseif ~isempty(ephysPropertiesExist) [paramEP, ephysProperties, ~] = bc_loadSavedProperties(savePath); end if ~isempty(region) - % classify cells + % classify striatal, GPe and cortical cells end end \ No newline at end of file diff --git a/ephysProperties/bc_getWaveformMaxChannelEP.m b/ephysProperties/bc_getWaveformMaxChannelEP.m index caeda8e3..412a7119 100644 --- a/ephysProperties/bc_getWaveformMaxChannelEP.m +++ b/ephysProperties/bc_getWaveformMaxChannelEP.m @@ -1,4 +1,4 @@ -function [maxChannels, templateWaveforms] = bc_getWaveformMaxChannelEP(templateWaveforms_whitened, winv) +function maxChannels = bc_getWaveformMaxChannelEP(templateWaveforms) % JF, Get the max channel for all templates (channel with largest amplitude) % ------ % Inputs @@ -11,11 +11,6 @@ % maxChannels: nTemplates × 1 vector of the channel with maximum amplitude % for each template -% Unwhiten template -templateWaveforms = nan(size(templateWaveforms_whitened)); -for t = 1:size(templateWaveforms_whitened,1) - templateWaveforms(t,:,:) = squeeze(templateWaveforms_whitened(t,:,:))*winv; -end % Get the waveform of all templates (channel with largest amplitude) [~,maxChannels] = max(max(abs(templateWaveforms),[],2),[],3); From 59171fe57290ae1b06ea427990c5b214afcad8e6 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 14:51:30 +0000 Subject: [PATCH 2/8] more ephys properties --- .../bc_computeAllEphysProperties.asv | 130 +++++---------- .../bc_computeAllEphysProperties.m | 152 +++++++----------- ephysProperties/bc_computeSpikeProp.asv | 21 +++ ephysProperties/bc_computeSpikeProp.m | 23 +++ ephysProperties/bc_computeWaveformProp.asv | 8 - ephysProperties/bc_computeWaveformProp.m | 40 ++++- ephysProperties/bc_ephysPropValues.m | 33 ++++ ephysProperties/bc_ephysPropertiesPipeline.m | 15 +- qualityMetrics/bc_waveformShape.m | 9 +- 9 files changed, 226 insertions(+), 205 deletions(-) create mode 100644 ephysProperties/bc_computeSpikeProp.asv create mode 100644 ephysProperties/bc_computeSpikeProp.m delete mode 100644 ephysProperties/bc_computeWaveformProp.asv diff --git a/ephysProperties/bc_computeAllEphysProperties.asv b/ephysProperties/bc_computeAllEphysProperties.asv index 0181ccea..be47f46a 100644 --- a/ephysProperties/bc_computeAllEphysProperties.asv +++ b/ephysProperties/bc_computeAllEphysProperties.asv @@ -1,19 +1,43 @@ -function [ephysProperties, unitClassif] = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, 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, templateAmplitudes, ... + pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... + qMetric.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 = bc_getWaveformMaxChannelEP(templateWaveforms); -%% 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); %% ACG-based metrics @@ -29,64 +53,20 @@ for iUnit = 1:length(uniqueTemplates) ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); %% Waveform-based metrics - bc_computeWaveformProp() - ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... - paramEP.ephys_sample_rate); - - %% Burstiness metrics + % 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, maxChannels(thisUnit), paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); + %% Burstiness metrics + % Work in progress %% compute spike metrics - % firing rate - ephysProperties.spike_rateSimple(iUnit) = bc_computeFR(theseSpikeTimes); - - % fano factor - spikeCounts = histcounts(spikeTimes, 'BinWidth', window); - FanoFactor = var(spikeCounts) / mean(spikeCounts); - - - % Assuming 'waveform' is your waveform data vector -% And 'time' is a corresponding time vector - -% Find the peak and trough -[peakAmplitude, peakIndex] = max(waveform); -[troughAmplitude, troughIndex] = min(waveform(peakIndex:end)); -troughIndex = troughIndex + peakIndex - 1; % Adjust index - -% Compute Peak-to-Trough Duration -peakToTroughDuration = time(troughIndex) - time(peakIndex); - -% Compute Half-Width -halfAmplitude = peakAmplitude / 2; -aboveHalfIndices = find(waveform >= halfAmplitude); -halfWidthStartIndex = aboveHalfIndices(find(aboveHalfIndices < peakIndex, 1, 'last')); -halfWidthEndIndex = aboveHalfIndices(find(aboveHalfIndices > peakIndex, 1)); -halfWidth = time(halfWidthEndIndex) - time(halfWidthStartIndex); - -% Compute Rise Time -riseTime = time(peakIndex) - time(halfWidthStartIndex); - -% Compute Decay Time -decayTime = time(halfWidthEndIndex) - time(peakIndex); - -% Compute Rise Slope (Max Slope during the Rising Phase) -riseSlope = max(diff(waveform(halfWidthStartIndex:peakIndex)) ./ diff(time(halfWidthStartIndex:peakIndex))); - -% Compute Decay Slope (Max Slope during the Falling Phase) -decaySlope = min(diff(waveform(peakIndex:halfWidthEndIndex)) ./ diff(time(peakIndex:halfWidthEndIndex))); - - - - - %% cv, cv2 - - %% Fano factor - - %% skewISI - - %% max firing rate + [ephysProperties.mean_firingRate(iUnit), FanoFactor, maxFiringRate] = 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 @@ -97,34 +77,12 @@ 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 \ No newline at end of file diff --git a/ephysProperties/bc_computeAllEphysProperties.m b/ephysProperties/bc_computeAllEphysProperties.m index f79ab203..c3a32259 100644 --- a/ephysProperties/bc_computeAllEphysProperties.m +++ b/ephysProperties/bc_computeAllEphysProperties.m @@ -1,101 +1,81 @@ -function ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, 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, templateAmplitudes, ... + pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... + qMetric.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 = bc_getWaveformMaxChannelEP(templateWaveforms); -%% 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); %% ACG-based metrics ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); - %units? -|> ms convert + [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); %% ISI-based metrics - ISIs = diff(spikeTimes); - - % prop long isi - bc_computePropLongISI + ISIs = diff(spikeTimes); - % Coefficient of Variation (CV) of ISI - ISI_CV = std(ISIs) / mean(ISIs); + [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),... + ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); - % Coefficient of Variation 2 (CV2) of ISI - ISI_CV2 = 2 * mean(abs(diff(ISIs))) / mean([ISIs(1:end-1); ISIs(2:end)]); + %% Waveform-based metrics + % Work in progress: add option to use mean raw waveform + [ephysProperties.waveformDuration_peakTrough_ms(iUnit), ephysProperties.halfWidth_ms(iUnit), ... + ephysProperties.peakTroughRatio(iUnit), ephysProperties.firstPeakTroughRatio(iUnit),... + ephysProperties.nPeaks(iUnit), ephysProperties.nTroughs(iUnit), ephysProperties.isSomatic(iUnit)] = bc_computeWaveformProp(templateWaveforms, ... + thisUnit, maxChannels(thisUnit), paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); - % ISI Skewness - ISI_Skewness = skewness(ISIs); - -% Fano Factor for Spike Counts -% Assuming 'window' is the time window over which you want to compute Fano Factor -spikeCounts = histcounts(spikeTimes, 'BinWidth', window); -FanoFactor = var(spikeCounts) / mean(spikeCounts); - - - %% Waveform-based metrics - ephysProperties.templateDuration(iUnit) = bc_computeTemplateWaveformDuration(templateWaveforms(thisUnit, :, maxChannels(iUnit)),... - paramEP.ephys_sample_rate); - %% Burstiness metrics + % Work in progress + %% compute spike metrics + % firing rate + ephysProperties.mean_firingRate(iUnit) = bc_computeFR(theseSpikeTimes); - %% compute firing rate - ephysProperties.spike_rateSimple(iUnit) = bc_computeFR(theseSpikeTimes); - - - % Assuming 'waveform' is your waveform data vector -% And 'time' is a corresponding time vector - -% Find the peak and trough -[peakAmplitude, peakIndex] = max(waveform); -[troughAmplitude, troughIndex] = min(waveform(peakIndex:end)); -troughIndex = troughIndex + peakIndex - 1; % Adjust index - -% Compute Peak-to-Trough Duration -peakToTroughDuration = time(troughIndex) - time(peakIndex); - -% Compute Half-Width -halfAmplitude = peakAmplitude / 2; -aboveHalfIndices = find(waveform >= halfAmplitude); -halfWidthStartIndex = aboveHalfIndices(find(aboveHalfIndices < peakIndex, 1, 'last')); -halfWidthEndIndex = aboveHalfIndices(find(aboveHalfIndices > peakIndex, 1)); -halfWidth = time(halfWidthEndIndex) - time(halfWidthStartIndex); - -% Compute Rise Time -riseTime = time(peakIndex) - time(halfWidthStartIndex); - -% Compute Decay Time -decayTime = time(halfWidthEndIndex) - time(peakIndex); - -% Compute Rise Slope (Max Slope during the Rising Phase) -riseSlope = max(diff(waveform(halfWidthStartIndex:peakIndex)) ./ diff(time(halfWidthStartIndex:peakIndex))); - -% Compute Decay Slope (Max Slope during the Falling Phase) -decaySlope = min(diff(waveform(peakIndex:halfWidthEndIndex)) ./ diff(time(peakIndex:halfWidthEndIndex))); - - - - - %% cv, cv2 - - %% Fano factor + % get spike counts + spikeCounts = histcounts(spikeTimes, 'BinWidth', window); + + % fano factor + FanoFactor = var(spikeCounts) / mean(spikeCounts); - %% skewISI + % max firing rate - %% max firing rate - %% bursting things + %% Progress info if ((mod(iUnit, 100) == 0) || iUnit == length(uniqueTemplates)) && paramEP.verbose fprintf(['\n Finished ', num2str(iUnit), ' / ', num2str(length(uniqueTemplates)), ' units.']); end @@ -106,34 +86,12 @@ 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 \ No newline at end of file diff --git a/ephysProperties/bc_computeSpikeProp.asv b/ephysProperties/bc_computeSpikeProp.asv new file mode 100644 index 00000000..c5b988e4 --- /dev/null +++ b/ephysProperties/bc_computeSpikeProp.asv @@ -0,0 +1,21 @@ +function [mean_firingRate, FanoFactor, max_firingRate] = bc_computeSpikeProp(theseSpikes) + spiking_stat_window = max(theseSpikes) - min(theseSpikes); + spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; + + % Get firing rate across the session + bin_spikes = histcounts(theseSpikes, spiking_stat_bins); + + min_spikes = 10; + use_spiking_stat_bins = bsxfun(@ge, bin_spikes, prctile(bin_spikes, 80, 2)) & bin_spikes > min_spikes; + + mean_firingRate = sum(bin_spikes.*use_spiking_stat_bins, 2) ./ ... + (sum(use_spiking_stat_bins, 2) * spiking_stat_window); + + % get spike counts + spikeCounts = histcounts(spikeTimes, min(theseSpikes):1:max(spikeTimes)); + + % fano factor + FanoFactor = var(spikeCounts) / mean(spikeCounts); + +end + diff --git a/ephysProperties/bc_computeSpikeProp.m b/ephysProperties/bc_computeSpikeProp.m new file mode 100644 index 00000000..74230918 --- /dev/null +++ b/ephysProperties/bc_computeSpikeProp.m @@ -0,0 +1,23 @@ +function [mean_firingRate, FanoFactor, max_firingRate] = bc_computeSpikeProp(theseSpikes) + spiking_stat_window = max(theseSpikes) - min(theseSpikes); + spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; + + % Get firing rate across the session + bin_spikes = ... + histcounts(theseSpikes, ... + spiking_stat_bins); + + min_spikes = 10; + use_spiking_stat_bins = bsxfun(@ge, bin_spikes, prctile(bin_spikes, 80, 2)) & bin_spikes > min_spikes; + + mean_firingRate = sum(bin_spikes.*use_spiking_stat_bins, 2) ./ ... + (sum(use_spiking_stat_bins, 2) * spiking_stat_window); + + % get spike counts + spikeCounts = histcounts(spikeTimes, min(spikeTimes):1:max(spikeTimes)); + + % fano factor + FanoFactor = var(spikeCounts) / mean(spikeCounts); + +end + diff --git a/ephysProperties/bc_computeWaveformProp.asv b/ephysProperties/bc_computeWaveformProp.asv deleted file mode 100644 index 8034beac..00000000 --- a/ephysProperties/bc_computeWaveformProp.asv +++ /dev/null @@ -1,8 +0,0 @@ -function bc_computeWaveformProp()- - -[nPeaks, nTroughs, isSomatic, peakLocs, troughLocs, waveformDuration_peakTrough, ... - spatialDecayPoints, spatialDecaySlope, ... - waveformBaseline, thisWaveform] = bc_waveformShape(templateWaveforms, ... - thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... - waveformBaselineWindow, minThreshDetectPeaksTroughs, plotThis) -end \ No newline at end of file diff --git a/ephysProperties/bc_computeWaveformProp.m b/ephysProperties/bc_computeWaveformProp.m index 40028426..7eeeaf8c 100644 --- a/ephysProperties/bc_computeWaveformProp.m +++ b/ephysProperties/bc_computeWaveformProp.m @@ -1,18 +1,46 @@ -function bc_computeWaveformProp(templateWaveforms, ... - thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... - waveformBaselineWindow, minThreshDetectPeaksTroughs) +function [waveformDuration_peakTrough, halfWidth, peakTroughRatio, firstPeakTroughRatio,... + nPeaks, nTroughs, isSomatic] = bc_computeWaveformProp(templateWaveforms, ... + thisUnit, maxChannel, ephys_sample_rate, channelPositions, minThreshDetectPeaksTroughs) + +% all waveform metrics based on template and not mean raw waveform for now % get waveform, waveform duration and peak/trough locations plotThis = false; +baselineThresh = NaN; +waveformBaselineWindow = NaN; [nPeaks, nTroughs, isSomatic, peakLocs, troughLocs, waveformDuration_peakTrough, ... ~, ~, ~, thisWaveform] = bc_waveformShape(templateWaveforms, ... thisUnit, maxChannel, ephys_sample_rate, channelPositions, baselineThresh, ... waveformBaselineWindow, minThreshDetectPeaksTroughs, plotThis); +% time +wvTime = 1e3 * ((0:size(thisWaveform, 2) - 1) / ephys_sample_rate); + +% Compute Half-Width +troughAmplitude = thisWaveform(troughLocs(1)); +halfAmplitude = troughAmplitude / 2; +aboveHalfIndices = find(thisWaveform >= halfAmplitude); +halfWidthStartIndex = aboveHalfIndices(find(aboveHalfIndices < peakLocs(end), 1, 'last')); +halfWidthEndIndex = aboveHalfIndices(find(aboveHalfIndices > peakLocs(end), 1)); +halfWidth = wvTime(halfWidthEndIndex) - wvTime(halfWidthStartIndex); -% get full width half max +% peak to trough ratio +peakAmplitude = thisWaveform(peakLocs(end)); +peakTroughRatio = abs(peakAmplitude/troughAmplitude); -% rise time +% 1rst peak to trough ratio +firstPeakAmplitude = max(thisWaveform(1:troughLocs(1))); +firstPeakTroughRatio = abs(firstPeakAmplitude/troughAmplitude); -% asymetry +% % Compute Rise Time +% riseTime = time(peakIndex) - time(halfWidthStartIndex); +% +% % Compute Decay Time +% decayTime = time(halfWidthEndIndex) - time(peakIndex); +% +% % Compute Rise Slope (Max Slope during the Rising Phase) +% riseSlope = max(diff(thisWaveform(halfWidthStartIndex:peakIndex)) ./ diff(time(halfWidthStartIndex:peakIndex))); +% +% % Compute Decay Slope (Max Slope during the Falling Phase) +% decaySlope = min(diff(thisWaveform(peakIndex:halfWidthEndIndex)) ./ diff(time(peakIndex:halfWidthEndIndex))); end \ No newline at end of file diff --git a/ephysProperties/bc_ephysPropValues.m b/ephysProperties/bc_ephysPropValues.m index fd62ab8b..095a86ea 100644 --- a/ephysProperties/bc_ephysPropValues.m +++ b/ephysProperties/bc_ephysPropValues.m @@ -17,13 +17,46 @@ paramEP.plotThis = 0; paramEP.verbose = 1; +% duplicate spikes parameters +paramEP.removeDuplicateSpikes = 1; +paramEP.duplicateSpikeWindow_s = 0.00001; % in seconds +paramEP.saveSpikes_withoutDuplicates = 1; +paramEP.recomputeDuplicateSpikes = 0; + +% amplitude / raw waveform parameters +paramEP.detrendWaveform = 1; % If this is set to 1, each raw extracted spike is + % detrended (we remove the best straight-fit line from the spike) + % using MATLAB's builtin function detrend. +paramEP.nRawSpikesToExtract = 100; % how many raw spikes to extract for each unit +paramEP.saveMultipleRaw = 0; % If you wish to save the nRawSpikesToExtract as well, + % currently needed if you want to run unit match https://github.com/EnnyvanBeest/UnitMatch + % to track chronic cells over days after this +paramEP.decompressData = 0; % whether to decompress .cbin ephys data +paramEP.spikeWidth = 82; % width in samples +paramEP.extractRaw = 1; %whether to extract raw waveforms or not +paramEP.probeType = 1; % if you are using spikeGLX and your meta file does + % not contain information about your probe type for some reason + % specify it here: '1' for 1.0 (3Bs) and '2' for 2.0 (single or 4-shanks) + % For additional probe types, make a pull request with more + % information. If your spikeGLX meta file contains information about your probe + % type, or if you are using open ephys, this paramater wil be ignored. +paramEP.detrendWaveforms = 0; +paramEP.reextractRaw = 0; % re extract raw waveforms or not + +% ephys properties paramEP.ephys_sample_rate = 30000; +% ACG paramEP.ACGbinSize = 0.001; paramEP.ACGduration = 1; +%Proportion Long ISI paramEP.longISI = 2; +% waveform parameters +paramEP.minThreshDetectPeaksTroughs = 0.2; % this is multiplied by the max value + % in a units waveform to give the minimum prominence to detect peaks using + % matlab's findpeaks function. % QQ set this the same as qmetrics (load useChunks Start and stop) %paramEP.computeTimeChunks = 1; % compute ephysProperties for different time chunks diff --git a/ephysProperties/bc_ephysPropertiesPipeline.m b/ephysProperties/bc_ephysPropertiesPipeline.m index 48551ffe..586b3edb 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.m +++ b/ephysProperties/bc_ephysPropertiesPipeline.m @@ -1,22 +1,25 @@ -function ephysProperties = bc_ephysPropertiesPipeline(ephysPath, savePath, rerunEP, region) +function [ephysProperties, unitClassif] = bc_ephysPropertiesPipeline(ephysPath, savePath, rerunEP, region) -%% load ephys data -paramEP = bc_ephysPropValues; %% compute ephys properties ephysPropertiesExist = dir(fullfile(savePath, 'templates._bc_ephysProperties.parquet')); if isempty(ephysPropertiesExist) || rerunEP - - [spikeTimes_samples, spikeTemplates, templateWaveforms,~, ~, ~, ~] = bc_loadEphysData(ephysPath); - ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms, paramEP, savePath); + paramEP = bc_ephysPropValues; + [spikeTimes_samples, spikeTemplates, templateWaveforms, templateAmplitudes, ... + pcFeatures, ~, channelPositions, ~] = bc_loadEphysData(ephysPath); + ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms,... + templateAmplitudes, pcFeatures, channelPositions, paramEP, savePath); elseif ~isempty(ephysPropertiesExist) [paramEP, ephysProperties, ~] = bc_loadSavedProperties(savePath); end +%% classify cells if ~isempty(region) % classify striatal, GPe and cortical cells +else + unitClassif = nan(length(),1); end end \ No newline at end of file diff --git a/qualityMetrics/bc_waveformShape.m b/qualityMetrics/bc_waveformShape.m index 949567fe..a6da4d22 100644 --- a/qualityMetrics/bc_waveformShape.m +++ b/qualityMetrics/bc_waveformShape.m @@ -153,8 +153,13 @@ spatialDecayPoints = nan(1, 6); end -waveformBaseline = max(abs(thisWaveform(waveformBaselineWindow(1): ... - waveformBaselineWindow(2)))) / max(abs(thisWaveform)); +if ~isnan(waveformBaselineWindow(1)) + waveformBaseline = max(abs(thisWaveform(waveformBaselineWindow(1): ... + waveformBaselineWindow(2)))) / max(abs(thisWaveform)); +else + waveformBaseline=NaN; +end + if plotThis From d70d50c1620c8772f3d9336638de8a67406bb719 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 16:33:14 +0000 Subject: [PATCH 3/8] debugging edge cases --- ephysProperties/bc_computeACGprop.m | 6 +++- .../bc_computeAllEphysProperties.asv | 23 +++++++------ .../bc_computeAllEphysProperties.m | 34 +++++++------------ ephysProperties/bc_computeSpikeProp.asv | 10 ++++-- ephysProperties/bc_computeSpikeProp.m | 14 ++++---- .../bc_ephysPropertiesPipeline.asv | 29 ++++++++++++++++ ephysProperties/bc_ephysPropertiesPipeline.m | 5 +++ ephysProperties/bc_saveEphysProperties.m | 16 ++++----- 8 files changed, 85 insertions(+), 52 deletions(-) create mode 100644 ephysProperties/bc_ephysPropertiesPipeline.asv diff --git a/ephysProperties/bc_computeACGprop.m b/ephysProperties/bc_computeACGprop.m index 9a4479f9..68c0e3d7 100644 --- a/ephysProperties/bc_computeACGprop.m +++ b/ephysProperties/bc_computeACGprop.m @@ -22,7 +22,11 @@ % 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)); -refractoryPeriod_ms = refractoryPeriod *1000; +if isempty(refractoryPeriod) + refractoryPeriod_ms = NaN; +else + refractoryPeriod_ms = refractoryPeriod *1000; +end % figure(); diff --git a/ephysProperties/bc_computeAllEphysProperties.asv b/ephysProperties/bc_computeAllEphysProperties.asv index be47f46a..dff97154 100644 --- a/ephysProperties/bc_computeAllEphysProperties.asv +++ b/ephysProperties/bc_computeAllEphysProperties.asv @@ -12,9 +12,8 @@ maxChannels = bc_getWaveformMaxChannel(templateWaveforms); % an average dataset, the first time it is run, <1min after that % remove any duplicate spikes -[uniqueTemplates, ~, spikeTimes_samples, spikeTemplates, templateAmplitudes, ... - pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... - qMetric.maxChannels] = ... +[uniqueTemplates, ~, spikeTimes_samples, spikeTemplates, ~, ~, ~, ~, ~, ... + ephysProperties.maxChannels] = ... bc_removeDuplicateSpikes(spikeTimes_samples, spikeTemplates, templateAmplitudes, ... pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... maxChannels, paramEP.removeDuplicateSpikes, paramEP.duplicateSpikeWindow_s, ... @@ -40,30 +39,32 @@ for iUnit = 1:length(uniqueTemplates) ephysProperties.clusterID(iUnit) = thisUnit; % this is the cluster ID as it appears in phy, 1-indexed (adding 1) theseSpikeTimes = spikeTimes(spikeTemplates == thisUnit); - %% ACG-based metrics + %% ACG-based properties ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); [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); - %% ISI-based metrics + %% ISI-based properties ISIs = diff(spikeTimes); [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),... ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); - %% Waveform-based metrics + %% 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, maxChannels(thisUnit), paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); + ephysProperties.nPeaks(iUnit), ephysProperties.nTroughs(iUnit), ephysProperties.isSomatic(iUnit)] =... + bc_computeWaveformProp(templateWaveforms,thisUnit, ephysProperties.maxChannels(thisUnit),... + paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); - %% Burstiness metrics + %% Burstiness properties % Work in progress - %% compute spike metrics - [ephysProperties.mean_firingRate(iUnit), FanoFactor, maxFiringRate] = bc_computeSpikeProp(theseSpikeTimes); + %% Spike properties + [ephysProperties.mean_firingRate(iUnit), ephysProperties.fanoFactor(iUnit),... + ephysProperties.max_FiringRate(iUnit), ephysProperties.min_FiringRate(iUnit)] = bc_computeSpikeProp(theseSpikeTimes); %% Progress info diff --git a/ephysProperties/bc_computeAllEphysProperties.m b/ephysProperties/bc_computeAllEphysProperties.m index c3a32259..dff97154 100644 --- a/ephysProperties/bc_computeAllEphysProperties.m +++ b/ephysProperties/bc_computeAllEphysProperties.m @@ -12,9 +12,8 @@ % an average dataset, the first time it is run, <1min after that % remove any duplicate spikes -[uniqueTemplates, ~, spikeTimes_samples, spikeTemplates, templateAmplitudes, ... - pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... - qMetric.maxChannels] = ... +[uniqueTemplates, ~, spikeTimes_samples, spikeTemplates, ~, ~, ~, ~, ~, ... + ephysProperties.maxChannels] = ... bc_removeDuplicateSpikes(spikeTimes_samples, spikeTemplates, templateAmplitudes, ... pcFeatures, rawWaveformsFull, rawWaveformsPeakChan, signalToNoiseRatio, ... maxChannels, paramEP.removeDuplicateSpikes, paramEP.duplicateSpikeWindow_s, ... @@ -40,40 +39,33 @@ ephysProperties.clusterID(iUnit) = thisUnit; % this is the cluster ID as it appears in phy, 1-indexed (adding 1) theseSpikeTimes = spikeTimes(spikeTemplates == thisUnit); - %% ACG-based metrics + %% ACG-based properties ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); [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); - %% ISI-based metrics + %% ISI-based properties ISIs = diff(spikeTimes); [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),... ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); - %% Waveform-based metrics + %% Waveform-based properties % Work in progress: add option to use mean raw waveform - [ephysProperties.waveformDuration_peakTrough_ms(iUnit), ephysProperties.halfWidth_ms(iUnit), ... + [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, maxChannels(thisUnit), paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); + ephysProperties.nPeaks(iUnit), ephysProperties.nTroughs(iUnit), ephysProperties.isSomatic(iUnit)] =... + bc_computeWaveformProp(templateWaveforms,thisUnit, ephysProperties.maxChannels(thisUnit),... + paramEP.ephys_sample_rate, channelPositions, paramEP.minThreshDetectPeaksTroughs); - %% Burstiness metrics + %% Burstiness properties % Work in progress - %% compute spike metrics - % firing rate - ephysProperties.mean_firingRate(iUnit) = bc_computeFR(theseSpikeTimes); - - % get spike counts - spikeCounts = histcounts(spikeTimes, 'BinWidth', window); + %% Spike properties + [ephysProperties.mean_firingRate(iUnit), ephysProperties.fanoFactor(iUnit),... + ephysProperties.max_FiringRate(iUnit), ephysProperties.min_FiringRate(iUnit)] = bc_computeSpikeProp(theseSpikeTimes); - % fano factor - FanoFactor = var(spikeCounts) / mean(spikeCounts); - - % max firing rate - %% Progress info if ((mod(iUnit, 100) == 0) || iUnit == length(uniqueTemplates)) && paramEP.verbose diff --git a/ephysProperties/bc_computeSpikeProp.asv b/ephysProperties/bc_computeSpikeProp.asv index c5b988e4..5f3acf6a 100644 --- a/ephysProperties/bc_computeSpikeProp.asv +++ b/ephysProperties/bc_computeSpikeProp.asv @@ -1,4 +1,4 @@ -function [mean_firingRate, FanoFactor, max_firingRate] = bc_computeSpikeProp(theseSpikes) +function [mean_firingRate, fanoFactor, max_firingRate, min_firingRate] = bc_computeSpikeProp(theseSpikes) spiking_stat_window = max(theseSpikes) - min(theseSpikes); spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; @@ -12,10 +12,14 @@ function [mean_firingRate, FanoFactor, max_firingRate] = bc_computeSpikeProp(the (sum(use_spiking_stat_bins, 2) * spiking_stat_window); % get spike counts - spikeCounts = histcounts(spikeTimes, min(theseSpikes):1:max(spikeTimes)); + spikeCounts = histcounts(theseSpikes, min(theseSpikes):1:max(theseSpikes)); % fano factor - FanoFactor = var(spikeCounts) / mean(spikeCounts); + fanoFactor = var(spikeCounts) / mean(spikeCounts); + + % max and min firing rate + max_firingRate = prctile(spikeCounts, 95); + min_firingRate = prctile(spikeCounts, 5); end diff --git a/ephysProperties/bc_computeSpikeProp.m b/ephysProperties/bc_computeSpikeProp.m index 74230918..5f3acf6a 100644 --- a/ephysProperties/bc_computeSpikeProp.m +++ b/ephysProperties/bc_computeSpikeProp.m @@ -1,11 +1,9 @@ -function [mean_firingRate, FanoFactor, max_firingRate] = bc_computeSpikeProp(theseSpikes) +function [mean_firingRate, fanoFactor, max_firingRate, min_firingRate] = bc_computeSpikeProp(theseSpikes) spiking_stat_window = max(theseSpikes) - min(theseSpikes); spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; % Get firing rate across the session - bin_spikes = ... - histcounts(theseSpikes, ... - spiking_stat_bins); + bin_spikes = histcounts(theseSpikes, spiking_stat_bins); min_spikes = 10; use_spiking_stat_bins = bsxfun(@ge, bin_spikes, prctile(bin_spikes, 80, 2)) & bin_spikes > min_spikes; @@ -14,10 +12,14 @@ (sum(use_spiking_stat_bins, 2) * spiking_stat_window); % get spike counts - spikeCounts = histcounts(spikeTimes, min(spikeTimes):1:max(spikeTimes)); + spikeCounts = histcounts(theseSpikes, min(theseSpikes):1:max(theseSpikes)); % fano factor - FanoFactor = var(spikeCounts) / mean(spikeCounts); + fanoFactor = var(spikeCounts) / mean(spikeCounts); + + % max and min firing rate + max_firingRate = prctile(spikeCounts, 95); + min_firingRate = prctile(spikeCounts, 5); end diff --git a/ephysProperties/bc_ephysPropertiesPipeline.asv b/ephysProperties/bc_ephysPropertiesPipeline.asv new file mode 100644 index 00000000..280b82fd --- /dev/null +++ b/ephysProperties/bc_ephysPropertiesPipeline.asv @@ -0,0 +1,29 @@ +function [ephysProperties, unitClassif] = bc_ephysPropertiesPipeline(ephysPath, savePath, rerunEP, region) + + +%% compute ephys properties +ephysPropertiesExist = dir(fullfile(savePath, 'templates._bc_ephysProperties.parquet')); + +if isempty(ephysPropertiesExist) || rerunEP + paramEP = bc_ephysPropValues; + [spikeTimes_samples, spikeTemplates, templateWaveforms, templateAmplitudes, ... + pcFeatures, ~, channelPositions, ~] = bc_loadEphysData(ephysPath); + ephysProperties = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms,... + templateAmplitudes, pcFeatures, channelPositions, paramEP, savePath); + +elseif ~isempty(ephysPropertiesExist) + [paramEP, ephysProperties, ~] = bc_loadSavedProperties(savePath); +end + +%% classify cells +if ~isempty(region) + % classify striatal, GPe and cortical cells + if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum + elseif ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum + elseif ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum + end +else + unitClassif = nan(length(),1); +end + +end \ No newline at end of file diff --git a/ephysProperties/bc_ephysPropertiesPipeline.m b/ephysProperties/bc_ephysPropertiesPipeline.m index 586b3edb..6f97e9f5 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.m +++ b/ephysProperties/bc_ephysPropertiesPipeline.m @@ -17,7 +17,12 @@ %% classify cells if ~isempty(region) + unitClassif = bc_classifyCells(ephysProperties, paramEP); % classify striatal, GPe and cortical cells + if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum + elseif ismember(region, {'Ctx', 'Cortical'}) %striatum + elseif ismember(region, {'GPe', 'Globus Pallidus external'}) %striatum + end else unitClassif = nan(length(),1); end diff --git a/ephysProperties/bc_saveEphysProperties.m b/ephysProperties/bc_saveEphysProperties.m index cd5eb79e..d82ef63e 100644 --- a/ephysProperties/bc_saveEphysProperties.m +++ b/ephysProperties/bc_saveEphysProperties.m @@ -3,20 +3,18 @@ % ------ % Inputs % ------ -% param: matlab structure defining extraction and classification parameters -% (see bc_qualityParamValues for required fields +% paramEP: matlab structure defining extraction and classification parameters +% (see bc_ephysPropValues for required fields % and suggested starting values) -% qMetric: matlab structure computed in the main loop of -% bc_runAllQualityMetrics, each field is a nUnits x 1 vector containing +% ephysProperties: matlab structure computed in the main loop of +% bc_computeAllEphysProperties, each field is a nUnits x 1 vector containing % the quality metric value for that unit -% forGUI: matlab structure computed in the main loop of bc_runAllQualityMetrics, -% for use in bc_unitQualityGUI % savePath: character array defining the path where you want to save your % quality metrics and parameters % ------ % Outputs % ------ -% ephysProperties: reformated qMetric structure into a table array +% ephysProperties: reformated ephysProperties structure into a table array if ~exist(savePath, 'dir') mkdir(fullfile(savePath)) @@ -29,9 +27,7 @@ parquetwrite([fullfile(savePath, 'templates._bc_acg.parquet')], array2table(ephysProperties.acg)) ephysProperties = rmfield(ephysProperties, 'acg'); -% save the rest of quality metrics and fraction refractory period -% violations for each unit's estimated tauR -% make sure everything is a double first +% save rest of ephys properties FNames = fieldnames(ephysProperties ); for fid = 1:length(FNames) eval(['ephysProperties.', FNames{fid}, '=double(ephysProperties.', FNames{fid}, ');']) From ff6f147cf2621b384fce165d9b84636c0b4df90e Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 16:37:23 +0000 Subject: [PATCH 4/8] more debugging edge cases --- ephysProperties/bc_computeISIprop.m | 10 +++++++++- ephysProperties/bc_computeSpikeProp.asv | 25 ------------------------- 2 files changed, 9 insertions(+), 26 deletions(-) delete mode 100644 ephysProperties/bc_computeSpikeProp.asv diff --git a/ephysProperties/bc_computeISIprop.m b/ephysProperties/bc_computeISIprop.m index 1275165b..68e653f3 100644 --- a/ephysProperties/bc_computeISIprop.m +++ b/ephysProperties/bc_computeISIprop.m @@ -43,7 +43,15 @@ %isi_ratios = [isi_ratios; (2 * abs(curr_isi(2:end)-curr_isi(1:end-1))) ./ ... % (curr_isi(2:end) + curr_isi(1:end-1))]; %WRONG, see Holt 1996 end - +if isempty(curr_bin) + coefficient_variation = NaN; + + % Coefficient of Variation 2 (CV2) of ISI + coefficient_variation2 = NaN; + + % ISI Skewness + isi_skewness = NaN; +end propLongISI = long_isi_total / ... (sum(use_spiking_stat_bins(:)) * spiking_stat_window); diff --git a/ephysProperties/bc_computeSpikeProp.asv b/ephysProperties/bc_computeSpikeProp.asv deleted file mode 100644 index 5f3acf6a..00000000 --- a/ephysProperties/bc_computeSpikeProp.asv +++ /dev/null @@ -1,25 +0,0 @@ -function [mean_firingRate, fanoFactor, max_firingRate, min_firingRate] = bc_computeSpikeProp(theseSpikes) - spiking_stat_window = max(theseSpikes) - min(theseSpikes); - spiking_stat_bins = [min(theseSpikes), max(theseSpikes)]; - - % Get firing rate across the session - bin_spikes = histcounts(theseSpikes, spiking_stat_bins); - - min_spikes = 10; - use_spiking_stat_bins = bsxfun(@ge, bin_spikes, prctile(bin_spikes, 80, 2)) & bin_spikes > min_spikes; - - mean_firingRate = sum(bin_spikes.*use_spiking_stat_bins, 2) ./ ... - (sum(use_spiking_stat_bins, 2) * spiking_stat_window); - - % get spike counts - spikeCounts = histcounts(theseSpikes, min(theseSpikes):1:max(theseSpikes)); - - % fano factor - fanoFactor = var(spikeCounts) / mean(spikeCounts); - - % max and min firing rate - max_firingRate = prctile(spikeCounts, 95); - min_firingRate = prctile(spikeCounts, 5); - -end - From 90a75e2e01d876b24c283a81c8d61a9e2e863d14 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 16:40:48 +0000 Subject: [PATCH 5/8] more debugging edge cases --- ephysProperties/bc_computeAllEphysProperties.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ephysProperties/bc_computeAllEphysProperties.m b/ephysProperties/bc_computeAllEphysProperties.m index dff97154..e5e1b18c 100644 --- a/ephysProperties/bc_computeAllEphysProperties.m +++ b/ephysProperties/bc_computeAllEphysProperties.m @@ -74,6 +74,8 @@ end %% save ephys properties +ephysProperties.maxChannels = ephysProperties.maxChannels(uniqueTemplates)'; + fprintf('\n Finished extracting ephys properties') try bc_saveEphysProperties(paramEP, ephysProperties, savePath); From 3c7a60d3417fd2cc24314ad87aa1a8d1fc74ecad Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 17:31:57 +0000 Subject: [PATCH 6/8] cell type classification for striatum and cortex --- ephysProperties/bc_classifyCells.asv | 17 ++++ ephysProperties/bc_classifyCells.m | 32 +++++++ .../bc_computeAllEphysProperties.asv | 89 ------------------- ephysProperties/bc_ephysPropValues.m | 10 ++- .../bc_ephysPropertiesPipeline.asv | 11 +-- 5 files changed, 60 insertions(+), 99 deletions(-) create mode 100644 ephysProperties/bc_classifyCells.asv create mode 100644 ephysProperties/bc_classifyCells.m delete mode 100644 ephysProperties/bc_computeAllEphysProperties.asv diff --git a/ephysProperties/bc_classifyCells.asv b/ephysProperties/bc_classifyCells.asv new file mode 100644 index 00000000..e1ac182a --- /dev/null +++ b/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 diff --git a/ephysProperties/bc_classifyCells.m b/ephysProperties/bc_classifyCells.m new file mode 100644 index 00000000..d5581073 --- /dev/null +++ b/ephysProperties/bc_classifyCells.m @@ -0,0 +1,32 @@ +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', '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 +end + +end diff --git a/ephysProperties/bc_computeAllEphysProperties.asv b/ephysProperties/bc_computeAllEphysProperties.asv deleted file mode 100644 index dff97154..00000000 --- a/ephysProperties/bc_computeAllEphysProperties.asv +++ /dev/null @@ -1,89 +0,0 @@ -function [ephysProperties, unitClassif] = bc_computeAllEphysProperties(spikeTimes_samples, spikeTemplates, templateWaveforms,... - templateAmplitudes, pcFeatures, channelPositions, paramEP, savePath) - -ephysProperties = struct; - -% 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 - -% 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.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); - - %% ACG-based properties - ephysProperties.acg(iUnit, :) = bc_computeACG(theseSpikeTimes, paramEP.ACGbinSize, paramEP.ACGduration, paramEP.plotThis); - - [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); - - %% ISI-based properties - ISIs = diff(spikeTimes); - - [ephysProperties.proplongISI(iUnit), ephysProperties.coefficient_variation(iUnit),... - ephysProperties.coefficient_variation2(iUnit), ephysProperties.isi_skewness(iUnit)] = bc_computeISIprop(ISIs, theseSpikeTimes); - - %% 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); - - %% Burstiness properties - % Work in progress - - %% Spike properties - [ephysProperties.mean_firingRate(iUnit), ephysProperties.fanoFactor(iUnit),... - ephysProperties.max_FiringRate(iUnit), ephysProperties.min_FiringRate(iUnit)] = bc_computeSpikeProp(theseSpikeTimes); - - - %% 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 -fprintf('\n Finished extracting ephys properties') -try - bc_saveEphysProperties(paramEP, ephysProperties, savePath); - fprintf('\n Saved ephys properties to %s \n', savePath) - -catch - warning('\n Warning, ephys properties not saved! \n') -end - -%% get some summary plots - work in progress - - -end \ No newline at end of file diff --git a/ephysProperties/bc_ephysPropValues.m b/ephysProperties/bc_ephysPropValues.m index 095a86ea..44b2f9ba 100644 --- a/ephysProperties/bc_ephysPropValues.m +++ b/ephysProperties/bc_ephysPropValues.m @@ -63,8 +63,12 @@ %paramEP.deltaTimeChunk = 360; %time in seconds % cell classification parameters -paramEP.propISI = 0.1; -paramEP.templateDuration = 400; -paramEP.pss = 40; +% - striatum +paramEP.propISI_CP_threshold = 0.1; +paramEP.templateDuration_CP_threshold = 400; +paramEP.postSpikeSup_CP_threshold = 40; + +% - cortex +paramEP.templateDuration_Ctx_threshold = 400; end diff --git a/ephysProperties/bc_ephysPropertiesPipeline.asv b/ephysProperties/bc_ephysPropertiesPipeline.asv index 280b82fd..bb93dc39 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.asv +++ b/ephysProperties/bc_ephysPropertiesPipeline.asv @@ -16,14 +16,11 @@ elseif ~isempty(ephysPropertiesExist) end %% classify cells -if ~isempty(region) - % classify striatal, GPe and cortical cells - if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum - elseif ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum - elseif ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum - end +if ~isempty(region) && ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS','Ctx', 'Cortical','GPe', 'Globus Pallidus external'}) + unitClassif = bc_classifyCells(ephysProperties, paramEP); + else - unitClassif = nan(length(),1); + unitClassif = nan(size(ephysProperties,1),1); end end \ No newline at end of file From bff329ea08c2e4beea8f77641e1cd14945f2f144 Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 17:34:14 +0000 Subject: [PATCH 7/8] cell type classification for striatum and cortex --- bc_qualityMetrics_pipeline.m | 2 +- ephysProperties/bc_classifyCells.m | 2 +- ephysProperties/bc_ephysPropertiesPipeline.m | 12 +++++------- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/bc_qualityMetrics_pipeline.m b/bc_qualityMetrics_pipeline.m index 03f00a3d..24ef1336 100644 --- a/bc_qualityMetrics_pipeline.m +++ b/bc_qualityMetrics_pipeline.m @@ -107,7 +107,7 @@ %% optional: additionally compute ephys properties for each unit rerunEP = 0; -region = ''; +region = ''; % options include 'Striatum' and 'Cortex' [ephysProperties, unitClassif] = bc_ephysPropertiesPipeline(ephysKilosortPath, savePath, rerunEP, region); diff --git a/ephysProperties/bc_classifyCells.m b/ephysProperties/bc_classifyCells.m index d5581073..b4c99ea2 100644 --- a/ephysProperties/bc_classifyCells.m +++ b/ephysProperties/bc_classifyCells.m @@ -20,7 +20,7 @@ scatter3(ephysProperties.waveformDuration_peakTrough_us, ephysProperties.postSpikeSuppression_ms, ephysProperties.proplongISI, 4, 'filled'); hold on; set(gca, 'YDir', 'reverse' ); -elseif ismember(region, {'Ctx', 'Cortical'}) % cortex, classification as in Peters et al., Cell Reports, 2022 +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'}; diff --git a/ephysProperties/bc_ephysPropertiesPipeline.m b/ephysProperties/bc_ephysPropertiesPipeline.m index 6f97e9f5..cc2faf9f 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.m +++ b/ephysProperties/bc_ephysPropertiesPipeline.m @@ -16,15 +16,13 @@ end %% classify cells -if ~isempty(region) +if ~isempty(region) &&... + ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS',... + 'Ctx', 'Cortical', 'Cortex', 'GPe', 'Globus Pallidus external'}) % cortex, striaum and gpe spelled every possible way unitClassif = bc_classifyCells(ephysProperties, paramEP); - % classify striatal, GPe and cortical cells - if ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS'}) %striatum - elseif ismember(region, {'Ctx', 'Cortical'}) %striatum - elseif ismember(region, {'GPe', 'Globus Pallidus external'}) %striatum - end + else - unitClassif = nan(length(),1); + unitClassif = nan(size(ephysProperties,1),1); end end \ No newline at end of file From c4b0fe26b0db4c8828b51b2b1c7d70d49b036bde Mon Sep 17 00:00:00 2001 From: Julie-Fabre Date: Mon, 27 Nov 2023 17:34:52 +0000 Subject: [PATCH 8/8] cell type classification for striatum and cortex --- ephysProperties/bc_classifyCells.m | 3 ++- ephysProperties/bc_ephysPropertiesPipeline.m | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/ephysProperties/bc_classifyCells.m b/ephysProperties/bc_classifyCells.m index b4c99ea2..da59b866 100644 --- a/ephysProperties/bc_classifyCells.m +++ b/ephysProperties/bc_classifyCells.m @@ -26,7 +26,8 @@ unitClassif(ephysProperties.waveformDuration_peakTrough_us <= paramEP.templateDuration_Ctx_threshold) = {'Narrow-spiking'}; -elseif ismember(region, {'GPe', 'Globus Pallidus external'}) % GPe +% elseif ismember(region, {'GPe', 'Globus Pallidus external'}) % GPe - work +% in progress end end diff --git a/ephysProperties/bc_ephysPropertiesPipeline.m b/ephysProperties/bc_ephysPropertiesPipeline.m index cc2faf9f..8004edd2 100644 --- a/ephysProperties/bc_ephysPropertiesPipeline.m +++ b/ephysProperties/bc_ephysPropertiesPipeline.m @@ -18,7 +18,7 @@ %% classify cells if ~isempty(region) &&... ismember(region, {'CP', 'STR', 'Striatum', 'DMS', 'DLS', 'PS',... - 'Ctx', 'Cortical', 'Cortex', 'GPe', 'Globus Pallidus external'}) % cortex, striaum and gpe spelled every possible way + 'Ctx', 'Cortical', 'Cortex'}) % cortex and striaum spelled every possible way unitClassif = bc_classifyCells(ephysProperties, paramEP); else