preprocessCCEPForCCEPVisual.m

%% This script specifically preprocesses stim sites used for CCEPvisual and preprocesses them analogously to the ccepVisual data, then saving as output
%
% 1) Load CCEP data and highpass using ccep_PreprocessMef
% 2) Isolate trials corresponding to stim site of interest only
% 3) Identify SOZ electrodes
% 4) Nan-out interictal trials, stim-neighbor electrodes, bad channels, and perform CAR using bottom 25% by variance
%
%    If this code is used in a publication, please cite the manuscript:
%    "H Huang, KN Kay, NM Gregg, G Ojeda Valencia, M In, C Kapeller, Y Shu, GA Worrell, KJ Miller, and D Hermes.
%    Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex. (Under Review)"
%
%    A preprint is available currently at doi: https://doi.org/10.1101/2025.05.05.652264.
%
%    The dataset corresponding to this code and manuscript is in BIDS format (version 1.10.0) on OpenNeuro (ds006485),
%    and it will be made publicly available upon manuscript acceptance.
%
%    Copyright (C) 2025 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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% Configure paths

sub = upper(input('Subject? ', "s"));

dataPath = fullfile('data', sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg');
elecsPath = fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_electrodes.tsv', sub));

% determines which CCEP runs (not ccepvisual)
chReref = '';
switch upper(sub)

    case '2'
        runs = [1, 2]; % CCEP runs matching stimSites
        stimSites = {'ROC6-ROC7', 'RMO3-RMO4'};
    case '1'
        runs = [1, 1]; % CCEP runs matching stimSites
        stimSites = {'LOC4-LOC5', 'LG6-LG7'};
        chReref = 'LY5'; % perform monopolar rereferencing before CAR to match ccepVisual
    otherwise
        error('Error: choose subjects 1 or 2');
end

%% Iterate across runs

for ii = 1:length(runs)

    %%

    run = runs(ii);
    stimSite = stimSites{ii};
    
    mefPath = fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
    channelsPath = fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
    eventsPath = fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));

    % read events to remove bad events from annots
    eventsTemp = readtable(eventsPath, 'FileType', 'text', 'Delimiter', '\t');
    eventsBadIdx = ~strcmpi(eventsTemp.status, 'good');

    % Load CCEP data and highpass
    mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath);
    mef.loadMefAll;
    mef.highpass;
    mef.loadMefTrials([-1, 1.5]);
        
    tt = mef.tt;
    srate = mef.srate;
    channels = mef.channels;
    events = mef.evts;
    data = mef.data;
    
    % load and identify SOZ channels
    elecs = ieeg_readtableRmHyphens(elecsPath);
    elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
    disp('SOZs:'); disp(elecsSOZ');

    % Identify stim-neighboring channels
    stimChs = getNeighborChs(split(stimSite, '-'), 2);

    % load annots and process data according to annots (manually done by mnl_DataCuration/mnl_ccep_mefEpochReview)
    annotPath = fullfile('data',  'derivatives', 'event_annotations', sprintf('sub-%s', sub), sprintf('sub-%s_ses-ieeg01_task-ccep_run-%02d_events_trialStatus.tsv', sub, runs(ii)));
    annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
    annots(eventsBadIdx, :) = []; % remove rows that were deleted when ccep_PreprocessMef loaded events.
    assert(height(annots) == height(events), 'Error: mismatch between number of trials in events and number of rows in annots');
    data = rmBadTrialsAnnots(data, channels.name, annots.status_description);
    events(strcmp(annots.status_description, 'all'), :) = [];
    annots(strcmp(annots.status_description, 'all'), :) = [];

    clear eventsTemp eventsBadIdx

    % Keep only trials corresponding to stim site of interest
    trsStimSite = strcmp(events.electrical_stimulation_site, stimSite);
    events = events(trsStimSite, :);
    data = data(:, :, trsStimSite);
    annots = annots(trsStimSite, :);
    
    % monopolar rereferencing so CAR is more reliable, consistent with ccepVisual preprocessing
    if ~isempty(chReref)
        data = data - data(strcmp(channels.name, chReref), :, :);
        channels.status{strcmp(channels.name, chReref)} = 'bad'; % set the 0-ed out channel to bad so it isn't used in CAR
    end

    % CAR by variance
    grp = ones(height(events), 1); % group all trials together (already one stim site)
    badChs = find(strcmp(channels.status, 'bad') | ismember(channels.name, stimChs) | ismember(channels.name, elecsSOZ) | ~strcmp(channels.type, 'SEEG'));
    opts = struct();
    opts.vartype = 'var'; % use geometric mean variance as ranking criteria to be consistent with ccepVisual trials
    [Mcar, chsUsed] = CARVariance(tt, data, srate, badChs, grp, opts);

    % Save output matrix for use with ccepVisual
    ephysInds = find(strcmp(channels.type, 'SEEG')); assert(all(diff(ephysInds) == 1), 'Error: expecting contiguous ephys channels');
    MCcep = Mcar(ephysInds, :, :);
    ttCcep = tt;

    save(fullfile('output', sprintf('sub-%s', sub), sprintf('CCEP_%s_preproc.mat', stimSite)), 'ttCcep', 'MCcep', 'chsUsed');
    fprintf('\n');

end