From 58cc651d24c94497f485d58b98a3054b2be84101 Mon Sep 17 00:00:00 2001 From: Harvey Huang Date: Wed, 14 Feb 2024 11:40:18 -0600 Subject: [PATCH] Minor bug fixes without impact to results, added groupby applyCARLARealCCEPs.m - loading mef from absolute path (using fullfile(pwd, dataPath)) because of known crash issue with readMef3 - annotPath set to correctly load from derivatives - now annots are cleared before loading next subject, and if no annots variable exists, rmBadTrialsAnnots is not called. This affects only subject 4, and did not affect results after fixing, because the corresponding trial annots previously preserved from subject 3 were all n/a. applyCARLARealCCEPsLoop - mef path set to be absolute - annotPath correctly points to derivatives - clear annots added before loading annots Added missing function groupby.m Minor comment updates to main --- applyCARLARealCCEPs.m | 18 +++++++++++------- applyCARLARealCCEPsLoop.m | 8 +++++--- functions/groupby.m | 32 ++++++++++++++++++++++++++++++++ main.m | 4 ++-- 4 files changed, 50 insertions(+), 12 deletions(-) create mode 100644 functions/groupby.m diff --git a/applyCARLARealCCEPs.m b/applyCARLARealCCEPs.m index e179bfb..cda2ec2 100644 --- a/applyCARLARealCCEPs.m +++ b/applyCARLARealCCEPs.m @@ -30,10 +30,11 @@ elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub)); elecs = ieeg_readtableRmHyphens(elecsPath); -mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); +% known issue: Matlab crashing when readMef3 loads mef data from relative path. Temporary fix = fullfile with pwd. If you specified an absolute "dataPath", remove "pwd". +mefPath = fullfile(pwd, dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run)); eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run)); -annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); +annotPath = fullfile(dataPath, 'derivatives', 'event_annotations', sub, sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); % Load CCEP data mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath); @@ -41,7 +42,7 @@ mef.highpass; mef.loadMefTrials([-1, 1]); -mef.plotOutputs(1, [], 200); +%mef.plotOutputs(1, [], 200); tt = mef.tt; srate = mef.srate; @@ -54,7 +55,8 @@ outdir = fullfile('output', 'realCCEPs'); mkdir(outdir); -try +clear annots +try % subject 4 has no annots because no interictal activity, hence the try annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t'); annots = annots.status_description; assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height'); @@ -67,7 +69,7 @@ disp('SOZs:'); disp(SOZ'); -clear *Path* +clear elecsPath mefPath channelsPath eventsPath annotPath % cleaning %% Extract data from site of interest, NaN out individual bad trials @@ -82,8 +84,10 @@ chsSite = channels.name(goodChs); % remove bad trials and set channels to nan in individual trials, according to annotations -annotsSite = annots(idxesSite); % get annotations corresponding to this stim site -dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite); +if exist('annots', 'var') + annotsSite = annots(idxesSite); % get annotations corresponding to this stim site + dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite); +end figure('Position', [200, 200, 1200, 1200]); subplot(1, 2, 1); % first half of channels diff --git a/applyCARLARealCCEPsLoop.m b/applyCARLARealCCEPsLoop.m index 0bc346a..4da0acd 100644 --- a/applyCARLARealCCEPsLoop.m +++ b/applyCARLARealCCEPsLoop.m @@ -44,10 +44,11 @@ run = runs{ss}(rr); - mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); + % known issue: Matlab crashing when readMef3 loads mef data from relative path. Temporary fix = fullfile with pwd. If you specified an absolute "dataPath", remove "pwd". + mefPath = fullfile(pwd, dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run)); channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run)); eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run)); - annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); + annotPath = fullfile(dataPath, 'derivatives', 'event_annotations', sub, sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run)); elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub)); % Load CCEP data from each subject @@ -65,7 +66,8 @@ sites = groupby(events.electrical_stimulation_site); % read manual annotations on events - try + clear annots + try % subject 4 has no annots because no interictal activity, hence the try annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t'); annots = annots.status_description; catch diff --git a/functions/groupby.m b/functions/groupby.m new file mode 100644 index 0000000..89ea915 --- /dev/null +++ b/functions/groupby.m @@ -0,0 +1,32 @@ +%% Creates a 2-column cell array where the first column indicates unique occurreences in the input, and second col tracks indices of each occurrence +% Modeled after Python groupby function +% +% 2024/02/14 +% +% If this code is used in a publication, please cite the manuscript: +% "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data" +% by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya, +% GA Worrell, KJ Miller, and D Hermes. +% +% CARLA manuscript package. +% Copyright (C) 2023 Harvey Huang +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . +% +function groups = groupby(strs) + + gs = cellfun(@(x) find(strcmp(x, strs)), unique(strs), 'UniformOutput', false); + groups = [unique(strs), gs]; + +end \ No newline at end of file diff --git a/main.m b/main.m index 345b886..c453e6a 100644 --- a/main.m +++ b/main.m @@ -37,7 +37,7 @@ % Make sure you are currently working in the root manuscript directory (CARLA) addpath('functions'); -% Data should be downloaded from OpenNeuro and copied to a folder here (./data) +% Data should be downloaded from OpenNeuro and copied to a folder here in the root directory (./data) dataPath = 'data'; % Add paths to other dependencies. Change the dummy paths listed here to where they are located on your system. @@ -140,7 +140,7 @@ %% Figure S3: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA % This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR -% Run after the previous section, keeping all workspace variables +% Run after the previous section, keeping all workspace variables. Output also saved to .\output\realCCEPs\. figS3_compareCARvCARLA %% Figure S4: Apply CARLA to real data at another 4 stim sites in subjects 1-4