preprocessFromMef.m

%% This preprocessing script is adapated from the preprocess.m, and works on Mef-converted data
%
%    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 subjects and paths, 

addpath('functions');

sub = upper(input('Subject? ', "s")); % 1 or 2

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

chReref = ''; % whether to perform rereference before CAR (if original reference was inconsistent). Use an unimportant WM ch
switch upper(sub)
    case '2'
        runs = 1:3;
    case '1'
        runs = 1:4; % LOC4-5 for runs 1-2, LG6-7 for runs 3-4
        chReref = 'LY5';
    otherwise
        error('Error: enter subject 1 or 2');
end

srate = 4800; % hard coded cuz this never changes

mkdir('output', sprintf('sub-%s', sub));

%% Save ch x time x trial data across runs for one subject

for rr = 1:length(runs)
    
    %% Load channels and events. Convert to center on visual events with preceding rest and stim ISI
    
    channels = ieeg_readtableRmHyphens(fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccepVisual_run-%02d_channels.tsv', sub, rr)));
    events = ieeg_readtableRmHyphens(fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccepVisual_run-%02d_events.tsv', sub, rr)), 'electrical_stimulation_site', 1);
    
    % load electrodes and set SOZ channels to bad
    elecs = ieeg_readtableRmHyphens(fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_electrodes.tsv', sub)));
    elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
    channels.status(ismember(channels.name, elecsSOZ)) = {'bad'};
    fprintf('Setting %d SOZ channels to bad in channels.tsv\n', length(elecsSOZ));
    
    % display minimum rest and image durations -- should be ~1s
    minRestDur = min(events.duration(contains(events.trial_type, 'rest')));
    minImgDur = min(events.duration(contains(events.trial_type, {'img', 'zeroContrast'}))); % contrast type experiments did not have 'img' in name
    fprintf('Minimum rest dur: %0.03fs, image dur: %0.03fs\n', minRestDur, minImgDur);
    
    % separate visual and electrical events
    eventsVisual = events(contains(events.trial_type, {'img', 'zeroContrast'}), :);
    eventsElec = events(strcmp(events.trial_type, 'electrical_stimulation'), :); 
    
    % find e.stim corresponding to visual event, if it exists. Contatenate event annotations
    e2vs = nan(height(eventsVisual), 1); % isi from electrical -> visual
    for ii = 1:height(eventsVisual)
        e2vsCurr = eventsVisual.sample_start(ii) - eventsElec.sample_start;
        ind = find(e2vsCurr < 0.3*srate & e2vsCurr > -0.1*srate); % e2v ISI should be on range of [0, 200 ms]
        if ~isempty(ind)
            e2vs(ii) = e2vsCurr(ind); % found a matching stim trial
            eventsVisual(ii, 5:8) = eventsElec(ind, 5:8);
        end
    end
    eventsVisual.e2v = e2vs;
    
    % modes for the 200 ms, 100 ms, and 0 ms ISIs
    modes = [mode(e2vs(e2vs > 0.15*srate)), mode(e2vs(e2vs > 0.05*srate & e2vs < 0.15*srate)), mode(e2vs(e2vs < 0.05*srate))];
    fprintf('E2V modes at %0.03fs, %0.03fs, %0.03fs\n', modes(1)/srate, modes(2)/srate, modes(3)/srate);
    
    figure; histogram(eventsVisual.e2v/srate, 30); title('E2V distributions');
    
    % write prelim events based on each visual trial so that trials can be annotated in BIDS
    writetable(eventsVisual, fullfile('output', sprintf('sub-%s', sub), sprintf('eventsCcepVisualRun%02dPRELIM.tsv', rr)), 'FileType', 'text', 'Delimiter', '\t');
    
    %% Load all mef data and high pass
    
    mefPath = fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccepVisual_run-%02d_ieeg.mefd', sub, rr));
    [metadata, data] = readMef3(mefPath);
    
    % highpass SEEG channels
    data(strcmp(channels.type, 'SEEG'), :) = ieeg_highpass(data(strcmp(channels.type, 'SEEG'), :)', srate)';
    
     %% Calculate button presses for each trial. Done before removing trials
    
    buttons = data(ismember(channels.name, {'DigitalInput3', 'DigitalInput4'}), :);
    rightOn = find(diff(buttons(1, :) == 1)) + 1;
    leftOn = find(diff(buttons(2, :) == 1)) + 1;
    
    eventsVisual.respTime = nan(height(eventsVisual), 1);
    eventsVisual.respButton = nan(height(eventsVisual), 1);
    for ii = 1:height(eventsVisual)
        
        % Right button -- only keep ones that come after stim
        diffs = rightOn - eventsVisual.sample_start(ii); diffs(diffs <= 0) = inf;
        respCurrRight = min(diffs);

        % Left button
        diffs = leftOn - eventsVisual.sample_start(ii); diffs(diffs <= 0) = inf;
        respCurrLeft = min(diffs);

        try % time to next trial
            trialDur = eventsVisual.sample_start(ii + 1) - eventsVisual.sample_start(ii);
        catch
            trialDur = inf;
        end

        if isinf(respCurrRight) && isinf(respCurrLeft), continue; end % no more responses

        if respCurrRight < respCurrLeft
            respMin = respCurrRight; respCode = 3; % digital input 3
        elseif respCurrLeft < respCurrRight
            respMin = respCurrLeft; respCode = 4; % digital input 4
        else
            blah; % both response times are equal
        end

        % add response time and choice only if responses were before next sample start
        if respMin < trialDur
            eventsVisual.eventsrespTime(ii) = respMin/srate;
            eventsVisual.respButton(ii) = respCode;
        end
    end
    
    %% At this point we used mnl_DataCuration/mnl_mefEpochReview to annotate good/bad trials
    % The annotated trials are in derivatives/event_annotations/. See "status" and "status_description" columns

    %% Load visual trial annotations and remove bad or 'all' interictal trials
    
    try
        eventsAnnot = readtable(fullfile('output', sprintf('sub-%s', sub), sprintf('eventsCcepVisualRun%02dPRELIM_trialStatus.tsv', rr)), ...
                                    'FileType', 'text', 'Delimiter', '\t');
    catch
        warning('No event annots available in output for run %d -- generate annots', rr);
        continue
    end
    
    % separately load and use annotations from curated eventsAnnot.tsv files
    eventsVisual.status = eventsAnnot.status;
    eventsVisual.status_description = eventsAnnot.status_description;
    
    % Remove bad trials and those that are annotated 'all' for interictal
    trsRm = strcmp(eventsVisual.status, 'bad') | strcmp(eventsVisual.status_description, 'all');
    fprintf('Removing %d bad or interictal visual trials\n', sum(trsRm));
    eventsVisual(trsRm, :) = [];

    %% Convert to trial structure
    
    ephysInds = find(strcmp(channels.type, 'SEEG')); assert(all(diff(ephysInds) == 1), 'Error: expecting contiguous ephys channels');
    
    eyeInds = find(contains(channels.name, 'Eyetracker'));
    
    % sample range and time points for each trial
    samprange = -srate:2*srate-1;
    tt = samprange/srate;
    
    M = zeros(length(ephysInds), length(samprange), height(eventsVisual)); % the ephys data only
    Meye = zeros(length(eyeInds), length(samprange), height(eventsVisual)); % the eyetracker data only
    for ii = 1:height(eventsVisual)
        M(:, :, ii) = data(ephysInds, (eventsVisual.sample_start(ii)+samprange(1)+1) : (eventsVisual.sample_start(ii)+samprange(end)+1)); % convert to 1 indexing
        Meye(:, :, ii) = data(eyeInds, (eventsVisual.sample_start(ii)+samprange(1)+1) : (eventsVisual.sample_start(ii)+samprange(end)+1));
    end
    
    channelsEphys = channels(ephysInds, :);
    eyetrackerLabels = channels.name(eyeInds);
    
    % Plot first 20 good channels, mean trials to look for possible large stim artifact across channels. If present, should perform virtual rereferencing before CAR
    chsTest = find(strcmp(channelsEphys.status, 'good'), 20, 'first');
    figure('Position', [400, 400, 1000, 600]);
    plot(tt, mean(M(chsTest, :, :), 3)'); xlim([-0.5, 1]); ylim([-1000, 1000]);
    title('Mean trial across first 20 good channels'); xlabel('Time'); ylabel('\muV');
    
    %% Common average rereference
    
    % Perform a "virtual" monopolar rereferencing to remove stim artifacts so that CAR variance calcuation is more reliable
    % Do this if the plot above indicated large stim artifacts across multiple channels
    if ~isempty(chReref)
        M = M - M(strcmp(channelsEphys.name, chReref), :, :);
        channelsEphys.status{strcmp(channelsEphys.name, chReref)} = 'bad'; % set the 0-ed out channel to bad so it isn't used in CAR
    end
    
%     % ALTERNATIVE WAYS TO GROUP FOR CAR
%     % Indices corresponding to different E2V categories, to parse group for adjusted CAR, (if grouping by isi is desired). 1 = no stim
%     e2vTypes = ones(height(eventsVisual), 1);
%     e2vTypes(e2vs < 0.05*srate) = 2; % simultaneous
%     e2vTypes(e2vs > 0.05*srate & e2vs < 0.15*srate) = 3; % short e2v
%     e2vTypes(e2vs > 0.15*srate) = 4; % long
%     
%     % convert trial_type strings (image description) to numerical indices
%     grp = cellfun(@(x) sum(double(x)), eventsVisual.trial_type);
%     grp = grp.*e2vTypes; % create trial_type x stim conditions (3 for no stim, 6 for stim conditions)

%     % group trials for re-referencing only by image. This means the same channels are used in common average across coherence & e2v of the same image
%     % Why not just one group for both images? Rationale: don't care as much about directly comparing img1 to img2. Possibly img-dependent activation map in brain (e.g. ffa, ppa)
%     grp = zeros(height(eventsVisual), 1);
%     grp(contains(eventsVisual.trial_type, 'img1')) = 1;
%     grp(contains(eventsVisual.trial_type, 'img2')) = 2;
%     assert(~any(grp == 0), 'Error: Unassigned group identity for trial');
    
    % group all trials together for rereferencing, so that 0-coh conditions can be pooled between both images
    grp = ones(height(eventsVisual), 1);
    
    stimPair = unique(eventsElec.electrical_stimulation_site);
    assert(length(stimPair) == 1, 'Error: Expecting 1 stim site per run');
    stimChs = getNeighborChs(split(stimPair, '-'), 2);
    
    % identify SOZ channels
    elecsSOZ = elecs.name(contains(elecs.seizure_zone, 'SOZ'));
    
    % Nan-out bad trials at individual channels
    for tr = 1:height(eventsVisual)
        chsExclCurr = unique(upper(strip(split(eventsVisual.status_description(tr), ',')))); % some cleaning, make sure no duplicates
        chsExclCurr(cellfun(@isempty, chsExclCurr)) = []; % remove accidental channels after/before commas
        if strcmp(chsExclCurr, 'N/A'), continue; end
        
        lia = ismember(channelsEphys.name, chsExclCurr);
        assert(sum(lia) == length(chsExclCurr), 'Error: Mismatch between channels to be excluded and actual channels in channels table, at trial %d', tr);
        M(lia, :, tr) = nan; % bad channels set to nan for current trial. These will be ignored when computing CAR
    end
    
    % CAR, with bad channels being either bad in the channels table (including the SOZ set earlier) or is a stim neighbor, or is a SOZ
    badChs = find(strcmp(channelsEphys.status, 'bad') | ismember(channelsEphys.name, stimChs) | ismember(channelsEphys.name, elecsSOZ));
        
    opts = struct();
    opts.vartype = 'var'; % use geometric mean variance as ranking criteria because the trials don't necessarily have same shape (different e2v, coherence) 
    [Mcar, chsUsed] = CARVariance(tt, M, srate, badChs, grp, opts);
    
    %% Save trial data, final events table
    
    save(fullfile('output', sprintf('sub-%s', sub), sprintf('ccepVisualPreprocRun%02d.mat', rr)), 'tt', 'srate', 'Mcar', 'Meye', 'eyetrackerLabels', 'chsUsed', '-v7.3');
    writetable(eventsVisual, fullfile('output', sprintf('sub-%s', sub), sprintf('sub-%s_eventsCcepVisualRun%02d.tsv', sub, rr)), 'FileType', 'text', 'Delimiter', '\t');
    
    % Plot and save map of EPs
    clim = [-50, 50];
    evoked = mean(Mcar, 3, 'omitnan');
    evoked(strcmp(channelsEphys.status, 'bad'), :) = nan;
    figure('Position', [200, 200, 600, 1200]);
    imagescNaN(tt, 1:height(channelsEphys), evoked, 'CLim', clim);
    for ii = 1:length(modes)
        xline(-modes(ii)/srate, 'Color', 'r', 'LineWidth', 1); % mark the stim modes
    end
    xlim([-0.5, 1]);
    set(gca, 'YTick', 1:2:height(channels), 'YTickLabels', channels.name(1:2:end));
    title('Evoked potential heatmap');
    ylabel('Channels');
    xlabel('Time (s)');
    colorbar;
    
    saveas(gcf, sprintf('output/sub-%s/epCcepVisualRun%d.png', sub, rr));
    
end