diff --git a/README.md b/README.md index a7aaac0..e0af603 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,63 +64,46 @@ 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 -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) 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/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) 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 + load('ChanlocsMaster.mat'); %Load channel location file, please make sure the location of this file is in your path [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 + + %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 + %Baseline Correction + [EEG] = doBaseline(EEG,[-200/1000,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 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 channels to reject information into a mat file - save('Chans_rejected_auto500','p_chanreject'); + 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 - filenames = dir('*.vhdr'); %Compile list of all data +filenames = dir('../../data/*.vhdr') +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,101 +76,64 @@ 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 - %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 - - %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 - - %Save output - save(participant_varname,'EEG'); %Save the current output so that the lengthy ICA process can run without user intervention - catch + %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 + % 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 + 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 - 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'; + % 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 - 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 - + %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 + %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 @@ -276,44 +152,36 @@ %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'; - - 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 + 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; - 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 + +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,:)); + +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 diff --git a/src/eeglab/eeglab_preprocessing.m b/src/eeglab/eeglab_preprocessing.m index 644a233..11d0ccf 100644 --- a/src/eeglab/eeglab_preprocessing.m +++ b/src/eeglab/eeglab_preprocessing.m @@ -2,77 +2,90 @@ % % --- 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 - +cd('../../tools/eeglab2022.1'); %% === 1st pass: Detect bad channels ==== -filenames = dir('*.vhdr') +filenames = dir('../../data/*.vhdr') nb=500; for participant = 1:nb %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 - - 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'); + % 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 - + % 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 catch disp('Problem'); rm_channels{participant}=''; - continue; 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 -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'); +save('chans_eeglab','p_chanreject'); %% === 2nd pass: process data ==== - clear all; clc; -load('eeglab_removedchans500'); -cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data +load('chans_eeglab'); filenames = dir('*.vhdr') nb=500; for participant = 1:nb %Cycle through participants @@ -82,57 +95,60 @@ 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('../../data/', 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 ==== + + 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'}); + 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 +157,46 @@ [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] ,-100,100,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 +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 -tt1=squeeze(All_ERP(chanOfInterest,:,1,:)); -tt2=squeeze(All_ERP(chanOfInterest,:,2,:)); +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',[tt1,tt2]'); %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 diff --git a/src/eeglab/eeglab_stats.m b/src/eeglab/eeglab_stats.m index 41977f9..8c35611 100644 --- a/src/eeglab/eeglab_stats.m +++ b/src/eeglab/eeglab_stats.m @@ -1,98 +1,44 @@ 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); + +% % Section 2: Compute the cluster-based permutation maps for each +% condition- gain, loss, difference + +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)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));