From 883d27e530c2d0cf4e96a359efd8de04d2653f5a Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 7 Sep 2022 12:06:29 +0300 Subject: [PATCH 01/18] save more variables and adjust indents/outputs --- src/article/RewardsPreprocessing_withoutICA.m | 408 ++++++------------ 1 file changed, 143 insertions(+), 265 deletions(-) diff --git a/src/article/RewardsPreprocessing_withoutICA.m b/src/article/RewardsPreprocessing_withoutICA.m index 35263b4..5f79d55 100755 --- a/src/article/RewardsPreprocessing_withoutICA.m +++ b/src/article/RewardsPreprocessing_withoutICA.m @@ -1,161 +1,85 @@ - % %% PREPROCESSING - % %% Stage 1: Process data to determine noisy/faulty electrodes - % %% Step 1.1: Pre-ICA - clear all; - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data - filenames = dir('*.vhdr'); - filenames.name%Compile list of all data - - %for participant = 1:length(filenames) %Cycle through participants - for participant = 1:500 %Cycle through participants - - %Get participant name information - disp(['Participant: ', num2str(participant)]) %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - participant_varname = ['RewardProcessing_S1PostICA_',participant_number{2}]; %Create new file name - - %Load Data - EEG = []; %Clear past data - - [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data - - %Make it a 32 channels cap - reduces any participants with a 64 channel system - if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup - %a=EEG; - [EEG] = doRemoveChannels(EEG,{'AF3','AF4','AF7','AF8','C1','C2','C5','C6','CP3','CPz','F1','F2','F5','F6','FC3','FC4','FT10','FT7','FT8','FT9','Oz','P1','P2','P5','P6','PO3','PO4','PO7','PO8','TP7','TP8','CP4'},EEG.chanlocs); %Removes electrodes that are not part of the 32 channel system - %b=EEG - %AF3 AF4 AF7 AF8 C1 C2 C5 C6 CP3 CPz F1 F2 F5 F6 FC3 FC4 FT10 FT7 FT8 FT9 Oz P1 P2 P5 P6 PO3 PO4 PO7 PO8 TP7 TP8 CP4 - end - - %Re-Reference - load('/ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path - % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references - [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) - [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data - [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) - - %Filter - [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter data: Low cutoff 0.1, high cutoff 30, order 4, notch 60 - - %ICA - % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal - - %Save Output - save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention - end - %% Step 1.2: Post-ICA - clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data - %Find and change working folder to saved data from last for loop - filenames = dir('RewardProcessing_S1PostICA*'); %Compile list of all data - - %for participant = 1:length(filenames) %Cycle through participants - for participant = 1:length(filenames) - disp(['Participant: ', num2str(participant)]); %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-4),'_'); %Split filename into components - participant_varname = ['RewardProcessing_S1Final_',participant_number{3}]; %Create new file name - - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output - - %Inverse ICA - % ICAViewer %Script to navigate ICA loading topographic maps and loadings in order to decide which component(s) reflect eye blinks - % [EEG] = doICARemoveComponents(EEG,str2num(cell2mat(EEG.ICAcomponentsRemoved))); %Remove identified components and reconstruct EEG data - % - a2=EEG; - %Segment Data - [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) - c2=EEG; - %Baseline Correction - [EEG] = doBaseline(EEG,[-200,0]); %Baseline correction in ms - - %Artifact Rejection - [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 uV/ms gradient criteria - [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 uV max-min criteria - - save(participant_varname,'EEG'); %Save the current output - end - - %% Step 1.3: Determining faulty electrodes - % clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data - %Find and change working folder to saved data from last for loop - filenames = dir('RewardProcessing_S1Final_*'); %Compile list of all data - - for participant = 318:length(filenames) %Cycle through participants - disp(['Participant: ', num2str(participant)]); %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-4),'_'); %Split filename into components - - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output - - %Extract Artifact Rejection Output - AR_indx = []; trials_removed = []; %Clear past participant data - AR_indx = EEG.artifactPresent; %Reassign artifact rejection output - AR_indx(AR_indx > 1)=1; %Re-ID any indication of a rejected segment to be identified as the number 1 - trials_removed = sum(AR_indx')/size(AR_indx,2); %Determine the percentage of segments removed for each electrode - - %Determines the number of electrodes that surpassed different percentages of rejections - avg_removed(participant,1) = sum(trials_removed>.1); - avg_removed(participant,2) = sum(trials_removed>.1); - avg_removed(participant,3) = sum(trials_removed>.15); - avg_removed(participant,4) = sum(trials_removed>.2); - avg_removed(participant,5) = sum(trials_removed>.25); - avg_removed(participant,6) = sum(trials_removed>.3); - avg_removed(participant,7) = sum(trials_removed>.35); - avg_removed(participant,8) = sum(trials_removed>.4); - avg_removed(participant,9) = sum(trials_removed>.45); - avg_removed(participant,10) = sum(trials_removed>.5); - avg_removed(participant,11) = sum(trials_removed>.55); - avg_removed(participant,12) = sum(trials_removed>.6); - avg_removed(participant,13) = sum(trials_removed>.65); - avg_removed(participant,14) = sum(trials_removed>.7); - avg_removed(participant,15) = sum(trials_removed>.75); - avg_removed(participant,16) = sum(trials_removed>.8); - avg_removed(participant,17) = sum(trials_removed>.85); - avg_removed(participant,18) = sum(trials_removed>.9); - avg_removed(participant,19) = sum(trials_removed>.95); - - %Determine a level of rejection for each individual - % rejection_level = find(avg_removed(participant,:)<11); %Uses a percentage within which less than 11 electrodes have been removed - % rejection_level = rejection_level(1); %Uses the lowest rejection level within which less than 11 electrodes have been removed - % - %Save rejection information into a variable - numbers = [10,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95]/100; %First, determine the levels of rejection. 10 is repeated twice because it used ti be 5 but we decided this was too strict - p_chanreject{participant,1} = cellstr(participant_number{3}); %Insert participant number - p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed - p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed - end - - %Save channels to reject information into a mat file - save('Chans_rejected_auto500','p_chanreject'); +% % % %% PREPROCESSING +% % % %% Stage 1: Process data to determine noisy/faulty electrodes +% % % %% Step 1.1: Pre-ICA +% % clear all; +% % dirdata='data/Raw Data Part 1'; +% % cd(dirdata); %Find and change working folder to raw EEG data +% % filenames = dir('*.vhdr'); +% % +% % for participant = 1:500 %Cycle through participants +% % +% % %Get participant name information +% % disp(['Participant: ', num2str(participant)]) %Display current participant being processed +% % participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components +% % participant_varname = ['RewardProcessing_S1Post_',participant_number{2}]; %Create new file name +% % +% % %Load Data +% % EEG = []; %Clear past data +% % [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data +% % +% % %Make it a 32 channels cap - reduces any participants with a 64 channel system +% % if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup +% % [EEG] = doRemoveChannels(EEG,{'AF3','AF4','AF7','AF8','C1','C2','C5','C6','CP3','CPz','F1','F2','F5','F6','FC3','FC4','FT10','FT7','FT8','FT9','Oz','P1','P2','P5','P6','PO3','PO4','PO7','PO8','TP7','TP8','CP4'},EEG.chanlocs); %Removes electrodes that are not part of the 32 channel system +% % end +% % +% % %Re-Reference +% % load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path +% % % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references +% % [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) +% % [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data +% % [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) +% % +% % %Filter +% % [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter data: Low cutoff 0.1, high cutoff 30, order 4, notch 60 +% % +% % % %ICA +% % % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal +% % +% % %Segment Data +% % [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) +% % +% % %Baseline Correction +% % [EEG] = doBaseline(EEG,[-200,0]); %Baseline correction in ms +% % +% % %Artifact Rejection +% % [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 uV/ms gradient criteria +% % [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 uV max-min criteria +% % +% % %Extract Artifact Rejection Output +% % AR_indx = []; trials_removed = []; %Clear past participant data +% % AR_indx = EEG.artifactPresent; %Reassign artifact rejection output +% % AR_indx(AR_indx > 1)=1; %Re-ID any indication of a rejected segment to be identified as the number 1 +% % trials_removed = sum(AR_indx')/size(AR_indx,2); %Determine the percentage of segments removed for each electrode +% % +% % %Save rejection information into a variable: electrodes with a trial rejection rate greater than 40% were tagged for removal +% % p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number +% % p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed +% % p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed +% % +% % %Save Output +% % % save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention +% % end +% % %Save channels to reject information into a mat file +% % save('Chans_rejected500.mat','p_chanreject'); + %% Stage 2: Process data for analysis - %% Step 2.1: Pre-ICA clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data +% dirdata='data/Raw Data Part 1'; +% cd(dirdata); %Find and change working folder to raw EEG data filenames = dir('*.vhdr'); %Compile list of all data +load('Chans_rejected500.mat'); %Load list of electrodes to remove - for participant = 1:length(filenames) %Cycle through participants + for participant = 1:500 %Cycle through participants %Get participant name information disp(['Participant: ', num2str(participant)]); %Display current participant being processed participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - participant_varname = ['RewardProcessing_S2PostICA_',participant_number{2}]; %Create new file name + participant_varname = ['RewardProcessing_S2Post_',participant_number{2}]; %Create new file name %Load Data EEG = []; %Clear past data - [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data + [EEG] = doLoadBVData(filenames(participant).name); %Load raw EEG data %Make it a 32 channels cap - reduces any participants with a 64 channel system if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup @@ -163,47 +87,55 @@ end %Re-Reference - load('/ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path - % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references - try - %Note:Through visual inspection, some of the reference electrodes were deemed poor quality and were not used - if str2num(participant_number{2}) == 164 || str2num(participant_number{2}) == 255 || str2num(participant_number{2}) == 277 %Participants with one noisy mastoid electrode - [EEG] = doRereference(EEG,{'TP9'},{'ALL'},EEG.chanlocs); %Re-reference with only TP9 - elseif str2num(participant_number{2}) == 110 %Participants with one noisy mastoid electrode - [EEG] = doRereference(EEG,{'TP10'},{'ALL'},EEG.chanlocs); %Re-reference with only TP10 - else - [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) - end + load('ChanlocsMaster.mat'); + [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) - % %Remove Faulty Electrodes - % p_chans = []; %Clear past participants electrode indices - % p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove - % chans_to_remove = []; %Clear past participants electrode labels - % if ~isempty(p_chans{1}) %Determine whether electrodes need to be removed - % x = 1; %Begin index count - % for l = 1:size(p_chans,2) %Scroll through the electrodes - % if str2num(p_chans{1,l}) ~= 26 %Ensure that FCz id not removed - % chans_to_remove{1,x} = chanlocsMaster(str2num(p_chans{1,l})).labels; %Create list of electrode labels to remove - % x = x+1; %Increase index count - % end - % end - % [EEG] = doRemoveChannels(EEG,chans_to_remove,chanlocsMaster); %Remove electrodes - % else %If there is no electrodes to remove - % [EEG] = doRemoveChannels(EEG,{},chanlocsMaster); %Still run function, but do not remove any electrodes - % end + %Load channel location file, please make sure the location of this file is in your path + % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references + + % This manual removal has been removed from this code as well as in other software codes + %Note:Through visual inspection, some of the reference electrodes were deemed poor quality and were not used +% if str2num(participant_number{2}) == 164 || str2num(participant_number{2}) == 255 || str2num(participant_number{2}) == 277 %Participants with one noisy mastoid electrode +% [EEG] = doRereference(EEG,{'TP9'},{'ALL'},EEG.chanlocs); %Re-reference with only TP9 +% elseif str2num(participant_number{2}) == 110 %Participants with one noisy mastoid electrode +% [EEG] = doRereference(EEG,{'TP10'},{'ALL'},EEG.chanlocs); %Re-reference with only TP10 +% else +% [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) +% end + + %Remove faulty electrodes + p_chans = []; %Clear past participants electrode indices + p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove + chans_to_remove = []; %Clear past participants electrode labels + if ~isempty(p_chans{1}) %Determine whether electrodes need to be removed + x = 1; %Begin index count + for l = 1:size(p_chans,2) %Scroll through the electrodes + if str2num(p_chans{1,l}) ~= 26 %Ensure that FCz is not removed + chans_to_remove{1,x} = chanlocsMaster(str2num(p_chans{1,l})).labels; %Create list of electrode labels to remove + x = x+1; %Increase index count + end + end + [EEG] = doRemoveChannels(EEG,chans_to_remove,chanlocsMaster); %Remove electrodes + else %If there is no electrodes to remove + [EEG] = doRemoveChannels(EEG,{},chanlocsMaster); %Still run function, but do not remove any electrodes + end + + + % This manual removal has been removed from this code as well as in other software codes %Removed missed faulty electrodes on certain participants (determined on Step 1.2) - This was determined via visual inspection of the data after step 1.2 - if str2num(participant_number{2})==77 - [EEG] = doRemoveChannels(EEG,{'P8','CP6','AFz'},EEG.chanlocs); - elseif str2num(participant_number{2})==199 - [EEG] = doRemoveChannels(EEG,{'FC1','O2'},EEG.chanlocs); - elseif str2num(participant_number{2})==486 - [EEG] = doRemoveChannels(EEG,{'FC2'},EEG.chanlocs); - else - %Skip - end +% if str2num(participant_number{2})==77 +% [EEG] = doRemoveChannels(EEG,{'P8','CP6','AFz'},EEG.chanlocs); +% elseif str2num(participant_number{2})==199 +% [EEG] = doRemoveChannels(EEG,{'FC1','O2'},EEG.chanlocs); +% elseif str2num(participant_number{2})==486 +% [EEG] = doRemoveChannels(EEG,{'FC2'},EEG.chanlocs); +% else +% %Skip +% end + %Filter [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter da ta: Low cutoff 0.1, high cutoff 30, order 4, notch 60 @@ -211,56 +143,11 @@ %ICA % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal - %Save output - save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention - catch - end - end - %% Step 2.2: Post-ICA - clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data - %Find and change working folder to saved data from last for loop - filenames = dir('RewardProcessing_S2PostICA*'); %Compile list of all data - - for participant = 1:length(filenames) %Cycle through participants - disp(['Participant: ', num2str(participant)]); %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-4),'_'); %Split filename into components - participant_varname = ['RewardProcessing_S2PostInvICA_',participant_number{3}]; %Create new file name - - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output - - % %Remove ICA Component - % ICAViewer %Script to navigate ICA loading topographic maps and loadings in order to decide which component(s) reflect eye blinks - % [EEG] = doICARemoveComponents(EEG,str2num(cell2mat(EEG.ICAcomponentsRemoved))); %Remove identified components and reconstruct EEG data - % - save(participant_varname,'EEG'); %Save the current output - end - %% Step 2.3: Final - clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; - - cd(dirdata); %Find and change working folder to raw EEG data - filenames = dir('RewardProcessing_S2PostICA*'); %Compile list of all data - - for participant = 1:length(filenames) %Cycle through participants - - disp(['Participant: ', num2str(participant)]); %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-4),'_'); %Split filename into components - participant_varname = ['RewardProcessing_S2Final_',participant_number{4}]; %Create new file name - - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output - %Topograpic Interpolation - load('/ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path + load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate electrodes which were previously removed - %Determine markers of interest + %Determine markers of interest markers = {'S110','S111'}; %Loss, win %Segment Data @@ -276,44 +163,35 @@ %ERP try - [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses + [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses + All_ERP(:,:,:,participant) = EEG.ERP.data; %Store all the ERP data into a single variable + All_trials(1,participant) = EEG.ERP.epochCount(1); %Store the number of trials + All_trials(2,participant) = EEG.ERP.epochCount(2); + All_trials(3,participant) = EEG.ERP.totalEpochs; catch - msgbox('Hello'); - end - save(participant_varname,'EEG'); %Save the current output - end - %% EXTRACTION OF DATA %% - %% Aggregate data across participants - clear all; close all; clc; %First, clean the environment - dirdata='/Raw Data Part 1'; + msgbox('No remained data for this participant after trial rejection'); + All_ERP(:,:,:,participant) = []; + All_trials(1,participant) = 0; + All_trials(2,participant) = 0; + All_trials(3,participant) = 0; - cd(dirdata); %Find and change working folder to raw EEG data - %Find and change working folder to saved data from last for loop - filenames = dir('RewardProcessing_S2Final*'); %Compile list of all data - - for participant = 1:length(filenames) %Cycle through participants - disp(['Participant: ', num2str(participant)]); %Display current participant being processed - - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output - - All_ERP(:,:,:,participant) = EEG.ERP.data; %Store all the ERP data into a single variable + end + %Save output + save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention end - + save('Ref_All_ERP', 'All_ERP'); %Save ERP Data - All_ERP=All_ERP(:,151:750,:,:); - - % %% RewP_Waveforms - csvwrite('Ref_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(26,:,1,:)),2),nanmean(squeeze(All_ERP(26,:,2,:)),2),nanmean(squeeze(All_ERP(26,:,1,:)),2)-nanmean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. - % %% RewP_Waveforms_AllPs - tt1=squeeze(All_ERP(26,:,1,:)); - tt2=squeeze(All_ERP(26,:,2,:)); - - csvwrite('Ref_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. - % %% RewP_Latency - [~,peak_loc] = max(squeeze(All_ERP(26,226:276,1,:))-squeeze(All_ERP(26,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. - peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds - peak_loc(toberemoved)=[]; - csvwrite('Ref_RewP_Latency.csv',peak_loc'); %Export data +% All_ERP=All_ERP(:,151:750,:,:); +% % %% RewP_Waveforms +% csvwrite('Ref_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(26,:,1,:)),2),nanmean(squeeze(All_ERP(26,:,2,:)),2),nanmean(squeeze(All_ERP(26,:,1,:)),2)-nanmean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +% % %% RewP_Waveforms_AllPs +% tt1=squeeze(All_ERP(26,:,1,:)); +% tt2=squeeze(All_ERP(26,:,2,:)); +% +% csvwrite('Ref_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +% % %% RewP_Latency +% [~,peak_loc] = max(squeeze(All_ERP(26,226:276,1,:))-squeeze(All_ERP(26,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. +% peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds +% peak_loc(toberemoved)=[]; +% csvwrite('Ref_RewP_Latency.csv',peak_loc'); %Export data From b846713652d47446f0d1ca8875d9f94ab864b1ec Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 7 Sep 2022 15:33:52 +0300 Subject: [PATCH 02/18] add the code without channel detection step --- src/article/RewardsPreprocessing_withoutICA.m | 131 +++++++++--------- 1 file changed, 68 insertions(+), 63 deletions(-) diff --git a/src/article/RewardsPreprocessing_withoutICA.m b/src/article/RewardsPreprocessing_withoutICA.m index 5f79d55..c83f3e1 100755 --- a/src/article/RewardsPreprocessing_withoutICA.m +++ b/src/article/RewardsPreprocessing_withoutICA.m @@ -1,66 +1,66 @@ -% % % %% PREPROCESSING -% % % %% Stage 1: Process data to determine noisy/faulty electrodes -% % % %% Step 1.1: Pre-ICA -% % clear all; -% % dirdata='data/Raw Data Part 1'; -% % cd(dirdata); %Find and change working folder to raw EEG data -% % filenames = dir('*.vhdr'); -% % -% % for participant = 1:500 %Cycle through participants -% % -% % %Get participant name information -% % disp(['Participant: ', num2str(participant)]) %Display current participant being processed -% % participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components -% % participant_varname = ['RewardProcessing_S1Post_',participant_number{2}]; %Create new file name -% % -% % %Load Data -% % EEG = []; %Clear past data -% % [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data -% % -% % %Make it a 32 channels cap - reduces any participants with a 64 channel system -% % if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup -% % [EEG] = doRemoveChannels(EEG,{'AF3','AF4','AF7','AF8','C1','C2','C5','C6','CP3','CPz','F1','F2','F5','F6','FC3','FC4','FT10','FT7','FT8','FT9','Oz','P1','P2','P5','P6','PO3','PO4','PO7','PO8','TP7','TP8','CP4'},EEG.chanlocs); %Removes electrodes that are not part of the 32 channel system -% % end -% % -% % %Re-Reference -% % load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path -% % % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references -% % [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) -% % [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data -% % [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) -% % -% % %Filter -% % [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter data: Low cutoff 0.1, high cutoff 30, order 4, notch 60 -% % -% % % %ICA -% % % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal -% % -% % %Segment Data -% % [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) -% % -% % %Baseline Correction -% % [EEG] = doBaseline(EEG,[-200,0]); %Baseline correction in ms -% % -% % %Artifact Rejection -% % [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 uV/ms gradient criteria -% % [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 uV max-min criteria -% % -% % %Extract Artifact Rejection Output -% % AR_indx = []; trials_removed = []; %Clear past participant data -% % AR_indx = EEG.artifactPresent; %Reassign artifact rejection output -% % AR_indx(AR_indx > 1)=1; %Re-ID any indication of a rejected segment to be identified as the number 1 -% % trials_removed = sum(AR_indx')/size(AR_indx,2); %Determine the percentage of segments removed for each electrode -% % -% % %Save rejection information into a variable: electrodes with a trial rejection rate greater than 40% were tagged for removal -% % p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number -% % p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed -% % p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed -% % -% % %Save Output -% % % save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention -% % end -% % %Save channels to reject information into a mat file -% % save('Chans_rejected500.mat','p_chanreject'); + % %% PREPROCESSING + % %% Stage 1: Process data to determine noisy/faulty electrodes + % %% Step 1.1: Pre-ICA + clear all; + dirdata='data/Raw Data Part 1'; + cd(dirdata); %Find and change working folder to raw EEG data + filenames = dir('*.vhdr'); + + for participant = 1:500 %Cycle through participants + + %Get participant name information + disp(['Participant: ', num2str(participant)]) %Display current participant being processed + participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components + participant_varname = ['RewardProcessing_S1Post_',participant_number{2}]; %Create new file name + + %Load Data + EEG = []; %Clear past data + [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data + + %Make it a 32 channels cap - reduces any participants with a 64 channel system + if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup + [EEG] = doRemoveChannels(EEG,{'AF3','AF4','AF7','AF8','C1','C2','C5','C6','CP3','CPz','F1','F2','F5','F6','FC3','FC4','FT10','FT7','FT8','FT9','Oz','P1','P2','P5','P6','PO3','PO4','PO7','PO8','TP7','TP8','CP4'},EEG.chanlocs); %Removes electrodes that are not part of the 32 channel system + end + + %Re-Reference + load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path + % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references + [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) + [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data + [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) + + %Filter + [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter data: Low cutoff 0.1, high cutoff 30, order 4, notch 60 + +% %ICA + % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal + + %Segment Data + [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) + + %Baseline Correction + [EEG] = doBaseline(EEG,[-200,0]); %Baseline correction in ms + + %Artifact Rejection + [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 uV/ms gradient criteria + [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 uV max-min criteria + + %Extract Artifact Rejection Output + AR_indx = []; trials_removed = []; %Clear past participant data + AR_indx = EEG.artifactPresent; %Reassign artifact rejection output + AR_indx(AR_indx > 1)=1; %Re-ID any indication of a rejected segment to be identified as the number 1 + trials_removed = sum(AR_indx')/size(AR_indx,2); %Determine the percentage of segments removed for each electrode + + %Save rejection information into a variable: electrodes with a trial rejection rate greater than 40% were tagged for removal + p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number + p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed + p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed + + %Save Output +% save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention + end + %Save channels to reject information into a mat file + save('Chans_rejected500.mat','p_chanreject'); %% Stage 2: Process data for analysis clear all; close all; clc; %First, clean the environment @@ -109,6 +109,9 @@ %Remove faulty electrodes p_chans = []; %Clear past participants electrode indices p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove + % save the number of channels removed/interpolated + All_rejChan(participant)=length(p_chans); + chans_to_remove = []; %Clear past participants electrode labels if ~isempty(p_chans{1}) %Determine whether electrodes need to be removed x = 1; %Begin index count @@ -181,6 +184,8 @@ end save('Ref_All_ERP', 'All_ERP'); %Save ERP Data + save('All_trials', 'All_trials'); %Save ERP Data + save('All_rejChan', 'All_rejChan'); %Save ERP Data % All_ERP=All_ERP(:,151:750,:,:); % % %% RewP_Waveforms From 8e60a1af78c854cec67b0c968699794d71a25710 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 21:59:34 +0300 Subject: [PATCH 03/18] fix paths and indentations --- src/eeglab/eeglab_preprocessing.m | 273 +++++++++++++++--------------- 1 file changed, 137 insertions(+), 136 deletions(-) diff --git a/src/eeglab/eeglab_preprocessing.m b/src/eeglab/eeglab_preprocessing.m index 644a233..c44ab4f 100644 --- a/src/eeglab/eeglab_preprocessing.m +++ b/src/eeglab/eeglab_preprocessing.m @@ -1,78 +1,77 @@ -% % --- This is the preprocessing workflow reproduced by EEGLAB functions--- -% % For contact: aya.kabbara7@gmail.com - - cd('/Raw Data'); %Find and change working folder to raw EEG data - -%% === 1st pass: Detect bad channels ==== - -filenames = dir('*.vhdr') -nb=500; -for participant = 1:nb %Cycle through participants - - % Get participant name information - disp(['Participant: ', num2str(participant)]) %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - - EEG = pop_loadbv('/Raw Data Part 1/', filenames(participant).name); - % add the channel location file - EEG=pop_chanedit(EEG, 'lookup','/eeglab2021.1/functions/supportfiles/Standard-10-20-Cap81.ced'); - - % % reduce into 32 channels - - try - if(EEG.nbchan > 32) - EEG = pop_select(EEG, 'channel',{'Fp1' - 'Fz' - 'F3' - 'F7' - 'FC5' - 'FC1' - 'Cz' - 'C3' - 'T7' - 'CP5' - 'CP1' - 'Pz' - 'P3' - 'P7' - 'O1' - 'POz' - 'O2' - 'P4' - 'P8' - 'CP6' - 'CP2' - 'C4' - 'T8' - 'FC6' - 'FC2' - 'FCz' - 'F4' - 'F8' - 'Fp2' - 'TP10' - 'TP9'}); - end - catch - disp('Problem'); - rm_channels{participant}=''; - continue; - end - -EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',30); -[eeg2,HP,BUR,removed_channels] = clean_artifacts(EEG); -rm_channels{participant}=find(removed_channels); -end -%% === Save the detected bad channels ==== -save('eeglab_removedchans500','rm_channels'); +% % % --- This is the preprocessing workflow reproduced by EEGLAB functions--- +% % % For contact: aya.kabbara7@gmail.com + +% %% === 1st pass: Detect bad channels ==== +% +% filenames = dir('*.vhdr') +% nb=500; +% for participant = 1:nb %Cycle through participants +% +% % Get participant name information +% disp(['Participant: ', num2str(participant)]) %Display current participant being processed +% participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components +% % Load EEG data +% EEG = pop_loadbv(pwd, filenames(participant).name); +% +% % Add the channel location file +% EEG=pop_chanedit(EEG, 'lookup','eeglab2022.1/functions/supportfiles/Standard-10-20-Cap81.ced'); +% +% % Reduce into 32 channels +% try +% if(EEG.nbchan > 32) +% EEG = pop_select(EEG, 'channel',{'Fp1' +% 'Fz' +% 'F3' +% 'F7' +% 'FC5' +% 'FC1' +% 'Cz' +% 'C3' +% 'T7' +% 'CP5' +% 'CP1' +% 'Pz' +% 'P3' +% 'P7' +% 'O1' +% 'POz' +% 'O2' +% 'P4' +% 'P8' +% 'CP6' +% 'CP2' +% 'C4' +% 'T8' +% 'FC6' +% 'FC2' +% 'FCz' +% 'F4' +% 'F8' +% 'Fp2' +% 'TP10' +% 'TP9'}); +% end +% catch +% disp('Problem'); +% rm_channels{participant}=''; +% continue; +% end +% +% %Filter data +% EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',30); +% %Detect bad channels +% [eeg2,HP,BUR,removed_channels] = clean_artifacts(EEG); +% rm_channels{participant}=find(removed_channels); +% end +% %% === Save the detected bad channels ==== +% save('eeglab_removedchans500','rm_channels'); %% === 2nd pass: process data ==== clear all; clc; load('eeglab_removedchans500'); -cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data filenames = dir('*.vhdr') nb=500; for participant = 1:nb %Cycle through participants @@ -82,57 +81,55 @@ participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components %% === Read data ==== - - EEG = pop_loadbv('/Raw Data Part 1/', filenames(participant).name); + EEG = pop_loadbv(pwd, filenames(participant).name); %% === Add the channel location file ==== - EEG=pop_chanedit(EEG, 'lookup','/eeglab2021.1/functions/supportfiles/Standard-10-20-Cap81.ced'); - + EEG=pop_chanedit(EEG, 'lookup','eeglab2022.1/functions/supportfiles/Standard-10-20-Cap81.ced'); %% === reduce into 32 channels ==== try if(EEG.nbchan > 32) - EEG = pop_select(EEG, 'channel',{'Fp1' - 'Fz' - 'F3' - 'F7' - 'FC5' - 'FC1' - 'Cz' - 'C3' - 'T7' - 'CP5' - 'CP1' - 'Pz' - 'P3' - 'P7' - 'O1' - 'POz' - 'O2' - 'P4' - 'P8' - 'CP6' - 'CP2' - 'C4' - 'T8' - 'FC6' - 'FC2' - 'FCz' - 'F4' - 'F8' - 'Fp2' - 'TP10' - 'TP9'}); + EEG = pop_select(EEG, 'channel',{'Fp1' + 'Fz' + 'F3' + 'F7' + 'FC5' + 'FC1' + 'Cz' + 'C3' + 'T7' + 'CP5' + 'CP1' + 'Pz' + 'P3' + 'P7' + 'O1' + 'POz' + 'O2' + 'P4' + 'P8' + 'CP6' + 'CP2' + 'C4' + 'T8' + 'FC6' + 'FC2' + 'FCz' + 'F4' + 'F8' + 'Fp2' + 'TP10' + 'TP9'}); end - %% === Interpolate bad channels ==== - EEG = pop_interp(EEG, rm_channels{participant}, 'spherical'); +% % %% === Interpolate bad channels ==== +% % EEG = pop_interp(EEG, rm_channels{participant}, 'spherical'); +% % All_rejChan(participant)=length(rm_channels{participant}); %% === Re-referencing ==== - EEG=pop_chanedit(EEG, 'seteeglab',{'1:31','TP10 TP9'}); + EEG = pop_chanedit(EEG, 'seteeglab',{'1:31','TP10 TP9'}); EEG = pop_reref( EEG ,{'TP9','TP10'}); - %% === Filtering ==== EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',30); @@ -141,43 +138,47 @@ [EEG] = doSegmentData(EEG,markers,[-500 1298]); %Segment Data (S110 = Loss, S111 = Win). The segment window of interest is -200 to 1000ms, and we here add 300 ms before and after this because time-frequency edge artifacts (this is different than the first pass because we were being more conservative then) %% === Baseline adjustment ==== - EEG = pop_rmbase( EEG, [-200/1000,0]); + + %% === Bad trial identification and removal ==== + EEG = pop_eegthresh(EEG,1,[1:29] ,-50,50,0,1.2,0,1); + %% === ERP ==== try - %% === Bad trial identification and removal ==== - EEG = pop_eegthresh(EEG,1,[1:29] ,-50,50,-0.2,1.2,0,1); - %ERP - [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses - All_ERP_eeglab(:,:,:,participant) = EEG.ERP.data; %Store all the ERP data into a single variable -% save(participant_varname,'EEG'); %Save the current output - - catch - - continue + [EEG.ERP] = doERP(EEG,markers,0); + All_ERP(:,:,:,participant) = EEG.ERP.data; %Store all the ERP data into a single variable + All_trials(1,participant) = EEG.ERP.epochCount(1); %Store the number of trials + All_trials(2,participant) = EEG.ERP.epochCount(2); + All_trials(3,participant) = EEG.ERP.totalEpochs; + + catch + All_trials(1,participant) = 0; %Store the number of trials + All_trials(2,participant) = 0; + All_trials(3,participant) = 0; end catch - continue - + continue end end %% === Save variables and csv files ==== -save('All_ERP_eeglab', 'All_ERP_eeglab'); %Save ERP Data - -All_ERP=All_ERP_eeglab(:,151:750,:,:); - -chanOfInterest=17; -% channel of interest 17 is FCz -%% RewP_Waveforms_AllPs -tt1=squeeze(All_ERP(chanOfInterest,:,1,:)); -tt2=squeeze(All_ERP(chanOfInterest,:,2,:)); -% %% RewP_Waveforms -csvwrite('eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -csvwrite('eeglab_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. -% %% RewP_Latency -[~,peak_loc] = max(squeeze(All_ERP(chanOfInterest,226:276,1,:))-squeeze(All_ERP(chanOfInterest,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. -peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds -csvwrite('eeglab_RewP_Latency.csv',peak_loc'); %Export data +save('EEGLAB_All_ERP_withoutchans', 'All_ERP'); %Save ERP Data +save('EEGLAB_All_trials_withoutchans', 'All_trials'); %Save ERP Data +% save('EEGLAB_All_rejChan_withoutchans300', 'All_rejChan'); %Save ERP Data + +% All_ERP=All_ERP_eeglab(:,151:750,:,:); +% +% chanOfInterest=17; +% % channel of interest 17 is FCz +% %% RewP_Waveforms_AllPs +% tt1=squeeze(All_ERP(chanOfInterest,:,1,:)); +% tt2=squeeze(All_ERP(chanOfInterest,:,2,:)); +% % %% RewP_Waveforms +% csvwrite('eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +% csvwrite('eeglab_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +% % %% RewP_Latency +% [~,peak_loc] = max(squeeze(All_ERP(chanOfInterest,226:276,1,:))-squeeze(All_ERP(chanOfInterest,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. +% peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds +% csvwrite('eeglab_RewP_Latency.csv',peak_loc'); %Export data From 1a8a6a81f7b0dc200d8a47377f5d1d82bd4721b0 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 22:06:21 +0300 Subject: [PATCH 04/18] changed to permutation test --- results/README.md | 1 - src/eeglab/eeglab_stats.m | 266 ++++++-------------------------------- 2 files changed, 43 insertions(+), 224 deletions(-) delete mode 100644 results/README.md diff --git a/results/README.md b/results/README.md deleted file mode 100644 index c52f6ec..0000000 --- a/results/README.md +++ /dev/null @@ -1 +0,0 @@ -"2" at the end of a filename ==> results for the peak to peak 200 microV diff --git a/src/eeglab/eeglab_stats.m b/src/eeglab/eeglab_stats.m index 41977f9..5c62972 100644 --- a/src/eeglab/eeglab_stats.m +++ b/src/eeglab/eeglab_stats.m @@ -1,98 +1,40 @@ close all -clear - +% clear nbchan=29; nbpt=600; -P_Bonf=0.01/nbchan/nbpt; - -% % Section 1: Compute the statistical map(channels x time) showing the difference between EEGLAB and paper results in terms of conditional erp - -StatDiff=zeros(nbchan,nbpt); - -% orderEEGlab adjusts the order of channels to be like the one used by -% Williams et al - -orderEEglab=[8 -22 -7 -11 -21 -10 -20 -3 -27 -4 -28 -2 -6 -25 -5 -24 -26 -1 -29 -15 -17 -13 -18 -14 -19 -12 -16 -9 -23]; - - -load('All_ERP_ref.mat') -load('All_ERP_eeglab.mat') -All_ERP2=All_ERP_eeglab; - -% re-order channs for eeglab erps -All_ERP2=All_ERP2(orderEEglab,151:750,:,:); -All_ERP=All_ERP(:,151:750,:,:); - -tt1=squeeze(All_ERP2(26,:,1,:)); -tt2=squeeze(All_ERP2(26,:,2,:)); -% % remove zeros and nan -idx1 = isnan(tt1) ; -[r1,c1]=find(tt1==0); -[r1,c3]=find(idx1); - -idx2 = isnan(tt2) ; -[r2,c2]=find(tt2==0); -[r1,c4]=find(idx2); - -toberemoved=unique([unique(c1) ;unique(c2); unique(c3); unique(c4)]); -tKept1=[];tKept2=[]; -kept_eeglab=[]; -for su=1:500 - if(length(find(toberemoved==su))==0) -% % the subject su is present - kept_eeglab(end+1)=su; - tKept1(:,end+1)=tt1(:,su); - tKept2(:,end+1)=tt2(:,su); +% % Section 1: Compute the statistical map(channels x time) showing the difference between FieldTrip and paper results in terms of conditional erp +load('results/Ref_All_ERP.mat') +All_ERP1=All_ERP; +load('results/eeglab_All_ERP_samePipe.mat') +All_ERP2=All_ERP; +All_ERP1=All_ERP1(:,151:750,:,:); +All_ERP2=All_ERP2(:,151:750,:,:); + +condition_name={'Win','Loss','Difference'}; +matp=zeros(29,600,3); +for cond=1:3 + for ch=1:29 + ch + if(cond==3) + x1=squeeze(All_ERP1(ch,:,1,:))-squeeze(All_ERP1(ch,:,2,:)); + x2=squeeze(All_ERP2(ch,:,1,:))-squeeze(All_ERP2(ch,:,2,:)); + else + x1=squeeze(All_ERP1(ch,:,cond,:)); + x2=squeeze(All_ERP2(ch,:,cond,:)); + end + [clusters, p_values, t_sums, permutation_distribution ] = permutest(x1,x2,false,0.005,10000,true); + + for cl=1:length(clusters) + if(p_values(cl)<0.005) + significant_samples=clusters{cl}; + matp(ch,clusters{cl},cond)=1; + end + end end end -All_ERP2=All_ERP2(:,:,:,kept_eeglab); - - - -condition_name={'Win','Loss'}; - -for condition=1:2 - StatDiff=zeros(nbchan,nbpt); - - for channel=1:nbchan - for t=1:nbpt - Presults(channel,t)=ranksum(squeeze(All_ERP(channel,t,condition,:)),squeeze(All_ERP2(channel,t,condition,:))); - if(Presults(channel,t) Date: Tue, 11 Oct 2022 22:43:04 +0300 Subject: [PATCH 05/18] changed to permutation test --- src/eeglab/eeglab_stats.m | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/eeglab/eeglab_stats.m b/src/eeglab/eeglab_stats.m index 5c62972..5292a52 100644 --- a/src/eeglab/eeglab_stats.m +++ b/src/eeglab/eeglab_stats.m @@ -13,6 +13,10 @@ condition_name={'Win','Loss','Difference'}; matp=zeros(29,600,3); + +% % Section 2: Compute the cluster-based permutation maps for each +% condition- gain, loss, difference + for cond=1:3 for ch=1:29 ch @@ -74,6 +78,7 @@ 0 1 0 0 0 1 66/255 66/255 66/255]; +% % Section 2: Plot results for condition=1:3 StatDiff=matp(:,:,condition); From 96dd787a44290d3f91670e004347ce68a291808f Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 22:47:03 +0300 Subject: [PATCH 06/18] changed --- src/eeglab/eeglab_preprocessing.m | 202 ++++++++++++++++-------------- 1 file changed, 110 insertions(+), 92 deletions(-) diff --git a/src/eeglab/eeglab_preprocessing.m b/src/eeglab/eeglab_preprocessing.m index c44ab4f..6554b02 100644 --- a/src/eeglab/eeglab_preprocessing.m +++ b/src/eeglab/eeglab_preprocessing.m @@ -1,77 +1,91 @@ -% % % --- This is the preprocessing workflow reproduced by EEGLAB functions--- -% % % For contact: aya.kabbara7@gmail.com - -% %% === 1st pass: Detect bad channels ==== -% -% filenames = dir('*.vhdr') -% nb=500; -% for participant = 1:nb %Cycle through participants -% -% % Get participant name information -% disp(['Participant: ', num2str(participant)]) %Display current participant being processed -% participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components -% % Load EEG data -% EEG = pop_loadbv(pwd, filenames(participant).name); -% -% % Add the channel location file -% EEG=pop_chanedit(EEG, 'lookup','eeglab2022.1/functions/supportfiles/Standard-10-20-Cap81.ced'); -% -% % Reduce into 32 channels -% try -% if(EEG.nbchan > 32) -% EEG = pop_select(EEG, 'channel',{'Fp1' -% 'Fz' -% 'F3' -% 'F7' -% 'FC5' -% 'FC1' -% 'Cz' -% 'C3' -% 'T7' -% 'CP5' -% 'CP1' -% 'Pz' -% 'P3' -% 'P7' -% 'O1' -% 'POz' -% 'O2' -% 'P4' -% 'P8' -% 'CP6' -% 'CP2' -% 'C4' -% 'T8' -% 'FC6' -% 'FC2' -% 'FCz' -% 'F4' -% 'F8' -% 'Fp2' -% 'TP10' -% 'TP9'}); -% end -% catch -% disp('Problem'); -% rm_channels{participant}=''; -% continue; -% end -% -% %Filter data -% EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',30); -% %Detect bad channels -% [eeg2,HP,BUR,removed_channels] = clean_artifacts(EEG); -% rm_channels{participant}=find(removed_channels); -% end -% %% === Save the detected bad channels ==== -% save('eeglab_removedchans500','rm_channels'); +% % --- This is the preprocessing workflow reproduced by EEGLAB functions--- +% % For contact: aya.kabbara7@gmail.com -%% === 2nd pass: process data ==== +cd('tools/eeglab2022.1'); +%% === 1st pass: Detect bad channels ==== + +filenames = dir('../../data/*.vhdr') +nb=500; +for participant = 1:nb %Cycle through participants + + % Get participant name information + disp(['Participant: ', num2str(participant)]) %Display current participant being processed + participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components + % Load EEG data + EEG = pop_loadbv('../../data/', filenames(participant).name); + + % Add the channel location file + EEG=pop_chanedit(EEG, 'lookup','functions/supportfiles/Standard-10-20-Cap81.ced'); + + % Reduce into 32 channels + try + if(EEG.nbchan > 32) + EEG = pop_select(EEG, 'channel',{'Fp1' + 'Fz' + 'F3' + 'F7' + 'FC5' + 'FC1' + 'Cz' + 'C3' + 'T7' + 'CP5' + 'CP1' + 'Pz' + 'P3' + 'P7' + 'O1' + 'POz' + 'O2' + 'P4' + 'P8' + 'CP6' + 'CP2' + 'C4' + 'T8' + 'FC6' + 'FC2' + 'FCz' + 'F4' + 'F8' + 'Fp2' + 'TP10' + 'TP9'}); + end + catch + disp('Problem'); + rm_channels{participant}=''; + end + + %Filter data + EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',30); + + %% === Segmentation ==== + markers = {'S110','S111'}; %Loss, win + [EEG] = doSegmentData(EEG,markers,[-500 1298]); %Segment Data (S110 = Loss, S111 = Win). The segment window of interest is -200 to 1000ms, and we here add 300 ms before and after this because time-frequency edge artifacts (this is different than the first pass because we were being more conservative then) + %% === Baseline adjustment ==== + EEG = pop_rmbase( EEG, [-200/1000,0]); + + %% === Bad trial identification and removal ==== + for chan=1:29 + [EEG1 Indexes] = pop_eegthresh(EEG,1,chan,-150,150,0,1.2,0,0); + trials_removed(chan)=length(Indexes)/size(EEG.data,3); + end + %Save rejection information into a variable: electrodes with a trial rejection rate greater than 40% were tagged for removal + p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number + p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed + p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed + +end +%% === Save the detected bad channels ==== +save('chans_eeglab','p_chanreject'); + +%% === 2nd pass: process data ==== clear all; clc; -load('eeglab_removedchans500'); +load('chans_eeglab'); filenames = dir('*.vhdr') nb=500; for participant = 1:nb %Cycle through participants @@ -81,7 +95,7 @@ participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components %% === Read data ==== - EEG = pop_loadbv(pwd, filenames(participant).name); + EEG = pop_loadbv('../../data/', filenames(participant).name); %% === Add the channel location file ==== EEG=pop_chanedit(EEG, 'lookup','eeglab2022.1/functions/supportfiles/Standard-10-20-Cap81.ced'); @@ -122,9 +136,14 @@ 'TP9'}); end -% % %% === Interpolate bad channels ==== -% % EEG = pop_interp(EEG, rm_channels{participant}, 'spherical'); -% % All_rejChan(participant)=length(rm_channels{participant}); + %% === Interpolate bad channels ==== + + p_chans = []; %Clear past participants electrode indices + p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove + % save the number of channels removed/interpolated + All_rejChan(participant)=length(p_chans); + chans_to_remove = cellfun(@str2num,p_chans); %Clear past participants electrode labels + EEG = pop_interp(EEG, chans_to_remove, 'spherical'); %% === Re-referencing ==== EEG = pop_chanedit(EEG, 'seteeglab',{'1:31','TP10 TP9'}); @@ -141,7 +160,7 @@ EEG = pop_rmbase( EEG, [-200/1000,0]); %% === Bad trial identification and removal ==== - EEG = pop_eegthresh(EEG,1,[1:29] ,-50,50,0,1.2,0,1); + EEG = pop_eegthresh(EEG,1,[1:29] ,-100,100,0,1.2,0,1); %% === ERP ==== try @@ -164,21 +183,20 @@ %% === Save variables and csv files ==== -save('EEGLAB_All_ERP_withoutchans', 'All_ERP'); %Save ERP Data -save('EEGLAB_All_trials_withoutchans', 'All_trials'); %Save ERP Data -% save('EEGLAB_All_rejChan_withoutchans300', 'All_rejChan'); %Save ERP Data - -% All_ERP=All_ERP_eeglab(:,151:750,:,:); -% -% chanOfInterest=17; -% % channel of interest 17 is FCz -% %% RewP_Waveforms_AllPs -% tt1=squeeze(All_ERP(chanOfInterest,:,1,:)); -% tt2=squeeze(All_ERP(chanOfInterest,:,2,:)); -% % %% RewP_Waveforms -% csvwrite('eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -% csvwrite('eeglab_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. -% % %% RewP_Latency -% [~,peak_loc] = max(squeeze(All_ERP(chanOfInterest,226:276,1,:))-squeeze(All_ERP(chanOfInterest,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. -% peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds -% csvwrite('eeglab_RewP_Latency.csv',peak_loc'); %Export data +save('results/EEGLAB_All_ERP_samePipe', 'All_ERP'); %Save ERP Data +save('results/EEGLAB_All_trials_samePipe', 'All_trials'); %Save ERP Data +save('results/EEGLAB_All_rejChan_samePipe', 'All_rejChan'); %Save ERP Data + +All_ERP=All_ERP_eeglab(:,151:750,:,:); +chanOfInterest=17; +% channel of interest 17 is FCz +%% RewP_Waveforms_AllPs +win_erp=squeeze(All_ERP(chanOfInterest,:,1,:)); +loss_erp=squeeze(All_ERP(chanOfInterest,:,2,:)); +% %% RewP_Waveforms +csvwrite('eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +csvwrite('eeglab_RewP_Waveforms_AllPs.csv',[win_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +% %% RewP_Latency +[~,peak_loc] = max(squeeze(All_ERP(chanOfInterest,226:276,1,:))-squeeze(All_ERP(chanOfInterest,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. +peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds +csvwrite('eeglab_RewP_Latency.csv',peak_loc'); %Export data From 266bff267190c5ff86932c1686c790223c9ddf40 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 23:40:10 +0300 Subject: [PATCH 07/18] modified paths --- src/eeglab/eeglab_stats.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/eeglab/eeglab_stats.m b/src/eeglab/eeglab_stats.m index 5292a52..8c35611 100644 --- a/src/eeglab/eeglab_stats.m +++ b/src/eeglab/eeglab_stats.m @@ -4,9 +4,9 @@ nbpt=600; % % Section 1: Compute the statistical map(channels x time) showing the difference between FieldTrip and paper results in terms of conditional erp -load('results/Ref_All_ERP.mat') +load('../../results/Ref_All_ERP.mat') All_ERP1=All_ERP; -load('results/eeglab_All_ERP_samePipe.mat') +load('../../results/eeglab_All_ERP_samePipe.mat') All_ERP2=All_ERP; All_ERP1=All_ERP1(:,151:750,:,:); All_ERP2=All_ERP2(:,151:750,:,:); @@ -93,4 +93,4 @@ ylabel('Channels'); title([condition_name{condition} ' condition']); end -save('matp_eeglab','matp'); \ No newline at end of file +save('../../results/matp_eeglab','matp'); \ No newline at end of file From a82acbed9d4cf4ccbefdca79ce3b1231e50aedfc Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 23:40:51 +0300 Subject: [PATCH 08/18] changed path --- src/fieldtrip/ftpreprocessing.m | 83 ++++++++++++++++----------------- 1 file changed, 39 insertions(+), 44 deletions(-) diff --git a/src/fieldtrip/ftpreprocessing.m b/src/fieldtrip/ftpreprocessing.m index f687b79..8ca02ab 100644 --- a/src/fieldtrip/ftpreprocessing.m +++ b/src/fieldtrip/ftpreprocessing.m @@ -1,25 +1,22 @@ % % --- This is the preprocessing workflow reproduced by FieldTrip functions--- % % For contact: aya.kabbara7@gmail.com -%% === Regulate some parameters that will be used in interpolation ==== +cd('../../tools/fieldtrip'); -load('src/channelsTokeep.mat'); +%% Regulate some parameters that will be used in interpolation ==== % prepare neighborhood template for bad channels interpolation reduced_subjects=ss; cfg = []; -cfg.dataset = '/Raw Data Part 1/set/Subject001/set_001.set'; +cfg.dataset = '../../data/set/Subject001/set_001.set'; cfg.continuous = 'yes'; % force it to be continuous cfg.channel = ChannelsTokeep; data_prepare = ft_preprocessing(cfg); cfg=[]; cfg.method ='distance'; -cfg.neighbourdist=0.7; +cfg.neighbourdist = 0.7; [neighbours, cfg] = ft_prepare_neighbours(cfg, data_prepare) - -% % go to data folder -cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data -filenames = dir('*.vhdr') +filenames = dir('../../data/*.vhdr') nb=500; trials_loss=[]; trials_win=[]; @@ -28,17 +25,13 @@ for participant =1:500 %Cycle through participants -% Get participant name information + %% Get participant name information disp(['Participant: ', num2str(participant)]) %Display current participant being processed participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - RawFile = ['/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; + RawFile = ['../../data/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; SubjectName = ['participant_' participant_number{2}]; - -% % Check if the folder contains the required files -% if ~file_exist(RawFile) -% error(['The folder does not contain the folder from the file sample.']); -% end + %% === Read data ==== cfg = []; cfg.dataset = RawFile; @@ -55,7 +48,7 @@ %% === Remove the reference channels ==== cfg=[]; - if(length(find(reduced_subjects==participant))==0) + if(length(find(reduced_subjects==participant))==0) cfg.channel = [1:29]; else cfg.channel = [1:9 11:20 22:31]; @@ -88,14 +81,12 @@ end end - %% === Interpolate bad channels ==== cfg=[]; cfg.badchannel = bad_detected; cfg.neighbours = neighbours(index_bad_detected); try [data_interp] = ft_channelrepair(cfg, data_eeg_filtered) - %% === Segmentation into time-locked epochs ==== cfg = []; cfg.dataset= RawFile; @@ -126,69 +117,73 @@ cfg=[]; cfg.trl=cfg_loss.trl; cfg.continuous = 'no'; - cfg.artfctdef.threshold.range=200; -% cfg.artfctdef.threshold.min = -50; -% cfg.artfctdef.threshold.max = 50; + cfg.artfctdef.threshold.range=100; cfg.artfctdef.threshold.bpfilter = 'yes'; cfg.artfctdef.threshold.bpfreq = [0.1 30]; [cfg, artifact] = ft_artifact_threshold(cfg, data_loss_corrected) % % artifact rejection for loss condition try - data_loss_final = ft_rejectartifact(cfg, data_loss_corrected) + data_loss_final = ft_rejectartifact(cfg, data_loss_corrected) catch end % % win condition cfg=[]; cfg.trl=cfg_win.trl; cfg.continuous = 'no'; - cfg.artfctdef.threshold.range=100; -% cfg.artfctdef.threshold.min = -50; -% cfg.artfctdef.threshold.max = 50; + cfg.artfctdef.threshold.range=300; cfg.artfctdef.threshold.bpfilter = 'yes'; cfg.artfctdef.threshold.bpfreq = [0.1 30]; [cfg, artifact] = ft_artifact_threshold(cfg, data_win_corrected) + % % artifact rejection for win condition - try - data_win_final = ft_rejectartifact(cfg, data_win_corrected) - catch - end - %% === ERP calculation ==== - cfg=[]; try - [timelock] = ft_timelockanalysis(cfg, data_loss_final); - trials_loss(participant)=length(data_loss_final.trial); - All_ERP_ft(1,:,:,participant) = timelock.avg; + data_win_final = ft_rejectartifact(cfg, data_win_corrected) catch end - + %% === ERP calculation ==== + cfg=[]; try - [timelock] = ft_timelockanalysis(cfg, data_win_final); - trials_win(participant)=length(data_win_final.trial); - All_ERP_ft(2,:,:,participant) = timelock.avg; + [timelock] = ft_timelockanalysis(cfg, data_loss_final); + trials_loss(participant)=length(data_loss_final.trial); + All_ERP_ft(2,:,:,participant) = timelock.avg; + + [timelock] = ft_timelockanalysis(cfg, data_win_final); + trials_win(participant)=length(data_win_final.trial); + All_ERP_ft(1,:,:,participant) = timelock.avg; + + All_trials(1,participant) = trials_win(participant); %Store the number of trials + All_trials(2,participant) = trials_loss(participant); %Store the number of trials + All_trials(3,participant) = trials_loss(participant)+trials_win(participant); %Store the number of trials + catch end + catch end + end %% === Save variables and csv files ==== +All_ERP=All_ERP_ft; +save('../../results/ft_All_ERP', 'All_ERP'); %Save ERP Data +save('../../results/ft_All_trials', 'All_trials'); %Save ERP Data +save('../../results/ft_All_rejChan', 'All_rejChan'); %Save ERP Data - save('All_ERP_ft.mat','All_ERP_ft'); channelOfInterest=26; All_ERP_ft=All_ERP_ft(:,:,151:750,:); % %% RewP_Waveforms_AllPs -tt1=squeeze(All_ERP_ft(1,26,:,:)); -tt2=squeeze(All_ERP_ft(2,26,:,:)); +win_erp=squeeze(All_ERP_ft(1,26,:,:)); +loss_erp=squeeze(All_ERP_ft(2,26,:,:)); -csvwrite('ft_RewP_Waveforms_final22.csv',[(-200:2:998)',nanmean(squeeze(All_ERP_ft(1,26,:,:)),2),nanmean(squeeze(All_ERP_ft(2,26,:,:)),2),nanmean(squeeze(All_ERP_ft(1,26,:,:)),2)-nanmean(squeeze(All_ERP_ft(2,26,:,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +csvwrite('../../results/ft_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP_ft(1,26,:,:)),2),nanmean(squeeze(All_ERP_ft(2,26,:,:)),2),nanmean(squeeze(All_ERP_ft(1,26,:,:)),2)-nanmean(squeeze(All_ERP_ft(2,26,:,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -csvwrite('ft_RewP_Waveforms_AllPs_final22.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +csvwrite('../../results/ft_RewP_Waveforms_AllPs.csv',[win_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. % % %% RewP_Latency % [~,peak_loc] = max(squeeze(All_ERP_ft(1,26,226:276,:))-squeeze(All_ERP_ft(2,26,226:276,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds -csvwrite('ft_RewP_Latency.csv',peak_loc'); %Export data +csvwrite('../../results/ft_RewP_Latency.csv',peak_loc'); %Export data From 7683b9657c735b1f860db39b556d411b46e68250 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 23:41:11 +0300 Subject: [PATCH 09/18] paths changed --- src/fieldtrip/ft_stats.m | 243 ++++++++++----------------------------- 1 file changed, 58 insertions(+), 185 deletions(-) diff --git a/src/fieldtrip/ft_stats.m b/src/fieldtrip/ft_stats.m index fe3e847..02ca835 100644 --- a/src/fieldtrip/ft_stats.m +++ b/src/fieldtrip/ft_stats.m @@ -1,64 +1,51 @@ close all % clear - nbchan=29; nbpt=600; -P_Bonf=0.01/nbchan/nbpt; % % Section 1: Compute the statistical map(channels x time) showing the difference between FieldTrip and paper results in terms of conditional erp +load('../../results/Ref_All_ERP.mat') +All_ERP1=All_ERP; +load('../../results/ft_All_ERP_samePipe.mat') +All_ERP2=All_ERP; -StatDiff=zeros(nbchan,nbpt); - -load('All_ERP_ft.mat') -All_ERP2=All_ERP_ft; -load('All_ERP_ref.mat') - - -All_ERP=All_ERP(:,151:750,:,:); -tt1=squeeze(All_ERP2(1,26,:,:)); -tt2=squeeze(All_ERP2(2,26,:,:)); - -% % remove zeros and nan -idx1 = isnan(tt1) ; -[r1,c1]=find(tt1==0); -[r1,c3]=find(idx1); +All_ERP1=All_ERP1(:,151:750,:,:); +All_ERP2=All_ERP2(:,:,:,:); -idx2 = isnan(tt2) ; -[r2,c2]=find(tt2==0); -[r1,c4]=find(idx2); - -toberemoved=unique([unique(c1) ;unique(c2); unique(c3); unique(c4)]); -tKept1=[];tKept2=[]; -kept_ft=[]; -for su=1:500 - if(length(find(toberemoved==su))==0) -% % the subject su is present - kept_ft(end+1)=su; - tKept1(:,end+1)=tt1(:,su); - tKept2(:,end+1)=tt2(:,su); +condition_name={'Win','Loss'}; +matp=zeros(29,600,3); +for cond=1:3 + for ch=1:29 + ch + if(cond==3) + x1=squeeze(All_ERP1(ch,:,1,:))-squeeze(All_ERP1(ch,:,2,:)); + x2=squeeze(All_ERP2(2,ch,:,:))-squeeze(All_ERP2(1,ch,:,:)); + else + if(cond==1) + x1=squeeze(All_ERP1(ch,:,1,:)); + x2=squeeze(All_ERP2(2,ch,:,:)); + else + if(cond==2) + x1=squeeze(All_ERP1(ch,:,2,:)); + x2=squeeze(All_ERP2(1,ch,:,:)); + end + end + end + [clusters, p_values, t_sums, permutation_distribution ] = permutest(x1,x2,false,0.005,10000,true); + + for cl=1:length(clusters) + if(p_values(cl)<0.005) + significant_samples=clusters{cl}; + matp(ch,clusters{cl},cond)=1; + end + end end end -% % -All_ERP2=All_ERP2(:,:,151:750,kept_ft); -condition_name={'Win','Loss'}; - -for condition=1:2 - StatDiff=zeros(nbchan,nbpt); - - for channel=1:nbchan - for t=1:nbpt - Presults(channel,t)=ranksum(squeeze(All_ERP(channel,t,condition,:)),squeeze(All_ERP2(condition,channel,t,:))); - if(Presults(channel,t) Date: Tue, 11 Oct 2022 23:42:20 +0300 Subject: [PATCH 10/18] changed paths --- src/eeglab/eeglab_preprocessing.m | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/eeglab/eeglab_preprocessing.m b/src/eeglab/eeglab_preprocessing.m index 6554b02..11d0ccf 100644 --- a/src/eeglab/eeglab_preprocessing.m +++ b/src/eeglab/eeglab_preprocessing.m @@ -2,7 +2,7 @@ % % --- This is the preprocessing workflow reproduced by EEGLAB functions--- % % For contact: aya.kabbara7@gmail.com -cd('tools/eeglab2022.1'); +cd('../../tools/eeglab2022.1'); %% === 1st pass: Detect bad channels ==== filenames = dir('../../data/*.vhdr') @@ -183,9 +183,9 @@ %% === Save variables and csv files ==== -save('results/EEGLAB_All_ERP_samePipe', 'All_ERP'); %Save ERP Data -save('results/EEGLAB_All_trials_samePipe', 'All_trials'); %Save ERP Data -save('results/EEGLAB_All_rejChan_samePipe', 'All_rejChan'); %Save ERP Data +save('../../results/EEGLAB_All_ERP_samePipe', 'All_ERP'); %Save ERP Data +save('../../results/EEGLAB_All_trials_samePipe', 'All_trials'); %Save ERP Data +save('../../results/EEGLAB_All_rejChan_samePipe', 'All_rejChan'); %Save ERP Data All_ERP=All_ERP_eeglab(:,151:750,:,:); chanOfInterest=17; @@ -194,9 +194,9 @@ win_erp=squeeze(All_ERP(chanOfInterest,:,1,:)); loss_erp=squeeze(All_ERP(chanOfInterest,:,2,:)); % %% RewP_Waveforms -csvwrite('eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -csvwrite('eeglab_RewP_Waveforms_AllPs.csv',[win_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +csvwrite('../../results/eeglab_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2),nanmean(squeeze(All_ERP(chanOfInterest,:,1,:)),2)-nanmean(squeeze(All_ERP(chanOfInterest,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +csvwrite('../../results/eeglab_RewP_Waveforms_AllPs.csv',[win_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. % %% RewP_Latency [~,peak_loc] = max(squeeze(All_ERP(chanOfInterest,226:276,1,:))-squeeze(All_ERP(chanOfInterest,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds -csvwrite('eeglab_RewP_Latency.csv',peak_loc'); %Export data +csvwrite('../../results/eeglab_RewP_Latency.csv',peak_loc'); %Export data From 6a158052dabad6a93bba242fdbfffd08165b5677 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Tue, 11 Oct 2022 23:43:03 +0300 Subject: [PATCH 11/18] changed to permutation tests --- src/BST/bsPreprocessing.m | 147 ++++++++++++------------ src/BST/bs_stats.m | 232 ++++++++------------------------------ 2 files changed, 113 insertions(+), 266 deletions(-) diff --git a/src/BST/bsPreprocessing.m b/src/BST/bsPreprocessing.m index 68b4796..ea66650 100644 --- a/src/BST/bsPreprocessing.m +++ b/src/BST/bsPreprocessing.m @@ -1,7 +1,7 @@ % % --- This is the preprocessing workflow reproduced by Brainstorm functions--- % % For contact: aya.kabbara7@gmail.com - +cd('../../tools/brainstorm'); % ======= CREATE PROTOCOL ======= % The protocol name has to be a valid folder name (no spaces, no weird characters...) ProtocolName = 'Protocol_PreProc'; @@ -14,21 +14,23 @@ % Create new protocol gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); -cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data -filenames = dir('*.vhdr') +% cd('data/Raw Data Part 1'); %Find and change working folder to raw EEG data +filenames = dir('../../data/*.dat') nb=500; trials_1=0; trials_2=0; -load('/Raw Data Part 1/set/Subject001/channelsTokeep.mat'); -BS_db='Documents/brainstorm_db/'; %Changed this to the brainstorm_db directory +load('channelsTokeep.mat'); +BS_db='/Users/ayakabbara/Documents/brainstorm_db/'; %Changed this to the brainstorm_db directory for participant =1:nb %Cycle through participants + trials_1=0; + trials_2=0; % Get participant name information disp(['Participant: ', num2str(participant)]) %Display current participant being processed participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - RawFile = ['/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; + RawFile = ['/Users/ayakabbara/StageEEGpre/data/Raw Data Part 1/' filenames(participant).name]; SubjectName = ['participant_' participant_number{2}]; % Check if the folder contains the required files @@ -40,34 +42,36 @@ % Process: Create link to raw file sFiles = bst_process('CallProcess', 'process_import_data_raw', [], [], ... 'subjectname', SubjectName, ... - 'datafile', {RawFile,'EEG-EEGLAB'} , ... + 'datafile', {RawFile,'EEG-BRAINAMP'} , ... 'channelreplace', 0, ... 'channelalign', 0) % % Start a new report bst_report('Start', sFiles); - %% ===== EEG REFERENCE ===== - % Process: Re-reference EEG + %% ===== EEG REFERENCE ===== sFiles=bst_process('CallProcess', 'process_eegref', sFiles, [], ... 'eegref', 'TP9, TP10', ... 'sensortypes', 'EEG'); - %% ==== Bad channel identification ==== - sFiles=bst_process('CallProcess', 'process_import_data_time', sFiles, [], ... - 'subjectname', SubjectName, ... - 'timewindow', []); - - % detect flat channels only - sFilesInterp=bst_process('CallProcess', 'process_detectbad', sFiles, [], ... - 'timewindow', [], ... - 'eeg', [10,2000],'rejectmode',1); - - %% ==== Bad channel interpolation ==== - sFiles=bst_process('CallProcess', 'process_eeg_interpbad', [sFilesInterp], []); - + %% ==== Bad channels identification ==== + p_chans = []; %Clear past participants electrode indices + p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove + All_rejChan(participant)=length(p_chans); + try + chans_to_remove = cellfun(@str2num,p_chans); %Clear past participants electrode labels + a_t=load([BS_db ProtocolName '/data/' sFiles.ChannelFile]) + for ch=1:length(p_chans) + chans_labels{ch}=a_t.Channel(chans_to_remove(ch)).Name; + end + sFiles=bst_process('CallProcess', 'process_channel_setbad', sFiles, [], ... + 'subjectname', SubjectName, ... + 'sensortypes', chans_labels); + %% ==== Bad channel interpolation ==== + sFiles=bst_process('CallProcess', 'process_eeg_interpbad', [sFiles], []); + catch + end %% ==== Filtering ==== - % Process: Notch filter: 50Hz 100Hz 150Hz sFiles = bst_process('CallProcess', 'process_notch', sFiles, [], ... 'sensortypes', 'EEG', ... @@ -76,7 +80,7 @@ 'useold', 0, ... 'read_all', 0); - % % Process: Band-pass:0.1Hz-30Hz + % Process: Band-pass:0.1Hz-30Hz sFiles = bst_process('CallProcess', 'process_bandpass', sFiles, [], ... 'sensortypes', 'EEG', ... 'highpass', 0.1, ... @@ -90,102 +94,89 @@ %% ===== Epoching ===== markers = {'S110','S111'}; %Loss, win - % toremove={'AF3','AF4','AF7','AF8','C1','C2','C5','C6','CP3','CPz','F1','F2','F5','F6','FC3','FC4','FT10','FT7','FT8','FT9','Oz','P1','P2','P5','P6','PO3','PO4','PO7','PO8','TP7','TP8','CP4'}; - % sFiles=bst_process('CallProcess', 'process_channel_setbad', sFiles, [], ... - % 'sensortypes', toremove); - sFilesEpochs1 = bst_process('CallProcess', 'process_import_data_event', sFiles, [], ... + sFilesEpochs1 = bst_process('CallProcess', 'process_import_data_event', sFiles, [], ... 'subjectname', SubjectName, ... 'eventname', 'S110', ... 'timewindow', [], ... 'epochtime', [-0.5, 1.3], ... 'baseline', [-0.2, 0]); - sFilesEpochs2 = bst_process('CallProcess', 'process_import_data_event', sFiles, [], ... + sFilesEpochs2 = bst_process('CallProcess', 'process_import_data_event', sFiles, [], ... 'subjectname', SubjectName, ... 'eventname', 'S111', ... 'timewindow', [], ... 'epochtime', [-0.5, 1.3], ... 'baseline', [-0.2, 0]); - sFilesEpochs1 = bst_process('CallProcess', 'process_baseline', sFilesEpochs1, [],'baseline', [-0.2, 0]) - sFilesEpochs2 = bst_process('CallProcess', 'process_baseline', sFilesEpochs2, [],'baseline', [-0.2, 0]) + sFilesEpochs1 = bst_process('CallProcess', 'process_baseline', sFilesEpochs1, [],'baseline', [-0.2, 0]) + sFilesEpochs2 = bst_process('CallProcess', 'process_baseline', sFilesEpochs2, [],'baseline', [-0.2, 0]) %% ===== Bad Trials detection ===== - sFilesEpochs3=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs1], [], ... + sFilesEpochs3=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs1], [], ... 'timewindow', [], ... 'eeg', [0,100],'rejectmode',2); -% - sFilesEpochs4=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs2], [], ... - 'timewindow', [], ... - 'eeg', [0,100],'rejectmode',2); + sFilesEpochs4=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs2], [], ... + 'timewindow', [], ... + 'eeg', [0,100],'rejectmode',2); %% ===== ERP computation ===== - sFilesAvg = bst_process('CallProcess', 'process_average', [sFilesEpochs3, sFilesEpochs4], [], ... + sFilesAvg = bst_process('CallProcess', 'process_average', [sFilesEpochs3, sFilesEpochs4], [], ... 'avgtype', 5, ... % By trial groups (folder average) - 'avg_func', 1, ... % Arithmetic average: mean(x) - 'weighted', 0, ... - 'keepevents', 1); + 'avg_func', 1, ... % Arithmetic average: mean(x) + 'weighted', 0, ... + 'keepevents', 1); - sFilesAvg2 = bst_process('CallProcess', 'process_average', [sFilesEpochs32, sFilesEpochs42], [], ... - 'avgtype', 5, ... % By trial groups (folder average) - 'avg_func', 1, ... % Arithmetic average: mean(x) - 'weighted', 0, ... - 'keepevents', 1); - % try - BS_db='Documents/brainstorm_db/'; - dirr=[BS_db ProtocolName '/data/' sFilesAvg(1).FileName]; + dirr=[BS_db ProtocolName '/data/' sFilesAvg(1).FileName]; - FF=load(dirr); - if(size(FF,1)>32) - All_ERP_BS(1,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable - else - All_ERP_BS(1,:,:,participant) = FF.F([1:9 11:20 22:31],151:750); %Store all the ERP data into a single variable + FF=load(dirr); + if(size(FF,1)>32) + All_ERP_BS(1,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable + else + All_ERP_BS(1,:,:,participant) = FF.F([1:9 11:20 22:31],151:750); %Store all the ERP data into a single variable - end - trials_1=trials_1+sFilesAvg(1).iStudy; + end + trials_1=trials_1+sFilesAvg(1).iStudy; catch - % rem_part(end+1)=participant; end try - dirr=[BS_db ProtocolName '/data/' sFilesAvg(2).FileName]; - FF=load(dirr); - if(size(FF,1)>32) - All_ERP_BS(2,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable - else - All_ERP_BS(2,:,:,participant) = FF.F([1:9 11:20 22:31],151:750); %Store all the ERP data into a single variable - end - trials_2=trials_2+sFilesAvg(2).iStudy; + dirr=[BS_db ProtocolName '/data/' sFilesAvg(2).FileName]; + FF=load(dirr); + if(size(FF,1)>32) + All_ERP_BS(2,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable + else + All_ERP_BS(2,:,:,participant) = FF.F([1:9 11:20 22:31],151:750); %Store all the ERP data into a single variable + end + trials_2=trials_2+sFilesAvg(2).iStudy; catch % rem_part(end+1)=participant; end - - + All_trials(1,participant) = trials_1; %Store the number of trials + All_trials(2,participant) = trials_2; %Store the number of trials + All_trials(3,participant) = trials_1+trials_2; %Store the number of trials % delete subject for space issues % bst_process('CallProcess', 'process_delete','subjectname', SubjectName) end %% === Save variables and csv files ==== -save('All_ERP_BS','All_ERP_BS'); +save('../../results/bs_All_ERP', 'All_ERP_BS'); %Save ERP Data +save('../../results/bs_All_trials', 'All_trials'); %Save ERP Data +save('../../results/bs_All_rejChan', 'All_rejChan'); %Save ERP Data All_ERP=All_ERP(:,:,151:750,:).*1000000; % the unit in BS is in microVolts so it should be transfomed channelOfInterest=26; - -tt1=squeeze(All_ERP(1,channelOfInterest,:,:)); -tt2=squeeze(All_ERP(2,channelOfInterest,:,:)); - - - -csvwrite('bs_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(1,26,:,:)),2),nanmean(squeeze(All_ERP(2,26,:,:)),2),nanmean(squeeze(All_ERP(1,26,:,:)),2)-nanmean(squeeze(All_ERP(2,26,:,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -% %% RewP_Waveforms_AllPs -csvwrite('bs_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. -% %% RewP_Latency +gain_erp=squeeze(All_ERP(1,channelOfInterest,:,:)); +loss_erp=squeeze(All_ERP(2,channelOfInterest,:,:)); +csvwrite('../../results/bs_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(1,26,:,:)),2),nanmean(squeeze(All_ERP(2,26,:,:)),2),nanmean(squeeze(All_ERP(1,26,:,:)),2)-nanmean(squeeze(All_ERP(2,26,:,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +%% RewP_Waveforms_AllPs +csvwrite('../../results/bs_RewP_Waveforms_AllPs.csv',[gain_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +%% RewP_Latency [~,peak_loc] = max(squeeze(All_ERP(1,26,226:276,:))-squeeze(All_ERP(2,26,226:276,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds peak_loc(toberemoved)=[]; -csvwrite('bs_RewP_Latency.csv',peak_loc'); %Export data +csvwrite('../../results/bs_RewP_Latency.csv',peak_loc'); %Export data diff --git a/src/BST/bs_stats.m b/src/BST/bs_stats.m index ab647fa..7ead584 100644 --- a/src/BST/bs_stats.m +++ b/src/BST/bs_stats.m @@ -1,60 +1,40 @@ close all % clear - nbchan=29; nbpt=600; -P_Bonf=0.01/nbchan/nbpt; - -% % Section 1: Compute the statistical map(channels x time) showing the difference between Brainstorm and paper results in terms of conditional erp - -StatDiff=zeros(nbchan,nbpt); - -load('All_ERP_BS.mat') -All_ERP2=All_ERP; - -load('All_ERP_ref.mat') -All_ERP=All_ERP(:,151:750,:,:); -tt1=squeeze(All_ERP2(1,26,:,:)); -tt2=squeeze(All_ERP2(2,26,:,:)); -% % remove zeros and nan -idx1 = isnan(tt1) ; -[r1,c1]=find(tt1==0); -[r1,c3]=find(idx1); - -idx2 = isnan(tt2) ; -[r2,c2]=find(tt2==0); -[r1,c4]=find(idx2); - -toberemoved=unique([unique(c1) ;unique(c2); unique(c3); unique(c4)]); -tKept1=[];tKept2=[]; -kept_bs=[]; -for su=1:500 - if(length(find(toberemoved==su))==0) -% % the subject su is present - kept_bs(end+1)=su; - tKept1(:,end+1)=tt1(:,su); - tKept2(:,end+1)=tt2(:,su); - end -end - -All_ERP2=All_ERP2(:,:,:,kept_bs).*1000000; - - -%% === Map for each condition (gain, loss) ==== - -condition_name={'Win','Loss'}; - -for condition=1:2 - StatDiff=zeros(nbchan,nbpt); - - for channel=1:nbchan - for t=1:nbpt - Presults(channel,t)=ranksum(squeeze(All_ERP(channel,t,condition,:)),squeeze(All_ERP2(condition,channel,t,:))); - if(Presults(channel,t) Date: Tue, 11 Oct 2022 23:44:41 +0300 Subject: [PATCH 12/18] removed manual interventetion --- src/article/RewardsPreprocessing_withoutICA.m | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/article/RewardsPreprocessing_withoutICA.m b/src/article/RewardsPreprocessing_withoutICA.m index c83f3e1..c17ab85 100755 --- a/src/article/RewardsPreprocessing_withoutICA.m +++ b/src/article/RewardsPreprocessing_withoutICA.m @@ -39,7 +39,7 @@ [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) %Baseline Correction - [EEG] = doBaseline(EEG,[-200,0]); %Baseline correction in ms + [EEG] = doBaseline(EEG,[-200/1000,0]); %Baseline correction in ms %Artifact Rejection [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 uV/ms gradient criteria @@ -187,7 +187,6 @@ save('All_trials', 'All_trials'); %Save ERP Data save('All_rejChan', 'All_rejChan'); %Save ERP Data -% All_ERP=All_ERP(:,151:750,:,:); % % %% RewP_Waveforms % csvwrite('Ref_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(26,:,1,:)),2),nanmean(squeeze(All_ERP(26,:,2,:)),2),nanmean(squeeze(All_ERP(26,:,1,:)),2)-nanmean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. % % %% RewP_Waveforms_AllPs From 2692dc99a2b778d6bd22374352018428c530ee44 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 10:17:41 +0300 Subject: [PATCH 13/18] readme changes --- README.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index a7aaac0..a871a53 100644 --- a/README.md +++ b/README.md @@ -19,7 +19,9 @@ This study sheds light on how the software tool used to preprocess EEG signals i ## Abstract -As an active field of research, Electroencephalography (EEG) analysis workflow has increasing flexibility and complexity, with a great variety of methodological options and tools to be selected at each step. This high analytical flexibility can be problematic as it can yield variability in research outcomes. Therefore, growing attention has been recently paid to understand the potential of different methodological decisions to influence the reproducibility of results. In this paper, we aim to examine how sensitive the results of EEG analyses are to variations in software tools. We reanalyzed shared EEG data (N=500) previously used in a published task EEG study using three of the most commonly used software tools: EEGLAB, Brainstorm and FieldTrip. After reproducing the same original preprocessing workflow in each software, the resultant evoked-related potentials (ERPs) were qualitatively and quantitatively compared in order to examine the degree of consistency/discrepancy between softwares. Our findings show a good degree of convergence in terms of the general profile of ERP waveforms, peak latencies and effect size estimates related to specific signal features. However, considerable variability was also observed between software packages as reflected by the similarity values and observed statistical differences at particular channels and time instants. In conclusion, we believe that this study provides valuable clues to better understand the impact of the software tool on analysis results of EEG. +As an active field of research and with the development of state-of-the-art algorithms to analyze EEG datasets, the parametrization of Electroencephalography (EEG) analysis workflows has become increasingly flexible and complex, with a great variety of methodological options and tools to be selected at each step. This high analytical flexibility can be problematic as it can yield to variability in research outcomes. Therefore, growing attention has been recently paid to understand the potential impact of different methodological decisions on the reproducibility of results. +In this paper, we aim to examine how sensitive the results of EEG analyses are to variations in preprocessing with different software tools. We reanalyzed the shared EEG data (N=500) from (Williams et al., 2021) using three of the most commonly used open-source Matlab-based EEG software tools: EEGLAB, Brainstorm and FieldTrip. After reproducing the same original preprocessing workflow in each software, the resulting evoked-related potentials (ERPs) were qualitatively and quantitatively compared in order to examine the degree of consistency/discrepancy between software packages. Our findings show a good degree of convergence in terms of the general profile of ERP waveforms, peak latencies and effect size estimates related to specific signal features. However, considerable variability was also observed in the magnitude of the absolute voltage observed with each software package as reflected by the similarity values and observed statistical differences at particular channels and time instants. In conclusion, we believe that this study provides valuable clues to better understand the impact of the software tool on the analysis of EEG results. + ## Software implementation @@ -62,18 +64,16 @@ https://github.com/Neuro-Tools/MATLAB-EEG-icaTools. - FieldTrip script: This script uses FieldTrip functions available at https://www.fieldtriptoolbox.org/download/ +All the tools directories are downloaded in `tools`. + ## Reproducing the results ### Reproducing Figure 2A, Figure 3A using the reference script of (Williams et al, 2021) > We used the same preprocessing code provided by the original paper, by eliminating the step of manual ICA employed to detect the arttifactual components related to eye blinks -1. Add path to corresponding dependencies listed above using: -``` - addpath(genpath('dependencyDir')) -``` -2. Run the analysis with [RewardsPreprocessing_withoutICA.m](https://github.com/AyaKabbara/StageEEGpre/tree/main/src/article/RewardsPreprocessing_withoutICA.m). This script will store the ERPs of all participants in a .mat file. It will also generate the csv files to be used for plotting the figures and the calculation of the descriptive statistics. -3. Create figures 2A,3A and generate the related values in tables 1,2 using the following R script: [RewardProcessing_Plots_and_Statistics.R](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/RewardProcessing_Plots_and_Statistics.R). +1. Run the analysis with [RewardsPreprocessing_withoutICA.m](src/article/RewardsPreprocessing_withoutICA.m). This script will store the ERPs of all participants in a .mat file. It will also generate the csv files to be used for plotting the figures and the calculation of the descriptive statistics. +2. Create figures 2A,3A and generate the related values in tables 1,2 using the following R script: [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R). Note: The name of csv files should be changed accordingly. ### Reproducing Figure 2B, Figure 3B using EEGLAB From 3fddb65095df15376977ade4912ef58f8f8ed33d Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 10:39:48 +0300 Subject: [PATCH 14/18] paths fixed --- src/article/RewardsPreprocessing_withoutICA.m | 84 ++++++++----------- 1 file changed, 35 insertions(+), 49 deletions(-) diff --git a/src/article/RewardsPreprocessing_withoutICA.m b/src/article/RewardsPreprocessing_withoutICA.m index c17ab85..c60e342 100755 --- a/src/article/RewardsPreprocessing_withoutICA.m +++ b/src/article/RewardsPreprocessing_withoutICA.m @@ -1,11 +1,9 @@ % %% PREPROCESSING % %% Stage 1: Process data to determine noisy/faulty electrodes % %% Step 1.1: Pre-ICA - clear all; - dirdata='data/Raw Data Part 1'; - cd(dirdata); %Find and change working folder to raw EEG data - filenames = dir('*.vhdr'); - +cd('../../tools/reference'); +filenames = dir('../../data/*.vhdr') +dirdata='../../data/'; for participant = 1:500 %Cycle through participants %Get participant name information @@ -15,7 +13,7 @@ %Load Data EEG = []; %Clear past data - [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data + [EEG] = doLoadBVData([dirdata filenames(participant).name] ); %Load raw EEG data %Make it a 32 channels cap - reduces any participants with a 64 channel system if EEG.nbchan > 31 %Determine whether current participant is a 64 channel setup @@ -24,17 +22,13 @@ %Re-Reference load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path - % chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references [EEG] = doRereference(EEG,{'TP9','TP10'},{'ALL'},EEG.chanlocs); %Reference all channels to an averaged TP9/TP10 (mastoids) [EEG] = doRemoveChannels(EEG,{'TP9','TP10'},EEG.chanlocs); %Remove reference channels from the data [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate the electrode used as reference during recording (AFz) %Filter [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter data: Low cutoff 0.1, high cutoff 30, order 4, notch 60 - -% %ICA - % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal - + %Segment Data [EEG] = doSegmentData(EEG,{'S110','S111'},[-500 1498]); %Segment Data (S110 = Loss, S111 = Win) @@ -52,21 +46,16 @@ trials_removed = sum(AR_indx')/size(AR_indx,2); %Determine the percentage of segments removed for each electrode %Save rejection information into a variable: electrodes with a trial rejection rate greater than 40% were tagged for removal - p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number - p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed - p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed + p_chanreject{participant,1} = cellstr(participant_number{2}); %Insert participant number + p_chanreject{participant,2} = cellstr(num2str(0.4)); %insert lowest level within which less than ten electrodes have been removed + p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>0.4))); %Determine the indices of electrodes to be removed - %Save Output -% save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention end - %Save channels to reject information into a mat file +%Save channels to reject information into a mat file save('Chans_rejected500.mat','p_chanreject'); %% Stage 2: Process data for analysis - clear all; close all; clc; %First, clean the environment -% dirdata='data/Raw Data Part 1'; -% cd(dirdata); %Find and change working folder to raw EEG data - filenames = dir('*.vhdr'); %Compile list of all data +filenames = dir('../../data/*.vhdr') load('Chans_rejected500.mat'); %Load list of electrodes to remove @@ -107,13 +96,13 @@ % end %Remove faulty electrodes - p_chans = []; %Clear past participants electrode indices - p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove - % save the number of channels removed/interpolated - All_rejChan(participant)=length(p_chans); + p_chans = []; %Clear past participants electrode indices + p_chans = strsplit(p_chanreject{participant,3}{1},' '); %Extract the indices of electrodes to remove + % save the number of channels removed/interpolated + All_rejChan(participant)=length(p_chans); - chans_to_remove = []; %Clear past participants electrode labels - if ~isempty(p_chans{1}) %Determine whether electrodes need to be removed + chans_to_remove = []; %Clear past participants electrode labels + if ~isempty(p_chans{1}) %Determine whether electrodes need to be removed x = 1; %Begin index count for l = 1:size(p_chans,2) %Scroll through the electrodes if str2num(p_chans{1,l}) ~= 26 %Ensure that FCz is not removed @@ -122,9 +111,9 @@ end end [EEG] = doRemoveChannels(EEG,chans_to_remove,chanlocsMaster); %Remove electrodes - else %If there is no electrodes to remove + else %If there is no electrodes to remove [EEG] = doRemoveChannels(EEG,{},chanlocsMaster); %Still run function, but do not remove any electrodes - end + end % This manual removal has been removed from this code as well as in other software codes @@ -142,15 +131,12 @@ %Filter [EEG] = doFilter(EEG,0.1,30,4,60,EEG.srate); %Filter da ta: Low cutoff 0.1, high cutoff 30, order 4, notch 60 - - %ICA - % [EEG] = doICA(EEG,1); %Run ICA for use of eye blink removal - + %Topograpic Interpolation load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate electrodes which were previously removed - %Determine markers of interest + %Determine markers of interest markers = {'S110','S111'}; %Loss, win %Segment Data @@ -183,19 +169,19 @@ save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention end - save('Ref_All_ERP', 'All_ERP'); %Save ERP Data - save('All_trials', 'All_trials'); %Save ERP Data - save('All_rejChan', 'All_rejChan'); %Save ERP Data +save('../../results/Ref_All_ERP', 'All_ERP'); %Save ERP Data +save('../../results/Ref_All_trials', 'All_trials'); %Save ERP Data +save('../../results/Ref_All_rejChan', 'All_rejChan'); %Save ERP Data + +%% RewP_Waveforms +csvwrite('../../results/Ref_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(26,:,1,:)),2),nanmean(squeeze(All_ERP(26,:,2,:)),2),nanmean(squeeze(All_ERP(26,:,1,:)),2)-nanmean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. +%% RewP_Waveforms_AllPs +win_erp=squeeze(All_ERP(26,:,1,:)); +loss_erp=squeeze(All_ERP(26,:,2,:)); -% % %% RewP_Waveforms -% csvwrite('Ref_RewP_Waveforms.csv',[(-200:2:998)',nanmean(squeeze(All_ERP(26,:,1,:)),2),nanmean(squeeze(All_ERP(26,:,2,:)),2),nanmean(squeeze(All_ERP(26,:,1,:)),2)-nanmean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -% % %% RewP_Waveforms_AllPs -% tt1=squeeze(All_ERP(26,:,1,:)); -% tt2=squeeze(All_ERP(26,:,2,:)); -% -% csvwrite('Ref_RewP_Waveforms_AllPs.csv',[tt1,tt2]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. -% % %% RewP_Latency -% [~,peak_loc] = max(squeeze(All_ERP(26,226:276,1,:))-squeeze(All_ERP(26,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. -% peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds -% peak_loc(toberemoved)=[]; -% csvwrite('Ref_RewP_Latency.csv',peak_loc'); %Export data +csvwrite('../../results/Ref_RewP_Waveforms_AllPs.csv',[win_erp,loss_erp]'); %Export data. Conditions: Loss, Win. Electrode 26 is FCz. +%% RewP_Latency +[~,peak_loc] = max(squeeze(All_ERP(26,226:276,1,:))-squeeze(All_ERP(26,226:276,2,:))); %Determine where the peak amplitude is for each participant. Electrode 26 is FCz. +peak_loc = (((peak_loc+225)*2)-200)/1000; %Convert into seconds +peak_loc(toberemoved)=[]; +csvwrite('../../results/Ref_RewP_Latency.csv',peak_loc'); %Export data From 411c43f8724dbb9d6f031b8a8c825115045fb3b8 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 10:49:12 +0300 Subject: [PATCH 15/18] fixed errors --- src/graphiques/Fig4.py | 27 ++++++++++++------------ src/graphiques/fig4_preparation.m | 35 +++++++++++++++++++------------ 2 files changed, 36 insertions(+), 26 deletions(-) diff --git a/src/graphiques/Fig4.py b/src/graphiques/Fig4.py index a971296..b98a75d 100644 --- a/src/graphiques/Fig4.py +++ b/src/graphiques/Fig4.py @@ -10,8 +10,9 @@ colors_list = ['#505050', '#9DBCEA', '#F78A88', '#FFD384'] + # Figure 4 A -dataframe = pd.read_csv("gain_peaks.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("gain_peaks_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -27,10 +28,10 @@ plt.show() fig = aa.get_figure() -fig.savefig("gain-peaks.png") +fig.savefig("gainpeat100.png") # Figure 4 B -dataframe = pd.read_csv("lose_peaks.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("lose_peaks_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -46,10 +47,10 @@ plt.show() fig = aa.get_figure() -fig.savefig("lose_peaks.png") +fig.savefig("losepeaks_t100.png") # Figure 4 C -dataframe = pd.read_csv("mean_peaks.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("meanpeaks_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -65,10 +66,10 @@ plt.show() fig = aa.get_figure() -fig.savefig("mean_peaks.png") +fig.savefig("meanpeaks_t100.png") # Figure 4 D -dataframe = pd.read_csv("max_peaks.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("maxpeaks_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -84,10 +85,10 @@ plt.show() fig = aa.get_figure() -fig.savefig("max-peaks.png") +fig.savefig("maxpeaks_t100.png") # Figure 4 E -dataframe = pd.read_csv("basetopeaks_peaks.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("basetopeaks_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -103,10 +104,10 @@ plt.show() fig = aa.get_figure() -fig.savefig("basetopeaks_peaks.png") +fig.savefig("basetopeaks_t100.png") # Figure 4 F -dataframe = pd.read_csv("peak_loc.csv", error_bad_lines=False, encoding="ISO-8859-1") +dataframe = pd.read_csv("peak_loc_samepipe.csv", error_bad_lines=False, encoding="ISO-8859-1") print(dataframe.head()) print(dataframe.isnull().values.any()) @@ -117,10 +118,10 @@ sns.set_theme() aa=sns.violinplot(data=dataframe,palette=colors_list,inner="box") -plt.ylabel('Voltage(μV)', fontsize=16); +plt.ylabel('Time(ms)', fontsize=16); plt.tick_params(axis='both', which='major', labelsize=14) plt.show() fig = aa.get_figure() -fig.savefig("peak_loc.png") +fig.savefig("peak_loc_t100.png") diff --git a/src/graphiques/fig4_preparation.m b/src/graphiques/fig4_preparation.m index c0f89df..981ab0c 100644 --- a/src/graphiques/fig4_preparation.m +++ b/src/graphiques/fig4_preparation.m @@ -1,23 +1,32 @@ +cd('../../results'); % % reference results load('Ref_All_ERP.mat'); -[peak_loc_ref, mean_peaks_ref, gain_peaks_ref, lose_peaks_ref, peak_time_ref,base_time_ref,basetopeaks_peaks_ref]=extract_features(All_ERP); +All_ERP=All_ERP(:,151:750,:,:); +[peak_loc_ref, mean_peaks_ref, gain_peaks_ref, lose_peaks_ref, peak_time_ref,base_time_ref,basetopeaks_peaks_ref]=ext_features(All_ERP,26); % % bs results -load('All_ERP_BS.mat'); -[peak_loc_bs, mean_peaks_bs, gain_peaks_bs, lose_peaks_bs, peak_time_bs,base_time_bs,basetopeaks_peaks_bs]=extract_features(All_ERP_BS); +load('bs_All_ERP_samepipe.mat'); +All_ERP=All_ERP_BS; +All_ERP=All_ERP.*1000000; % the unit in BS is in microVolts so it should be transfomed +[peak_loc_bs, mean_peaks_bs, gain_peaks_bs, lose_peaks_bs, peak_time_bs,base_time_bs,basetopeaks_peaks_bs]=ext_features(All_ERP,26); % % eeglab results -load('All_ERP_eeglab.mat'); -[peak_loc_eeglab, mean_peaks_eeglab, gain_peaks_eeglab, lose_peaks_eeglab, peak_time_eeglab,base_time_eeglab,basetopeaks_peaks_eeglab]=extract_features(All_ERP_eeglab); +load('EEGLAB_All_ERP_samePipe.mat'); +All_ERP=All_ERP(:,151:750,:,:); +[peak_loc_eeglab, mean_peaks_eeglab, gain_peaks_eeglab, lose_peaks_eeglab, peak_time_eeglab,base_time_eeglab,basetopeaks_peaks_eeglab]=ext_features(All_ERP,17); % % ft results -load('All_ERP_ft.mat'); -[peak_loc_ft, mean_peaks_ft, gain_peaks_ft, lose_peaks_ft, peak_time_ft,base_time_ft,basetopeaks_peaks_ft]=extract_features(All_ERP_ft); +load('ft_All_ERP_samePipe.mat'); +All_ERP_ft=All_ERP(:,:,151:750,:); + All_ERP=[]; +All_ERP(1,:,:,:)=All_ERP_ft(2,:,:,:); +All_ERP(2,:,:,:)=All_ERP_ft(1,:,:,:); +[peak_loc_ft, mean_peaks_ft, gain_peaks_ft, lose_peaks_ft, peak_time_ft,base_time_ft,basetopeaks_peaks_ft]=ext_features(All_ERP,26); % % save results -csvwrite('peak_loc.csv',[peak_loc_ref peak_loc_eeglab peak_loc_bs peak_loc_ft]); -csvwrite('mean_peaks.csv',[mean_peaks_ref mean_peaks_eeglab mean_peaks_bs mean_peaks_ft]); -csvwrite('gain_peaks.csv',[gain_peaks_ref gain_peaks_eeglab gain_peaks_bs gain_peaks_ft]); -csvwrite('lose_peaks.csv',[lose_peaks_ref lose_peaks_eeglab lose_peaks_bs lose_peaks_ft]); -csvwrite('peak_time.csv',[peak_time_ref peak_time_eeglab peak_time_bs peak_time_ft]); -csvwrite('basetopeaks_peaks.csv',[basetopeaks_peak_ref basetopeaks_peak_eeglab basetopeaks_peak_bs basetopeaks_peak_ft]); \ No newline at end of file +csvwrite('peak_loc_samepipe.csv',[peak_loc_ref peak_loc_eeglab peak_loc_bs peak_loc_ft]); +csvwrite('mean_peaks_samepipe.csv',[mean_peaks_ref mean_peaks_eeglab mean_peaks_bs mean_peaks_ft]); +csvwrite('gain_peaks_samepipe.csv',[gain_peaks_ref gain_peaks_eeglab gain_peaks_bs gain_peaks_ft]); +csvwrite('lose_peaks_samepipe.csv',[lose_peaks_ref lose_peaks_eeglab lose_peaks_bs lose_peaks_ft]); +csvwrite('peak_time_samepipe.csv',[peak_time_ref peak_time_eeglab peak_time_bs peak_time_ft]); +csvwrite('basetopeaks_peaks_samepipe.csv',[basetopeaks_peak_ref basetopeaks_peak_eeglab basetopeaks_peak_bs basetopeaks_peak_ft]); \ No newline at end of file From 356b8f4c0eaaebd9b16ea7fa429e6d47d6574cbb Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 10:50:35 +0300 Subject: [PATCH 16/18] readme changes --- README.md | 37 +++++++++++-------------------------- 1 file changed, 11 insertions(+), 26 deletions(-) diff --git a/README.md b/README.md index a871a53..6cb1de3 100644 --- a/README.md +++ b/README.md @@ -78,47 +78,32 @@ Note: The name of csv files should be changed accordingly. ### Reproducing Figure 2B, Figure 3B using EEGLAB -1. Add path to corresponding dependencies listed above using: -``` - addpath(genpath('dependencyDir')) -``` - -2. Run the analysis [eeglab_preprocessing.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/eeglab/eeglab_preprocessing.m). -3. Create figures 2B,3B and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. +1. Run the analysis [eeglab_preprocessing.m](src/eeglab/eeglab_preprocessing.m). +2. Create figures 2B,3B and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. ### Reproducing Figure 2C,Figure 3C using Brainstorm -1. Add path to corresponding dependencies listed above using: -``` - addpath(genpath('dependencyDir')) -``` -2. Convert the dataset into EEGLAB set using the following functions: [toset.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/BST/toset.m) -3. Run the analysis with [bsPreprocessing.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/BST/bsPreprocessing.m). -4. Create figures 2C,3C and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. +1. Run the analysis with [bsPreprocessing.m](src/BST/bsPreprocessing.m). +2. Create figures 2C,3C and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. ### Reproducing Figure 2D,Figure 3D using FieldTrip -1. Add path to corresponding dependencies listed above using: -``` - addpath(genpath('dependencyDir')) -``` - -2. Run the analysis with [ftPreprocessing.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/BST/ftPreprocessing.m). -3. Create figures 2D,3D and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. +2. Run the analysis with [ftPreprocessing.m](src/fieldtrip/ftPreprocessing.m). +3. Create figures 2D,3D and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. ### Reproducing Figure 4 -1. Prepare the csv files containing the feature to be plotted using [fig4_preparation.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/fig4_preparation.m). This script calls the function [extract_features.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/extract_features.m). -2. Create the figure using [Fig4.py](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/Fig4.py) +1. Prepare the csv files containing the feature to be plotted using [fig4_preparation.m](src/graphiques/fig4_preparation.m). This script calls the function [extract_features.m](src/graphiques/extract_features.m). +2. Create the figure using [Fig4.py](src/graphiques/Fig4.py) ### Reproducing Figure 5 -1. Generate the similarity matrix (for both conditions) using the function [similarity_calc.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/similarity_calc.m) -2. Create the figure using [Fig5.py](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/graphiques/Fig5.py) +1. Generate the similarity matrix (for both conditions) using the function [similarity_calc.m](src/graphiques/similarity_calc.m) +2. Create the figure using [Fig5.py](src/graphiques/Fig5.py) ### Reproducing Figure 6 -Create Figure 6.A using [eeglab_stats.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/eeglab/eeglab_stats.m), Figure 6.B using [bs_stats.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/BST/bs_stats.m) and Figure 6.C using [ft_stats.m](https://github.com/AyaKabbara/StageEEGpre/blob/main/src/fieldtrip/ft_stats.m) +Create Figure 6.A using [eeglab_stats.m](src/eeglab/eeglab_stats.m), Figure 6.B using [bs_stats.m](src/BST/bs_stats.m) and Figure 6.C using [ft_stats.m](src/fieldtrip/ft_stats.m) From f6ccd5dc5cb3f8963ef6af0b63c916cd193f0518 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 11:01:49 +0300 Subject: [PATCH 17/18] channeltokeep --- README.md | 2 +- src/BST/channelsTokeep.mat | Bin 0 -> 218 bytes src/graphiques/Fig5.py | 23 +++++++++++++++++++++-- src/graphiques/similarity_calc.m | 12 +++++++----- 4 files changed, 29 insertions(+), 8 deletions(-) create mode 100644 src/BST/channelsTokeep.mat diff --git a/README.md b/README.md index 6cb1de3..c1b5dcb 100644 --- a/README.md +++ b/README.md @@ -89,7 +89,7 @@ Note: The name of csv files should be changed accordingly. ### Reproducing Figure 2D,Figure 3D using FieldTrip -2. Run the analysis with [ftPreprocessing.m](src/fieldtrip/ftPreprocessing.m). +2. Run the analysis with [ftPreprocessing.m](src/fieldtrip/ftpreprocessing.m). 3. Create figures 2D,3D and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. ### Reproducing Figure 4 diff --git a/src/BST/channelsTokeep.mat b/src/BST/channelsTokeep.mat new file mode 100644 index 0000000000000000000000000000000000000000..8d78fd54558197873af8f92e4b5f36f5085569bc GIT binary patch literal 218 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2c0*X0%nwjV*I2WZRmZYXAaTNzs@7#SEDDG&=7V1UunmmkOu0^*7}kCPJ;JQ$M96rMF) zQaHf%^fa@fvdCw~jm2(BadBCGg&uWfSy_zC%)%kPo3^Z2(!^uFdn!YJ=e&u_mQLHZ YQ9-|b>!LN2W=!eowPj|QQ^GY302%B($N&HU literal 0 HcmV?d00001 diff --git a/src/graphiques/Fig5.py b/src/graphiques/Fig5.py index 890ba38..a9bbcd7 100644 --- a/src/graphiques/Fig5.py +++ b/src/graphiques/Fig5.py @@ -1,5 +1,24 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 29 21:58:03 2022 + +@author: ayakabbara +""" +import scipy.io import seaborn as sns -import matplotlib.pyplot as plt % matplotlib inline +import matplotlib.pyplot as plt + +mat = scipy.io.loadmat('mat_dif.mat') +aa=mat.get("matt_corr_diff") +names_list=['Reference','EEEGLAB','Brainstorm','FieldTrip'] ax = sns.heatmap(aa,linewidths=.5,square=True,xticklabels=names_list,yticklabels=names_list) fig = ax.get_figure() -fig.savefig("corrstudy_matrix.png",dpi=300) +fig.savefig("corrstudy_diff.png",dpi=300) + +mat = scipy.io.loadmat('mat_cor.mat') +aa=mat.get("mat_cor") +names_list=['Reference','EEEGLAB','Brainstorm','FieldTrip'] +ax = sns.heatmap(aa,linewidths=.5,square=True,xticklabels=names_list,yticklabels=names_list) +fig = ax.get_figure() +fig.savefig("corrstudy_gainloss.png",dpi=300) diff --git a/src/graphiques/similarity_calc.m b/src/graphiques/similarity_calc.m index 9ddd58e..8723351 100644 --- a/src/graphiques/similarity_calc.m +++ b/src/graphiques/similarity_calc.m @@ -1,24 +1,26 @@ function [loss_corr, gain_corr, diff_corr]=similarity_calc(All_ERP11,All_ERP12) nbchan=29; -nbsubjects=498; +nbsubjects=500; gain_cor=[]; loss_cor=[]; diff_cor=[]; +count=0; for subject=1:nbsubjects tt1=squeeze(All_ERP12(1,26,:,subject)); tt2=squeeze(All_ERP11(26,:,1,subject)); if(sum(tt1)*sum(tt2)>0) +count=count+1; for chan=1:nbchan cc=corrcoef(squeeze(All_ERP2(1,chan,:,subject)),squeeze(All_ERP1(chan,:,1,subject))); - gain_cor(chan,count+1)=cc(1,2); + gain_cor(chan,count)=cc(1,2); cc=corrcoef(squeeze(All_ERP2(2,chan,:,subject)),squeeze(All_ERP1(chan,:,2,subject))); - loss_cor(chan,count+1)=cc(1,2); + loss_cor(chan,count)=cc(1,2); cc=corrcoef(squeeze(All_ERP2(1,chan,:,subject))-squeeze(All_ERP2(2,chan,:,subject)),squeeze(All_ERP1(chan,:,1,subject))-squeeze(All_ERP1(chan,:,2,subject))); - diff_cor(chan,count+1)=cc(1,2); + diff_cor(chan,count)=cc(1,2); end end @@ -30,4 +32,4 @@ loss_corr = loss_corr(~isnan(loss_corr)); diff_corr = diff_corr(~isnan(diff_corr)); -gain_corr = gain_corr(~isnan(gain_corr)); \ No newline at end of file +gain_corr = gain_corr(~isnan(gain_corr)); From d35cf3f5242bdaeb26d8171412015a90cd27fa73 Mon Sep 17 00:00:00 2001 From: AyaKabbara Date: Wed, 12 Oct 2022 11:02:29 +0300 Subject: [PATCH 18/18] link ftpreprocessing.m --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c1b5dcb..e0af603 100644 --- a/README.md +++ b/README.md @@ -89,7 +89,7 @@ Note: The name of csv files should be changed accordingly. ### Reproducing Figure 2D,Figure 3D using FieldTrip -2. Run the analysis with [ftPreprocessing.m](src/fieldtrip/ftpreprocessing.m). +2. Run the analysis with [ftpreprocessing.m](src/fieldtrip/ftpreprocessing.m). 3. Create figures 2D,3D and generate the related values in tables 1,2 using the R script [RewardProcessing_Plots_and_Statistics.R](src/graphiques/RewardProcessing_Plots_and_Statistics.R).Note: the name of csv files should be changed accordingly. ### Reproducing Figure 4