123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426 |
- %%% 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
|