loadPreprocessedMef.m

%% Load preprocessed data from across runs for subject, after running preprocessFromMef
% This script loads the output from preprocessFromMef, sets SOZ and stimulated neighbor channels to bad, saves an output ERP heatmap for all channels.
% Then use either analysisERP or analysisBB to do further analysis
%
%    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/>.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

set(0, 'DefaultFigureRenderer', 'painters');

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

localPath = 'data';
dataPath = fullfile(localPath, sprintf('sub-%s', sub), 'ses-ieeg01', 'ieeg');

switch upper(sub)
    case '2'
        runs = 1:3;
        imgType = 'Coh';
    case '1'
        runs = 1:4;
        imgType = 'Coh';
    otherwise
        error('Error: choose subjects 1 or 2');
end

temp = struct();
M = [];
Meye = [];
events = [];
for ii = 1:length(runs)
    fprintf('Loading run %d... ', runs(ii));
    temp = load(fullfile('output', sprintf('sub-%s', sub), 'ccepVisualPreproc', sprintf('ccepVisualPreprocRun%02d', runs(ii))));
    
    M = cat(3, M, temp.Mcar);
    Meye = cat(3, Meye, temp.Meye);
    
    channels = ieeg_readtableRmHyphens(fullfile(dataPath, sprintf('sub-%s_ses-ieeg01_task-ccepVisual_run-%02d_channels.tsv', sub, ii)));
    ephysInds = find(strcmp(channels.type, 'SEEG')); assert(all(diff(ephysInds) == 1), 'Error: expecting contiguous ephys channels');
    channelsEphysCurr = channels(ephysInds, :); % ephys channels for current run
    
    eventsCurr = ieeg_readtableRmHyphens(fullfile('data', 'derivatives', 'events_annotations', sprintf('sub-%s', sub), sprintf('sub-%s_eventsCcepVisualRun%02d.tsv', sub, runs(ii))), 'electrical_stimulation_site', 1);
    eventsCurr.run = repmat(ii, height(eventsCurr), 1);
    
    events = [events; eventsCurr];
    
    if ii == 1
        srate = temp.srate;
        tt = temp.tt;
        channelsEphys = channelsEphysCurr;
    else % check that variables are the same for other runs
        assert(all(temp.tt == tt), 'tt is different for run %d', runs(ii));
        assert(temp.srate == srate, 'Srate is different for run %d', runs(ii));
        channelsEphys.status(strcmp(channelsEphysCurr.status, 'bad')) = {'bad'}; % set bad channels to be those bad at all runs
    end
    
    fprintf('%d trials\n', size(temp.Mcar, 3));
end

% 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'));
channelsEphys.status(ismember(channelsEphys.name, elecsSOZ)) = {'bad'};

% set electrodes neighboring stim sites to bad
stimSites = unique(events(~isnan(events.e2v), :).electrical_stimulation_site);
for ii = 1:length(stimSites)
    stimChs = getNeighborChs(split(stimSites{ii}, '-'), 2);
    channelsEphys.status(ismember(channelsEphys.name, stimChs)) = {'bad'};
end

clear temp channelsEphysCurr eventsCurr channels stimChs

%% Load preprocessed CCEPs for each stim site

dataCCEP = struct();
for ii = 1:size(stimSites, 1)
    fprintf('Loading CCEPs for stim site %s\n', stimSites{ii});
    temp = load(fullfile('output', sprintf('sub-%s', sub), sprintf('CCEP_%s_preproc.mat', stimSites{ii})));
    dataCCEP(ii).site = stimSites{ii};
    dataCCEP(ii).MCcep = temp.MCcep;
    dataCCEP(ii).chsUsed = temp.chsUsed{1};
    dataCCEP(ii).ttCcep = temp.ttCcep;
end
ttCcep = dataCCEP(1).ttCcep;

clear temp

%% Summary of E2V intervals, get logical trial indices corresponding to each image condition

figure; histogram(events.e2v/srate, 20); xlabel('Elec -> Visual Time (s)'); ylabel('Frequency'); title('All E2Vs');
fprintf('%d trials with sham stim\n', sum(isnan(events.e2v)));

% get indices with a tolerance around the desired E2V conditions
e2vs = struct();
e2vProg = [nan, 0, 0.1, 0.2]; % Programmed E2V durations, in s
e2vs(1).mode = nan; e2vs(1).idx = isnan(events.e2v); % the sham condition
for ii = 2:4
    [e2vs(ii).mode, e2vs(ii).idx] = getE2vs(e2vProg(ii), events.e2v/srate, 0.05);
end

cmE2v = parula(5); cmE2v = [0, 0, 0; cmE2v(2:4, :)]; % colormap
e2vCaptions = arrayfun(@(x) sprintf('%0.0fms', 1e3*x), [e2vs.mode], 'UniformOutput', false); % string for each e2v

% Plot evoked potential map
figure('Position', [200, 200, 600, 1200]);
evokedMap = mean(M(:, :, contains(events.trial_type, sprintf('high%s', imgType)) & isnan(events.e2v)), 3); % evoked map of strongest visual responses (100% coherence) but no stim
%evokedMap = mean(M, 3); % evoked map of all trials
evokedMap(strcmp(channelsEphys.status, 'bad'), :) = nan;
imagescNaN(tt, 1:height(channelsEphys), evokedMap, 'CLim', [-50, 50]);
xline(0, 'Color', 'r', 'LineWidth', 1);
set(gca, 'YTick', 1:2:height(channelsEphys), 'YTickLabels', channelsEphys.name(1:2:end));
xlim([-0.5, 1]); xlabel('time (s)'); ylabel('channels'); title('Evoked potential heatmap');
colorbar;

% image categories, in same format as the e2v struct
imgs = struct();
imglabs = cellfun(@(x) sprintf('%s%s', x, imgType), ...
    {'zero', 'img1-low', 'img1-med', 'img1-high', 'img2-low', 'img2-med', 'img2-high'}, 'UniformOutput', false);
for ii = 1:7
    imgs(ii).label = imglabs{ii};
    imgs(ii).idx = contains(events.trial_type, imglabs{ii});
end
clear imglabs;

%% Bipolar rereference data for broadband analysis (used by analysisBBFIR1 when bipolar referencing specifically examined)
% Note that data has been previously CAR-rereferenced and that we aren't considering lead segment info here

badChans = find(strcmpi(channelsEphys.status, 'bad'));

% rereference first trial to get number of bipolar channels, names, and bad channels
[tempBip, bipolarNames, bipolarChans, badChansBip] = ieeg_bipolarSEEG(M(:, :, 1)', channelsEphys.name, badChans, [], [], false);

% rereference all trials
Mbip = nan(size(tempBip, 2), size(M, 2), size(M, 3));
fprintf('Bipolar re-referencing');
for ii = 1:height(events)
    temp = ieeg_bipolarSEEG(M(:, :, ii)', channelsEphys.name, badChans, [], [], false);
    Mbip(:, :, ii) = temp';
    showdot;
end
fprintf('\n')

% same bipolar processing for CCEP. Channels should be the same as bipolar
for ii = 1:2
    dataCCEP(ii).Mbip = nan(size(tempBip, 2), size(dataCCEP(ii).MCcep, 2), size(dataCCEP(ii).MCcep, 3));
    for jj = 1:size(dataCCEP(ii).MCcep, 3)
        temp = ieeg_bipolarSEEG(dataCCEP(ii).MCcep(:, :, jj)', channelsEphys.name, badChans, [], [], false);
        dataCCEP(ii).Mbip(:, :, jj) = temp';
        showdot;
    end
end
fprintf('\n');

return;