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

@@ -30,18 +30,19 @@
 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);
 mef.loadMefAll;
 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

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

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 <https://www.gnu.org/licenses/>.
+%
+function groups = groupby(strs)
+
+    gs = cellfun(@(x) find(strcmp(x, strs)), unique(strs), 'UniformOutput', false);
+    groups = [unique(strs), gs];
+
+end

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