option to output .tsv files for phy

bc_qualityMetrics_pipeline.m

@@ -29,7 +29,7 @@
 rawFile = bc_manageDataCompression(ephysRawDir, decompressDataLocal);
 
 %% which quality metric parameters to extract and thresholds 
-param = bc_qualityParamValues(ephysMetaDir, rawFile); 
+param = bc_qualityParamValues(ephysMetaDir, rawFile, ephysKilosortPath); 
 % param = bc_qualityParamValuesForUnitMatch(ephysMetaDir, rawFile) % Run this if you want to use UnitMatch after
 
 %% compute quality metrics 

qualityMetrics/bc_qualityParamValues.asv

@@ -0,0 +1,123 @@
+function param = bc_qualityParamValues(ephysMetaDir, rawFile, ephysKilosortPath)
+% JF, Load a parameter structure defining extraction and
+% classification parameters
+% ------
+% Inputs
+% ------
+% ephysMetaDir: dir() structure of the path to your .meta or .oebin meta
+%   file
+% rawFile: character array defining the path where your uncompressed raw
+%   ephys data is
+% ------
+% Outputs
+% ------
+% param: matlab structure defining extraction and
+% classification parameters (see bc_qualityParamValues for required fields
+% and suggested starting values)
+% 
+
+param = struct; %initialize structure 
+
+%% calculating quality metrics parameters 
+% plotting parameters 
+param.plotDetails = 0; % generates a lot of plots, 
+% mainly good if you running through the code line by line to check things,
+% to debug, or to get nice plots for a presentation
+param.plotGlobal = 1; % plot summary of quality metrics 
+param.verbose = 1; % update user on progress
+param.reextractRaw = 0; % re extract raw waveforms or not 
+
+% saving parameters 
+param.saveAsTSV = 1; % additionally save outputs at .tsv file - this is 
+% useful if you want to use phy after bombcell: each quality metric value
+% will appear as a column in the Cluster view
+if nargin < 3
+    warning('no ephys kilosort path defined in bc_qualityParamValues, will save output tsv file in the save P')
+else
+    param.ephysKilosortPath = ephysKilosortPath;
+end
+param.saveMatFileForGUI = 1; % save certain outputs at .mat file - useful for GUI
+
+% amplitude parameters
+param.nRawSpikesToExtract = 100; % how many raw spikes to extract for each unit 
+param.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
+param.decompressData = 0; % whether to decompress .cbin ephys data 
+param.spikeWidth = 82; % width in samples 
+param.extractRaw = 1; %whether to extract raw waveforms or not 
+
+% signal to noise ratio
+param.waveformBaselineNoiseWindow = 20; %time in samples at beginning of times
+% extracted to computer the mean raw waveform - this needs to be before the
+% waveform starts 
+
+% refractory period parameters
+param.tauR_valuesMin = 2/1000; % refractory period time (s), usually 0.0020. 
+% If this value is different than param.tauR_valuesMax, bombcell will
+% estimate the tauR value taking possible values between :
+% param.tauR_valuesMin:param.tauR_valuesStep:param.tauR_valuesMax
+param.tauR_valuesStep = 0.5/1000; % refractory period time (s) steps. Only 
+% used if param.tauR_valuesMin is different from param.tauR_valuesMax
+param.tauR_valuesMax = 2/1000; % refractory period time (s), usually 0.0020
+param.tauC = 0.1/1000; % censored period time (s)
+
+% percentage spikes missing parameters 
+param.computeTimeChunks = 1; % compute fraction refractory period violations 
+% and percent sp[ikes missing for different time chunks 
+param.deltaTimeChunk = 360; %time in seconds 
+
+% presence ratio 
+param.presenceRatioBinSize = 60; % in seconds 
+
+% drift estimate
+param.driftBinSize = 60; % in seconds
+param.computeDrift = 0; % whether to compute each units drift. this is a 
+% critically slow step that takes around 2seconds per unit 
+
+% waveform parameters
+param.waveformBaselineWindowStart = 20;
+param.waveformBaselineWindowStop = 30; % in samples 
+param.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.
+
+% recording parametrs
+param.ephys_sample_rate = 30000; % samples per second
+param.nChannels = 385; %number of recorded channels (including any sync channels)
+% recorded in the raw data. This is usually 384 or 385 for neuropixels
+% recordings
+param.nSyncChannels = 1;
+param.ephysMetaFile = [ephysMetaDir.folder, filesep, ephysMetaDir.name];
+param.rawFile = rawFile;
+
+% distance metric parameters
+param.computeDistanceMetrics = 0; % whether to compute distance metrics - this can be time consuming 
+param.nChannelsIsoDist = 4; % number of nearby channels to use in distance metric computation 
+
+
+%% classifying units into good/mua/noise parameters 
+param.minAmplitude = 20; % in uV
+param.maxRPVviolations = 0.1; % fraction
+param.maxPercSpikesMissing = 20; % in percentage
+param.minNumSpikes = 300; % number of spikes
+
+param.maxDrift = 100;
+param.minPresenceRatio = 0.7;
+param.minSNR = 0.1;
+
+%waveform 
+param.maxNPeaks = 2; % maximum number of peaks
+param.maxNTroughs = 1; % maximum number of troughs
+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 V/um
+param.maxWvBaselineFraction = 0.3; % maximum absolute value in waveform baseline
+% should not exceed this fraction of the waveform's abolute peak value
+
+%distance metrics
+param.isoDmin = 20; % minimum isolation distance value
+param.lratioMax = 0.1; % maximum l-ratio value
+param.ssMin = NaN; % minimum silhouette score 
+end

qualityMetrics/bc_qualityParamValues.m

@@ -1,4 +1,4 @@
-function param = bc_qualityParamValues(ephysMetaDir, rawFile)
+function param = bc_qualityParamValues(ephysMetaDir, rawFile, ephysKilosortPath)
 % JF, Load a parameter structure defining extraction and
 % classification parameters
 % ------
@@ -28,7 +28,14 @@
 param.reextractRaw = 0; % re extract raw waveforms or not 
 
 % saving parameters 
-param.saveAsParquet = 1; % save outputs at .parquet file 
+param.saveAsTSV = 1; % additionally save outputs at .tsv file - this is 
+% useful if you want to use phy after bombcell: each quality metric value
+% will appear as a column in the Cluster view
+if nargin < 3
+    warning('no ephys kilosort path defined in bc_qualityParamValues, will save output tsv file in the savePath location')
+else
+    param.ephysKilosortPath = ephysKilosortPath;
+end
 param.saveMatFileForGUI = 1; % save certain outputs at .mat file - useful for GUI
 
 % amplitude parameters

qualityMetrics/bc_saveQMetrics.asv

@@ -0,0 +1,71 @@
+function qMetric = bc_saveQMetrics(param, qMetric, forGUI, savePath)
+% JF, Reformat and save ephys properties
+% ------
+% Inputs
+% ------
+% paramEP: matlab structure defining extraction and classification parameters 
+%   (see bc_ephysProperties Values for required fields
+%   and suggested starting values)
+% ephysProperties: matlab structure computed in the main loop of
+%   bc_computeAllEphysProperties
+% savePath: character array defining the path where you want to save your
+%   quality metrics and parameters 
+% ------
+% Outputs
+% ------
+% ephysProperties: reformated ephysProperties structure into a table array
+
+if ~exist(savePath, 'dir')
+    mkdir(fullfile(savePath))
+end
+
+% save parameters
+if ~istable(param)
+    parquetwrite([fullfile(savePath, '_bc_parameters._bc_qMetrics.parquet')], struct2table(param))
+end
+% save quality metrics
+if param.saveMatFileForGUI
+    save(fullfile(savePath, 'templates.qualityMetricDetailsforGUI.mat'), 'forGUI', '-v7.3')
+end
+
+% save fraction refractory period violations for all different tauR times
+parquetwrite([fullfile(savePath, 'templates._bc_fractionRefractoryPeriodViolationsPerTauR.parquet')], array2table(qMetric.fractionRPVs))
+qMetric.fractionRPVs_estimatedTauR = arrayfun(@(x) qMetric.fractionRPVs(x, qMetric.RPV_tauR_estimate(x)), 1:size(qMetric.fractionRPVs,1));
+qMetric = rmfield(qMetric, 'fractionRPVs');
+
+% save the rest of quality metrics and fraction refractory period
+% violations for each unit's estimated tauR
+% make sure everything is a double first
+FNames = fieldnames(qMetric);
+for fid = 1:length(FNames)
+    eval(['qMetric.', FNames{fid}, '=double(qMetric.', FNames{fid}, ');'])
+end
+qMetricArray = double(squeeze(reshape(table2array(struct2table(qMetric, 'AsArray', true)), size(qMetric.maxChannels, 2), ...
+    length(fieldnames(qMetric)))));
+qMetricTable = array2table(qMetricArray);
+qMetricTable.Properties.VariableNames = fieldnames(qMetric);
+
+parquetwrite([fullfile(savePath, 'templates._bc_qMetrics.parquet')], qMetricTable)
+
+% overwrite qMetric with the table, to be consistent with it for next steps
+% of the pipeline
+qMetric = qMetricTable;
+
+% optionally, also save output as a .tsv file
+if isfield(param,'saveAsTSV') % ensure back-compatibility if users have a previous version of param
+    fieldsToSave = {'percentageSpikesMissing_gaussian', ...
+        'presenceRatio', 'maxDriftEstimate', 'nPeaks', 'nTroughs', 'isSomatic','waveformDuration_peakTrough', 'spatialDecaySlope','waveformBaselineFlatness',...
+        'signalToNoiseRatio','fractionRPVs_estimatedTauR' };
+    fieldsRename = {'%_spikes_missing', 'presence_ratio', 'max_drift', 'n_peaks', 'n_troughs', ...
+        'is_somatic','waveform_dur', 'spatial_decay_slope','wv_baseline_flatness',...
+        'SNR','frac_RPVs' };
+    if param.saveAsTSV == 1 
+        cluster_id_vector = qMetric.clusterID - 1;
+        for fid = 1:length(fieldsToSave)
+            cluster_table = table(cluster_id_vector, qMetric.(fieldsToSave{fid}), 'VariableNames', {'cluster_id', fieldsRename{fid}});
+            writetable(cluster_table,[savePath 'cluster_' fieldsRename{fid} '.tsv'],'FileType', 'text','Delimiter','\t');  
+        end
+    end
+end
+
+end

qualityMetrics/bc_saveQMetrics.m

@@ -51,4 +51,26 @@
 % of the pipeline
 qMetric = qMetricTable;
 
+% optionally, also save output as a .tsv file
+if isfield(param,'saveAsTSV') % ensure back-compatibility if users have a previous version of param
+    fieldsToSave = {'percentageSpikesMissing_gaussian', ...
+        'presenceRatio', 'maxDriftEstimate', 'nPeaks', 'nTroughs', 'isSomatic','waveformDuration_peakTrough', 'spatialDecaySlope','waveformBaselineFlatness',...
+        'signalToNoiseRatio','fractionRPVs_estimatedTauR' };
+    fieldsRename = {'%_spikes_missing', 'presence_ratio', 'max_drift', 'n_peaks', 'n_troughs', ...
+        'is_somatic','waveform_dur', 'spatial_decay_slope','wv_baseline_flatness',...
+        'SNR','frac_RPVs' };
+    if param.saveAsTSV == 1 
+        cluster_id_vector = qMetric.clusterID - 1;
+        if isfield(param,'ephysKilosortPath')
+            saveTSV_path = param.ephysKilosortPath;
+        else
+            saveTSV_path = savePath;
+        end
+        for fid = 1:length(fieldsToSave)
+            cluster_table = table(cluster_id_vector, qMetric.(fieldsToSave{fid}), 'VariableNames', {'cluster_id', fieldsRename{fid}});
+            writetable(cluster_table,[saveTSV_path filesep 'cluster_' fieldsRename{fid} '.tsv'],'FileType', 'text','Delimiter','\t');  
+        end
+    end
+end
+
 end