Guidali_et_al_2023_EJN_RR/Script preprocessing EMG-EEG/preprocessingTEP_pipeline.m

 %%% Pipeline case gtec files

clear all
close all

%main_dir = fileparts(mfilename('fullpath')); % path of the script folder

%% paths for PLUG-INs functions

addpath([main_dir '\PLUG-INs\SOUND']);%SOUND
addpath([main_dir '\PLUG-INs\TESA1.1.1']); %TESA
addpath([main_dir '\PLUG-INs\M_functions']); % specific functions for the script
%path of the raw data
%addpath('/raw_data');



% initialize eeglab;
%addpath ('eeglab2022.1/') %% add path eeglab
%eeglab; close; %initialize eeg lab

%initialize fieldtrip
ft_defaults

%% set folders for input and output: Y
%%%% folder with data to be processed
% datafolder_parent = [main_dir '\raw_data'];  % parent data folder includes all the subdir of all dataset to be processed
% mkdir(main_dir, '\analysis_MOBI'); % build output folder
% output = [main_dir, '\analysis_MOBI'];

%% Variables

%%control to overwrite steps
global_overwrite=0 ;% 1 to overwrite all steps
prep1_1_overwrite=0; % 1 to overwrite step 1.1
prep1_2_overwrite=0; % 1 to overwrite step 1.2
sound_overwrite=0; %1 to overwrite step 2 SOUND
rej_overwrite=0; % 1 to overwrite step 3 artifact rejection
ICA_overwrite=0; % 1 to overwrite step 4.1 ICA
ICA_remove=0; %  1 to overwrite step 4.2 ICA selection components
SSP_SIR_step1 =0; % 1 to overwrite step 5.1 SSP-SIR
SSP_SIR_step2 =0; % 1 to overwrite step 5.2 SSP-SIR components selection

%%step 1.1
chan_loc_st = 'standard-10-5-cap385.elp';% standard 10-20 channel location
load([main_dir '\PLUG-INs\M_functions\EEG_chanlocs_73' ]); chan_loc=EEG_chanlocs_73;%custom channel location
ref_chan = 2;%row of the reference channel
emg_refchan= [76, 78];
emg_chan = [77, 79];
eeg_channels=73;
tot_cond=7;

epoching_long=[-0.7 0.7];% interval to epoch data, latency in sec relative to time-locking event
rmBase = [-200 -5]; % interval for baseline correction

%%step1.2
interp = [0 3]; %interval to cut and interpolate TMS pulse for TESA method, time in ms related to the event
interp_Win = [0.35,0.35];%times before and after artifact window for fitting cubic function(TESA) 
HP_f = 1; % % HP filter frequency
fftflag=1; % filterinf in the frequency domain, faster for high-order filters
DS= 4800; %down sampling frequency


%%step2 SOUND
ref_before_sound= 'FPz';
lambda_value = 0.1; % SOUND parameter
iter = 5;%SOUND parameter

%%step 3 artifact rejection
SD_first=5; %number of stnd dev to consider when computing joint probability and remove trial (first trial rejection)
SD_last=5; %number of stnd dev to consider when computing joint probability and remove trial(last trial rejection)

%%step 4.2 ICA
threshold = [0 0;0 0; 0.5 1; 0 0; 0 0; 0 0; 0 0];% threshold to label eye components (50%)

%%step 5 SSP-SIR
SSP_SIR_window = [0,50];%the time window on which SSP-SIR is applied 
limx=[-100 350];% xlim for figure with TEPs
limy = [-30 30]; %ylim for figure with TEPs

%%step 6 final preprocesing
LP_f = 70;% frequency in Hz for lowpass filter
%chan_loc2 = 'standard-10-5-cap385.elp';% standard 10-20 channel location
old_ref = 'FPz';
reref = {'TP9' 'TP10'};
epoching_short = [-0.2 0.4];
latency = [-0.1 0.4]; %epoch limit in final fieldtrip set
ylimit=[-15 15];%ylimit for final figure


%% 0. Get file ready for processing

% % function to build the mat file row data
[subjects,n_subjects]= get_data_MOBI(datafolder_parent,output,ref_chan,emg_refchan,tot_cond);
save('subjects.mat', 'subjects')
save('n_subjects.mat', 'n_subjects')

%% 1. BASIC PREPROCESSING
% Step 1.1
% load set, add channel location, remove unwanted electrodes, define
% reference, save the set

