eegmp_preproc/code/a7_prep_data_for_analysis/a7_prep_data_for_analysis.m

function a7_prep_data_for_analysis(id, rootpath)

if ismac % run if function is not pre-compiled
    id = '11'; % test for example subject
    currentFile = mfilename('fullpath');
    [pathstr,~,~] = fileparts(currentFile);
    cd(fullfile(pathstr,'..', '..'))
    rootpath = pwd;
end

% inputs
pn.eeg_BIDS = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'eeg_BIDS');
pn.channel_locations = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'channel_locations');
pn.events = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'events');
pn.tools = fullfile(rootpath, 'tools');
% outputs
pn.eeg_ft = fullfile(rootpath, 'data', 'outputs', 'eeg');
pn.history = fullfile(rootpath, 'data', 'outputs', 'history');

if ismac % run if function is not pre-compiled
     addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
     addpath(fullfile(pn.tools, 'helpers'));
end

% N = 33
IDs = tdfread(fullfile(pn.eeg_BIDS, 'participants.tsv'));
IDs = cellstr(IDs.participant_id);

id = str2num(id);

display(['processing ID ' num2str(IDs{id})]);

%% load files

% load data
load(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_seg.mat']),'data_eeg', 'events');

data_ica = load(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_ica.mat']),'data');
data = data_eeg; clear data_eeg;
data.elec = data_ica.data.elec;
data.chanlocs = data_ica.data.chanlocs; clear data_ica;

% load config
load(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

% cache original number of trials
Ntrial = numel(data.trial);

%% define parameter

% ICs to reject
iclabels    = config.ica1.iclabels.manual;
ics2reject = [iclabels.bli(:); iclabels.mov(:); iclabels.hrt(:); iclabels.ref(:); iclabels.art(:); iclabels.emg(:)];
clear iclabels

% trials to remove
trls2remove = config.ArtDect.trials;

% channels to interpolate
chns2interp = config.ArtDect.channels;

%% ICA (from weights)

% ica config
cfg.method           = 'runica';
cfg.channel          = {'all'};
cfg.trials           = 'all';
cfg.numcomponent     = 'all';
cfg.demean           = 'yes';
cfg.runica.extended  = 1;

% ICA solution
cfg.unmixing     = config.ica1.unmixing;
cfg.topolabel    = config.ica1.topolabel;

% components
comp = ft_componentanalysis(cfg,data);

% clear cfg
clear cfg data

%% remove components

% cfg for rejecting components (reject: blinks, eye movements, ecg, ref)
cfg.component = sortrows(ics2reject)';
cfg.demean    = 'yes';

% reject components
data = ft_rejectcomponent(cfg,comp);

% clear cfg
clear cfg comp ics2reject

%% remove trials

% define trials to keep
trials               = 1:Ntrial;
trials(trls2remove)  = [];

% config for deleting trials
cfg.trials   = trials;

% remove trials
data   = ft_preprocessing(cfg,data);
nartfree_trials = trials;

config.nartfree_trials = nartfree_trials;

% clear variables
clear cfg trials trls2remove

%% remove eye & reference channels (before interpolation)

cfg.channel     = {'all','-IO1','-IO2','-M1', '-M2','-VEOG','-HEOG'};
cfg.demean      = 'yes';

% remove channels
data = ft_preprocessing(cfg,data);

% clear cfg
clear cfg

if ~isempty(chns2interp)
    % interpolation cfg
    cfg            = [];
    cfg.method     = 'spline';
    cfg.badchannel = data.label(chns2interp);
    cfg.trials     = 'all';
    cfg.lambda     = 1e-5; 
    cfg.order      = 4; 
    cfg.elec       = config.elec;

    % interpolate
    data = ft_channelrepair(cfg,data);

    % clear cfg
    clear cfg chns2interp
end

%% keep channel x trial artifacts  

data.ArtDect.channels_x_trials = config.ArtDect.channels_x_trials;

%% create overview of percentage of excluded trials, update config.trl

tmp_trialNumbers(1,1) = length(data.trial);
tmp_trialNumbers(2,1) = 1200; % maximum amount of trials
tmp_trialNumbers(3,1) = (100-((tmp_trialNumbers(1,1)/tmp_trialNumbers(2,1))*100)); % percentage of excluded trials
config.overview_trial_numbers = tmp_trialNumbers;

%% update config.trl with trigger codes

% add info whether trial is included in final set; 1 = trial containing no artifact
config.trl(1:Ntrial,4)                 = 0;
config.trl(nartfree_trials, 4)  = 1;
% add trial number
config.trl(1:Ntrial,5) = 0;
config.trl(nartfree_trials, 5)  = nartfree_trials;
% load marker file
mrk = config.mrk;
for i = 1:size(mrk,2)
   mrk_val{i,1} = mrk(1,i).value;
end

% save config
save(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

%% update event table

% remove trials that were not retained during preprocessing

events(config.trl(:,5)==0,:) = [];

%% save data
data_eeg = data; clear data;
save(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_art.mat']),'data_eeg', 'events');