%uncomment these two lines if n_sibjects.mat already exists (and skip/comment
%section 0.):
%load n_subjects.mat 
%load subjects.mat 

for iSub = 1:n_subjects
    subj =  subjects{iSub}
    load([output '\' subj '\0_Data\Block_cond.mat'])
    
    for n_cond= 1:tot_cond
               
        name_condition = [subj '_' Block_cond{n_cond,2}];
        name_folder= [subj '\' name_condition];
        
        mkdir([output, '\' name_folder '\1_Preprocessing\step1_1']);
        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_1'];
    
        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_1.set']) ...
                && isfile([current_output_folder '\' Block_cond{n_cond,1} '_EMG_step1_1.set'])...
                && global_overwrite==0 && prep1_1_overwrite==0
                warning(['skipping step1_1 for ' name_condition ' , file already present'])
        else
            %load row data
            load([output '\' subj '\0_Data\' subj '_' Block_cond{n_cond,2}]);
            data = eval([Block_cond{n_cond,2}]);          

            % load raw data,discard EMG channels,read the trigger,load channel names
            [EEG, EMG]= load_rawdata_MOBI(data,name_folder,Block_cond{n_cond,3},Block_cond{n_cond,5},Block_cond{n_cond,4},eeg_channels,chan_loc);

            %plot before starting preprocessing signals. Epoching and baseline correction only
            %for visual purpose
            myplot_F(EEG,[-100 500],[-1000 45000],['Step 1.1 signals'],epoching_long);
            saveas(gcf, [current_output_folder '\step1_1'] )
            close
            
            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1} '_step1_1'],'filepath', current_output_folder);
            EMG = pop_saveset(EMG, 'filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'], 'filepath',current_output_folder );

        end

        disp('Step 1.1 ended')
   

%% STEP 1.2 Remove TMS-pulse artifact and interpolate, HP filter, dwnsampling,epoching 
        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_1'];
        mkdir([output, '\' name_folder '\1_Preprocessing\step1_2']);
        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_2'];

        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2.set']) && global_overwrite==0 && prep1_2_overwrite ==0
                warning(['skipping step1_2 for ' name_condition ' , file already present'])       
        else
            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_1.set'],'filepath', previous_folder);

            % function to remove TMS pulse with TESA Method
            % inputs: EEG, interpolation interval,times before and after artifact window for fitting cubic function 
           EEG= remove_TMS_TESA_G(EEG,interp,interp_Win);
           %plot with epoching and baseline correction only for visual purpose
           myplot_F(EEG,[-100 300],[-300 300],['After interpolation TESA'],epoching_long,rmBase);
           saveas(gcf, [current_output_folder '\after_interpolation_TESA'])
           close

            %%%%% High pass filter data 
            tic
            EEG =  high_pass_filter(EEG,HP_f,fftflag);
            toc

            %plot after filtering without baseline correction
            myplot_F(EEG,[-100 300],[-500 500],['After filtering'],epoching_long);
            saveas(gcf, [current_output_folder '\after_filtering'])
            close

             % DownSample 
            EEG = pop_resample(EEG, DS);

            myplot_F(EEG,[-100 300],[-500 500],['After resampling'],epoching_long);
            saveas(gcf, [current_output_folder '\after_resampling'])
            close

            % Epoching data 
            EEG = pop_epoch(EEG, {}, epoching_long);

            %save figure
            myplot_F(EEG,[-100 300],[-500 500],['Step 1.2']);
            saveas(gcf,[current_output_folder '\step1_2'])
            close

            savename = [Block_cond{n_cond,1} '_step1_2'];
            EEG = pop_saveset( EEG, 'filename',savename,'filepath',current_output_folder);
    
        end
        disp('Step 1.2 ended')


%% STEP 2 SOUND
% 2.1 - prepare for Sound
        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2'];
        mkdir([output, '\' name_folder '\2_SOUND']);
        current_output_folder = [output, '\' name_folder '\2_SOUND'];

        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step2_2.set']) && global_overwrite==0 && sound_overwrite==0
                warning(['skipping step 2 for ' name_condition ' , file already present'])
        else
            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder);

            %function to construct leadfield and  to re-referencing 
            [EEG,LFM_sphere]= before_SOUND_G(EEG,ref_before_sound,chan_loc_st);

            % plot butterfly and save before SOUND     
            myplot_F(EEG,[-100 300],[-50 50],['Before SOUND']);
            saveas(gcf,[current_output_folder '\' name_condition '_before_SOUND'])
            close 

            %save dateset and leadfield matrix
            %savename = [name_condition '_beforeSound'];
            %EEG = pop_saveset(EEG, 'filename',savename,'filepath', current_output_folder);
            % save leadfield
            save([current_output_folder '\LFM_sphere.txt'], 'LFM_sphere', '-ascii'); 

            disp('Step 2.1  pre-SOUND ended');

            %%%% 2.2 SOUND

            %function for SOUND algorithm
            EEG_clean = SOUND_to_samples(EEG,LFM_sphere,lambda_value,iter);
            delete(gcp)
            close

            EEG = EEG_clean;

            % plot butterfly and save after SOUND
            myplot_F(EEG,[-100 300],[-50 50],['After SOUND']);
            saveas(gcf,[current_output_folder '\step2_2'])
            close

            %save set after SOUND
            EEG.setname= [Block_cond{n_cond,1} '_afterSOUND'];
            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1} '_step2_2'],'filepath', current_output_folder);
        end

disp(' step 2.2 SOUND ended')

%% 3 Automated Artifact rejection

        previous_folder= [output, '\' name_folder '\2_SOUND'];
        mkdir([output, '\' name_folder  '\3_Artifact_rejection']);
        current_output_folder = [output, '\' name_folder  '\3_Artifact_rejection'];

        if  isfile([ current_output_folder '\' Block_cond{n_cond,1}  '_step3.set']) && global_overwrite==0 && rej_overwrite==0
                warning(['skipping step 3 for ' name_condition  ' , file already present'])

        else
            EEG = pop_loadset('filename',[Block_cond{n_cond,1}  '_step2_2.set'],'filepath',previous_folder);

            % Automated Trial rejetion
            [EEG,rejected_epochs] = reject_artifact(EEG,SD_first,SD_last);  

            %save rejected epochs
            save([current_output_folder '\rejected_epochs'], 'rejected_epochs');

            %plot
            myplot_F(EEG,[-100 300],[-50 50],['After trial rejction']);
            saveas(gcf,[current_output_folder '\after_trialrejction'])
            close

            % save data after analysisG rejection:
            EEG.setname= 'trialrejction';
            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1}  '_step3'],'filepath', current_output_folder);
        end

disp('step 3 Artifact rejection ended')

%% 4 ICA 
% 4.1 RUN ICA
        previous_folder= [output, '\' name_folder  '\3_Artifact_rejection'];
        mkdir([output,'\' name_folder '\4_ICA\4_1_RUNICA']);
        current_output_folder = [output,'\' name_folder '\4_ICA\4_1_RUNICA'];

        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step4_1.set']) && global_overwrite==0 && ICA_overwrite==0
                warning(['skipping step 4.1 for ' name_condition ' , file already present'])

        else
                EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step3.set'],'filepath',  previous_folder);
                
                %EEG = pop_tesa_removedata( EEG,interp);

                % Ica decomposition 
                EEG = pop_runica(EEG, 'chanind', [1: size(EEG.data,1)], 'interupt','on');

                %save set with ICA all components
                EEG.setname= 'After ICA';
                savename = [Block_cond{n_cond,1} '_step4_1'];
                EEG = pop_saveset(EEG, 'filename',savename,'filepath', current_output_folder);    
        end
disp('Step 4.1 ICA ended')
    end
end
%%  4.2 ICA - Remove ocular artifat with ICA
for iSub = 1:n_subjects
    subj =  subjects{iSub}
    for n_cond =1:tot_cond
        
        load([output '\' subj '\0_Data\Block_cond.mat']);
        name_condition = [subj '_' Block_cond{n_cond,2}];
        name_folder= [subj '\' name_condition]

            previous_folder= [output,'\' name_folder '\4_ICA\4_1_RUNICA'];
            mkdir([output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact']);
            current_output_folder = [output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact'];

        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step4_2.set']) && global_overwrite==0 && ICA_remove==0
                warning(['skipping step 4.2 for ' name_condition ' , file already present'])
            else
             EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step4_1.set'],'filepath',  previous_folder);
             
             Select_ICA_components       
        end
    disp('Step 4.2 ICA selection components ended')
    end
end
%% 5. Use SSP-SIR to clean muscle artifacts
%  SSP-SIR
for iSub = 1:n_subjects
    subj =  subjects{iSub}
    for n_cond= 1:tot_cond
        load([output '\' subj '\0_Data\Block_cond.mat']);
        name_condition = [subj '_' Block_cond{n_cond,2}];
        name_folder= [subj '\' name_condition]

        previous_folder= [output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact'];
        mkdir([output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations']);
        current_output_folder = [output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations'];

        if  size(dir([current_output_folder '/*.set' ]),1)== 6 && global_overwrite==0 && SSP_SIR_step1==0
                warning(['skipping step 5.1 for ' name_condition ' , files already present'])
        else
            % re-load  LFM_sphere and EEG from previuos step
            LFM_sphere= load([ output '\' name_folder '\2_SOUND\LFM_sphere.txt']);
            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step4_2.set'],'filepath',  previous_folder);    
            
            EEG = pop_tesa_removedata( EEG,interp);

            % SSP-SIR with TEP with incremental number of removed components, up to 5
            SSP_SIR_components(EEG,LFM_sphere,SSP_SIR_window, current_output_folder,name_condition,rmBase,limx,limy)
        end

disp('Step 5.1 SSP-SIR and TEPs with incremental number of removed components ended ')
    end
end
  %%  5.2 SSP-SIR manual components selection
for iSub = 1:n_subjects
    subj =  subjects{iSub} 
    for n_cond= 1:tot_cond
        load([output '\' subj '\0_Data\Block_cond.mat']);
        name_condition = [subj '_' Block_cond{n_cond,2}];
        name_folder= [subj '\' name_condition]

        previous_folder= [output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations'];
        mkdir([output, '\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components']);
        current_output_folder = [output,'\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components'];

%         if  size(dir([current_output_folder '/*.mat' ]),1)==1 && global_overwrite==0 && SSP_SIR_step2==0
%                      warning(['skipping step 5.2 for ' subj ' , files already present'])
%         else
            decide1_SSPSIR_comp_MOBI
%         end
disp('Step 5.2 SSP-SIR removed components ended ')
    end
end
%% 6. Last step of processing after SSP-SIR
for iSub = 1:n_subjects
    subj =  subjects{iSub};
    for n_cond= 1: tot_cond
        load([output '\' subj '\0_Data\Block_cond.mat']);
        name_condition = [subj '_' Block_cond{n_cond,2}];
        name_folder= [subj '\' name_condition]
        savename= [Block_cond{n_cond,1}];

        previous_folder= [output,'\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components'];
        mkdir([output, '\' name_folder '\6_final_preprocessing_new_interp']);
        current_output_folder = [output,'\' name_folder '\6_final_preprocessing_new_interp'];
        
         if  isfile([current_output_folder '\' savename '_TEPs_def.set'])&& ...
                isfile([current_output_folder '\' savename  '_TEPs_def.mat'])
                     warning(['skipping step 6 for ' name_condition ' , files already present'])
        else        
            % to know the number of removed components in previuos step
            load([previous_folder '\SSPSIR_delcomponents.mat'])
            k = SSPSIR_delcomponents(1,1); 

            % load the set after SSP-SIR with correct number of removed variables
            EEG = pop_loadset('filename',[name_condition '_5SSP-SIR_' num2str(k) '.set'],'filepath', [output,'\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations']);
            EEG.comments=pop_comments(EEG.comments,'','SSP-SIR components evaluated by 2 scientists',1);

            %EEG = pop_tesa_interpdata( EEG, 'cubic',interp_Win);

            EEG= last_steps_MOBI_new(EEG,savename,LP_f,old_ref,chan_loc_st,rmBase,current_output_folder,epoching_short,reref);
            %myplot_F(EEG,[-100 400],[-15 15],['Final']);

            %artefact rejection 2 done in folder 6_final_preprocessing
            load([output,'/' name_folder '/6_final_preprocessing/artrej2_rejepochs.mat'])
            trialrej=zeros([1 size(EEG.data,3)]);
            trialrej(artrej2_rejepochs)=1;
            EEG = pop_rejepoch(EEG,trialrej,0);
            
            %remove data between 0 and 3
            EEG = pop_tesa_removedata( EEG,interp);
            myplot_F(EEG,[-100 400],[-20 40],['final remove']);
            saveas(gcf,[current_output_folder '/final_remove'])
            close

            %save
            EEG.setname= 'TEPs_def';
            EEG = pop_saveset( EEG, 'filename',[savename '_TEPs_def.set'],'filepath',current_output_folder);

            
            % convert to fieldtrip structure and save figure TEPs def
            TEP_final2=convert_to_fieldtrip(EEG,latency,current_output_folder,savename,ylimit);
        
            end
disp('Step 6 ended ')
    end
    

end