diff --git a/README.md b/README.md index c3af393..7c4d8c2 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,28 @@ +# Successful reproduction of a large EEG study across software packages + +by +Aya Kabbara +Nina Forde +Camille Maumet +Mahmoud Hassan + +> This is a guide for researchers to reproduce the results and figures published in the paper + +This paper has been submitted for publication in NeuroImage Reports. + +> This study sheds light on how the software tool used to preprocess EEG signals impacts the analysis results and conclusions. EEGLAB, Brainstorm and FieldTrip were used to reproduce the same preprocessing pipeline of a published EEG study performed on 500 participants. + +![](figures/figure1_preprocess.001.jpeg) + +*The full pipeline of the study. (a) EEGs were acquired from 500 participants performing a simple gambling task of six blocks composed of 20 trials. (b) The same dataset was then preprocessed using the different software tools: Reference (using the code published by the original paper (Williams et al. 2021)), EEGLAB, Brainstorm and FieldTrip. The preprocessing steps to be performed in each tool includes: the reduction to 32 electrodes, the re-referencing, the automatic detection of bad electrodes, the band-pass filtering, the interpolation of bad channels, the segmentation into time-locked epochs and the removal of artifactual trials. (c) The preprocessed signals derived from the four preprocessing codes are used to validate and reproduce the reference statistics and hypotheses. A quantitative comparison between the resulting signals is also conducted in terms of signal features* + + +## 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. + + + # StageEEGpre dependencies: diff --git a/figures/figure1_preprocess.001.jpeg b/figures/figure1_preprocess.001.jpeg new file mode 100644 index 0000000..adcba7c Binary files /dev/null and b/figures/figure1_preprocess.001.jpeg differ diff --git a/src/BST/bsPreprocessing.m b/src/BST/bsPreprocessing.m index e9cc039..d91a736 100644 --- a/src/BST/bsPreprocessing.m +++ b/src/BST/bsPreprocessing.m @@ -1,6 +1,10 @@ +% % --- This is the preprocessing workflow reproduced by Brainstorm functions--- +% % For contact: aya.kabbara7@gmail.com + + % ======= CREATE PROTOCOL ======= % The protocol name has to be a valid folder name (no spaces, no weird characters...) -ProtocolName = 'Protocol_PreProc1'; +ProtocolName = 'Protocol_PreProc'; % Start brainstorm without the GUI if ~brainstorm('status') brainstorm nogui @@ -10,20 +14,21 @@ % Create new protocol gui_brainstorm('CreateProtocol', ProtocolName, 0, 0); -cd('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1'); %Find and change working folder to raw EEG data +cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data filenames = dir('*.vhdr') nb=500; trials_1=0; trials_2=0; -load('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/Subject001/channelsTokeep.mat'); +load('/Raw Data Part 1/set/Subject001/channelsTokeep.mat'); +BS_db='Documents/brainstorm_db/'; %Changed this to the brainstorm_db directory -for participant =421:500 %Cycle through participants +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 - RawFile = ['/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; + RawFile = ['/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; SubjectName = ['participant_' participant_number{2}]; % Check if the folder contains the required files @@ -39,12 +44,6 @@ 'channelreplace', 0, ... 'channelalign', 0) - - - % % Input files - % sFiles = {... - % 'Subject01/@rawNDARAA075AMK/data_0raw_NDARAA075AMK.mat'}; - % % Start a new report bst_report('Start', sFiles); @@ -54,7 +53,7 @@ 'eegref', 'TP9, TP10', ... 'sensortypes', 'EEG'); - % % % ==== bad channel identification === + % % % ==== Bad channel identification === sFiles=bst_process('CallProcess', 'process_import_data_time', sFiles, [], ... 'subjectname', SubjectName, ... 'timewindow', []); @@ -64,10 +63,11 @@ 'timewindow', [], ... 'eeg', [10,2000],'rejectmode',1); - % interpolate the flat channels + % % % ==== Bad channel interpolation === sFiles=bst_process('CallProcess', 'process_eeg_interpbad', [sFilesInterp], []); - % % ====filtering==== + % % ====Filtering==== + % Process: Notch filter: 50Hz 100Hz 150Hz sFiles = bst_process('CallProcess', 'process_notch', sFiles, [], ... 'sensortypes', 'EEG', ... @@ -85,7 +85,7 @@ 'attenuation', 'strict', ... % 60dB 'ver', '2019', ... % 2019 'mirror', 0, ... - 'read_all', 0); + 'read_all', 0); % % ===== Epoching ===== @@ -99,54 +99,44 @@ 'eventname', 'S110', ... 'timewindow', [], ... 'epochtime', [-0.5, 1.3], ... - 'baseline', [-0.2, 0]); + 'baseline', [-0.2, 0]); sFilesEpochs2 = bst_process('CallProcess', 'process_import_data_event', sFiles, [], ... 'subjectname', SubjectName, ... 'eventname', 'S111', ... 'timewindow', [], ... 'epochtime', [-0.5, 1.3], ... - 'baseline', [-0.2, 0]); + '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]) - % % ===== trials detection ===== + % % ===== Bad Trials detection ===== sFilesEpochs3=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs1], [], ... 'timewindow', [], ... 'eeg', [0,100],'rejectmode',2); -% % for the second threshold - - sFilesEpochs32=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs1], [], ... - 'timewindow', [], ... - 'eeg', [0,200],'rejectmode',2); - - + +% sFilesEpochs4=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs2], [], ... 'timewindow', [], ... - 'eeg', [0,100],'rejectmode',2); - % % for the second threshold - - sFilesEpochs42=bst_process('CallProcess', 'process_detectbad', [sFilesEpochs2], [], ... - 'timewindow', [], ... - 'eeg', [0,200],'rejectmode',2); - - % % ===== erp computation ===== - % + 'eeg', [0,100],'rejectmode',2); + + % % ===== ERP computation ===== + 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); -% % for the second threshold + 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='/Users/ayakabbara/Documents/brainstorm_db/'; + BS_db='Documents/brainstorm_db/'; dirr=[BS_db ProtocolName '/data/' sFilesAvg(1).FileName]; FF=load(dirr); @@ -160,25 +150,7 @@ catch % rem_part(end+1)=participant; end - -% % second threshold -try - BS_db='/Users/ayakabbara/Documents/brainstorm_db/'; - dirr=[BS_db ProtocolName '/data/' sFilesAvg2(1).FileName]; - - FF=load(dirr); - if(size(FF,1)>32) - All_ERP_BS2(1,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable - else - All_ERP_BS2(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; - catch - % rem_part(end+1)=participant; -end - - BS_db='/Users/ayakabbara/Documents/brainstorm_db/'; + try dirr=[BS_db ProtocolName '/data/' sFilesAvg(2).FileName]; FF=load(dirr); @@ -193,26 +165,26 @@ % rem_part(end+1)=participant; end -% % second threshold -BS_db='/Users/ayakabbara/Documents/brainstorm_db/'; - try - dirr=[BS_db ProtocolName '/data/' sFilesAvg2(2).FileName]; - FF=load(dirr); - if(size(FF,1)>32) - All_ERP_BS2(2,:,:,participant) = FF.F(ChannelsTokeep(1:29),151:750); %Store all the ERP data into a single variable - else - All_ERP_BS2(2,:,:,participant) = FF.F([1:9 11:20 22:31],151:750); %Store all the ERP data into a single variable - end - catch - % rem_part(end+1)=participant; - end % % ===delete subject for space issues % bst_process('CallProcess', 'process_delete','subjectname', SubjectName) end -% save('first500AllERP_BS','All_ERP_BS'); -save('first500AllERP_BS2','All_ERP_BS2'); +save('AllERP_BS','All_ERP_BS'); + +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,:,:)); + + -% save('fourT1','trials_1'); -% save('fourT2','trials_2'); \ No newline at end of file +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 +[~,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 diff --git a/src/BST/toset.m b/src/BST/toset.m index ce09c2e..cd069af 100644 --- a/src/BST/toset.m +++ b/src/BST/toset.m @@ -1,7 +1,7 @@ -% EEGLAB history file generated on the 04-Jan-2022 -% ------------------------------------------------ -% -cd('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1'); %Find and change working folder to raw EEG data +% % code to save .set file for each vhdr file +% % Contact: aya.kabbara7@gmail.com + +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 @@ -10,8 +10,8 @@ disp(['Participant: ', num2str(participant)]) %Display current participant being processed participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components participant_varname = ['set_',participant_number{2}]; %Create new file name - EEG = pop_loadbv('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/', filenames(participant).name); + EEG = pop_loadbv('/Raw Data Part 1/', filenames(participant).name); %ajouter chanlocs - EEG= pop_chanedit(EEG, 'lookup','/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/eeglab2021.1/functions/supportfiles/Standard-10-20-Cap81.ced'); - save(['/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/' participant_varname '.set'],'EEG'); -end \ No newline at end of file + EEG= pop_chanedit(EEG, 'lookup','/eeglab2021.1/functions/supportfiles/Standard-10-20-Cap81.ced'); + save(['/Raw Data Part 1/set/' participant_varname '.set'],'EEG'); +end diff --git a/src/article/RewardsPreprocessing_withoutICA.m b/src/article/RewardsPreprocessing_withoutICA.m index c6f0740..35263b4 100755 --- a/src/article/RewardsPreprocessing_withoutICA.m +++ b/src/article/RewardsPreprocessing_withoutICA.m @@ -1,330 +1,319 @@ - % %% PREPROCESSING -% %% Stage 1: Process data to determine noisy/faulty electrodes -% %% Step 1.1: Pre-ICA -clear all; -dirdata='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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 + % %% 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 - 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='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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 + %Load Data + EEG = []; %Clear past 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 + [EEG] = doLoadBVData(filenames(participant).name ); %Load raw EEG data - %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'); -%% Stage 2: Process data for analysis -%% Step 2.1: Pre-ICA -clear all; close all; clc; %First, clean the environment -dirdata='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1'; - -cd(dirdata); %Find and change working folder to raw EEG data -filenames = dir('*.vhdr'); %Compile list of all data - -%Load list of electrodes to remove -for participant = 1:length(filenames) %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 - - %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('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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 - [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 - end -end -%% Step 2.2: Post-ICA -clear all; close all; clc; %First, clean the environment -dirdata='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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 + %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 - %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='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/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 - markers = {'S110','S111'}; %Loss, win - - %Segment Data - [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 Correction - [EEG] = doBaseline(EEG,[-200/1000,0]); %Baseline correction in ms - - %Artifact Rejection - [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 microvolt/ms gradient criteria - [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 microvolt max-min criteria - [EEG] = doRemoveEpochs(EEG,EEG.artifactPresent,0); %Remove segments that violated artifact rejection criteria - - %ERP - try - [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses - 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='/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1'; - -cd(dirdata); %Find and change working folder to raw EEG data + %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_S2Final*'); %Compile list of all data + 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 -for participant = 1:length(filenames) %Cycle through participants - disp(['Participant: ', num2str(participant)]); %Display current participant being processed + %% Step 1.3: Determining faulty electrodes + % clear all; close all; clc; %First, clean the environment + dirdata='/Raw Data Part 1'; - %Load Data - EEG = []; %Clear past data - load(filenames(participant).name); %Load participant output + 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); - All_ERP(:,:,:,participant) = EEG.ERP.data; %Store all the ERP data into a single variable -end - -save('Ref_All_ERP', 'All_ERP'); %Save ERP Data + %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'); + %% 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 + + %Load list of electrodes to remove + for participant = 1:length(filenames) %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 + + %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 + 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 + [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 + 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 + [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate electrodes which were previously removed -All_ERP=All_ERP(:,151:750,:,:); + %Determine markers of interest + markers = {'S110','S111'}; %Loss, win + + %Segment Data + [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 Correction + [EEG] = doBaseline(EEG,[-200/1000,0]); %Baseline correction in ms + + %Artifact Rejection + [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 microvolt/ms gradient criteria + [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 microvolt max-min criteria + [EEG] = doRemoveEpochs(EEG,EEG.artifactPresent,0); %Remove segments that violated artifact rejection criteria + + %ERP + try + [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses + 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 + + 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 -% %% 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,:)); + save('Ref_All_ERP', 'All_ERP'); %Save ERP Data -idx1 = isnan(tt1) ; -[r1,c1]=find(idx1); -tt1(:,unique(c1))=[]; -tt2(:,unique(c1))=[]; + All_ERP=All_ERP(:,151:750,:,:); -idx2 = isnan(tt2) ; -[r2,c2]=find(idx2); -tt2(:,unique(c2))=[]; -tt1(:,unique(c2))=[]; -toberemoved=unique([unique(c1) ;unique(c2)]); + % %% 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('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 diff --git a/src/article/icaeeglab.m b/src/article/icaeeglab.m deleted file mode 100644 index 7f7f9ef..0000000 --- a/src/article/icaeeglab.m +++ /dev/null @@ -1,398 +0,0 @@ -%% PREPROCESSING -%% Stage 1: Process data to determine noisy/faulty electrodes -%% Step 1.1: Pre-ICA -clear all; -restoredefaultpath; -clear all; -% close all; -clc; %First, clean the environment - - -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -addpath(genpath('/home/StageEEGpre/src')) - -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-fileIO-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-icaTools-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-preProcessing-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-timeFrequencyAnalysis-master')) - - -% addpath('/home/StageEEGpre/dependencies/eeglab_current/eeglab2021.1') - -addpath(genpath('/home/dependencies/eeglab_current/eeglab2021.1')) - - -cd(dossier10); %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 - %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 - [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 - - %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 -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %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 - 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 -% - [EEG]= pop_clean_rawdata(EEG, 'FlatlineCriterion',5,'ChannelCriterion',0.8,'LineNoiseCriterion',4,'Highpass','off','BurstCriterion',20,'WindowCriterion',0.25,'BurstRejection','on','Distance','Euclidian','WindowCriterionTolerances',[-Inf 7] ); - - - - %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 - - save(participant_varname,'EEG'); %Save the current output -end -%% Step 1.3: Determining faulty electrodes -clear all; close all; clc; %First, clean the environment -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %Find and change working folder to saved data from last for loop -filenames = dir('RewardProcessing_S1Final_*'); %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 - - %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(numbers(rejection_level))); %insert lowest level within which less than ten electrodes have been removed - p_chanreject{participant,3} = cellstr(num2str(find(trials_removed>numbers(rejection_level)==1))); %Determine the indices of electrodes to be removed -end - -%Save channels to reject information into a mat file -save('Chans_rejected_auto','p_chanreject'); -%% Stage 2: Process data for analysis -%% Step 2.1: Pre-ICA -clear all; close all; clc; %First, clean the environment -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %Find and change working folder to raw EEG data -filenames = dir('*.vhdr'); %Compile list of all data - -%Load list of electrodes to remove -load('Chans_rejected_auto.mat'); - -for participant = 1:length(filenames) %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 - - %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 - chanlocsMaster([10,21]) = []; %Remove channels TP9 and TP10 as they will be the references - - %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 - [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 - -end -%% Step 2.2: Post-ICA -clear all; close all; clc; %First, clean the environment -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %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 -% - [EEG]= pop_clean_rawdata(EEG, 'FlatlineCriterion',5,'ChannelCriterion',0.8,'LineNoiseCriterion',4,'Highpass','off','BurstCriterion',20,'WindowCriterion',0.25,'BurstRejection','on','Distance','Euclidian','WindowCriterionTolerances',[-Inf 7] ); - - save(participant_varname,'EEG'); %Save the current output -end -%% Step 2.3: Final -clear all; close all; clc; %First, clean the environment -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %Find and change working folder to saved data from last for loop -filenames = dir('RewardProcessing_S2PostInvICA_*'); %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{3}]; %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 - [EEG] = doInterpolate(EEG,chanlocsMaster,'spherical'); %Interpolate electrodes which were previously removed - - %Determine markers of interest - markers = {'S110','S111'}; %Loss, win - - %Segment Data - [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 Correction - [EEG] = doBaseline(EEG,[-200/1000,0]); %Baseline correction in ms - - %Artifact Rejection - [EEG] = doArtifactRejection(EEG,'Gradient',10); %Use a 10 microvolt/ms gradient criteria - [EEG] = doArtifactRejection(EEG,'Difference',100); %Use a 100 microvolt max-min criteria - [EEG] = doRemoveEpochs(EEG,EEG.artifactPresent,0); %Remove segments that violated artifact rejection criteria - - %Save backup data because different window lengths for each processing type - EEG.backupdata = EEG.data; - - %Wavelet - [EEG.WAV] = doWAV(EEG,markers,[],1,30,30,6); %Conduct Time-Frequency Analyses: no baseline, ranging from frequencies 1 to 30 in 30 linear steps, using a Morlet parameter of 6 - EEG.WAV.eegdata = EEG.data; %Copy WAV data - - %Reduce Data Length for FFT - EEG.data = EEG.backupdata(:,1:750,:); %-500 to 1000ms - - %FFT - [EEG.FFT] = doFFT(EEG,markers); %Conduct FFT Analyses - - save(participant_varname,'EEG'); %Save the current output - - %Reduce Data Length for ERP - EEG.data = EEG.backupdata(:,151:750,:); %-200 to 1000ms - - %ERP - [EEG.ERP] = doERP(EEG,markers,0); %Conduct ERP Analyses - - 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 -dossier10='/home/StageEEGpre/data/Rawdata10icaauto/' - -cd(dossier10); %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 - - %FFT frequency resolution is 0.67, so must extract proper frequencies - All_FFT(:,1,:,participant) = mean(EEG.FFT.data(:,1:2,:),2); %Store all the FFT data into a single variable (1 Hz) - All_FFT(:,2,:,participant) = EEG.FFT.data(:,3,:); %Store all the FFT data into a single variable (2 Hz) - All_FFT(:,3,:,participant) = mean(EEG.FFT.data(:,4:5,:),2); %Store all the FFT data into a single variable (3 Hz) - All_FFT(:,4,:,participant) = EEG.FFT.data(:,6,:); %Store all the FFT data into a single variable (4 Hz) - All_FFT(:,5,:,participant) = mean(EEG.FFT.data(:,7:8,:),2); %Store all the FFT data into a single variable (5 Hz) - All_FFT(:,6,:,participant) = EEG.FFT.data(:,9,:); %Store all the FFT data into a single variable (6 Hz) - All_FFT(:,7,:,participant) = mean(EEG.FFT.data(:,10:11,:),2); %Store all the FFT data into a single variable (7 Hz) - All_FFT(:,8,:,participant) = EEG.FFT.data(:,12,:); %Store all the FFT data into a single variable (8 Hz) - All_FFT(:,9,:,participant) = mean(EEG.FFT.data(:,13:14,:),2); %Store all the FFT data into a single variable (9 Hz) - All_FFT(:,10,:,participant) = EEG.FFT.data(:,15,:); %Store all the FFT data into a single variable (10 Hz) - - All_WAV(:,:,:,:,participant) = EEG.WAV.data; %%Store all the WAV data into a single variable -end - -save('All_ERP', 'All_ERP'); %Save ERP Data -save('All_FFT', 'All_FFT'); %Save FFT Data -%Unfortunately, the WAV file is unable to save as it is too large, thus when using the All_WAV variable, this section must be run first -%% Load Data -load('All_ERP.mat'); %Load saved ERP data -load('All_FFT.mat'); %Load saved FFT data -%% RewP_Waveforms -csvwrite('RewP_Waveforms.csv',[(-200:2:998)',mean(squeeze(All_ERP(26,:,1,:)),2),mean(squeeze(All_ERP(26,:,2,:)),2),mean(squeeze(All_ERP(26,:,1,:)),2)-mean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_Waveforms_AllPs -csvwrite('RewP_Waveforms_AllPs.csv',[squeeze(All_ERP(26,:,1,:)),squeeze(All_ERP(26,:,2,:))]'); %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 -csvwrite('RewP_Latency.csv',peak_loc'); %Export data -%% RewP_FFT -csvwrite('RewP_FFT.csv',[(1:1:10)',squeeze(mean(All_FFT(26,:,1,:),4))',squeeze(mean(All_FFT(26,:,2,:),4))',(squeeze(mean(All_FFT(26,:,1,:),4))-squeeze(mean(All_FFT(26,:,2,:),4)))']); %Export data. Conditions: Frequency, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_FFT_AllPs -csvwrite('RewP_FFT_AllPs.csv',[squeeze(All_FFT(26,:,1,:)),squeeze(All_FFT(26,:,2,:))]'); %Export data. Conditions: Loss, Win. -%% RewP_WAV_Stats -delta_extract=zeros(30,600); %Create empty matrix for delta extraction -delta_extract(1:2,:)=(squeeze(mean(All_WAV(26,1:2,151:750,:,:),[4,5])))>5.5; %Determine delta effect via the collapsed localizer method with a power threshhold of 5.5 -theta_extract=zeros(30,600); %Create empty matrix for theta extraction -theta_extract(3:7,:)=(squeeze(mean(All_WAV(26,3:7,151:750,:,:),[4,5])))>5.5; %Determine theta effect via the collapsed localizer method with a power threshhold of 5.5 -both_extract = (delta_extract+theta_extract>0); %Determine all effects via the collapsed localizer method - -WAV_data1 = permute(squeeze(All_WAV(26,:,151:750,1,:)),[3,1,2]); %Extract participants time-frequency condition 1 -WAV_data2 = permute(squeeze(All_WAV(26,:,151:750,2,:)),[3,1,2]); %Extract participants time-frequency condition 2 -WAV_diff = WAV_data1-WAV_data2; %Create difference wave -nb=10; -for participant = 1:nb %Cycle through participants - WAV_Delta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*delta_extract; %Confine data to significant delta activity for difference WAV - WAV_Theta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*theta_extract; %Confine data to significant theta activity for difference WAV - gain_WAV_Delta(participant,:,:) = squeeze(WAV_data1(participant,:,:)).*delta_extract; %Confine data to significant delta activity for condition 1 - loss_WAV_Theta(participant,:,:) = squeeze(WAV_data2(participant,:,:)).*theta_extract; %Confine data to significant theta activity for condition 2 - WAV_Both(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*both_extract; %Determine all frequencies -end - -WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -gain_WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -loss_WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -WAV_Both(WAV_Both==0) = nan; %Remove non-significant data - -WAV_Extract(:,1) = nanmean(WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of difference WAV -WAV_Extract(:,2) = nanmean(WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of difference WAV -WAV_Extract(:,3) = nanmean(loss_WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of condition 2 -WAV_Extract(:,4) = nanmean(gain_WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of condition 1 -csvwrite('RewP_WAV_Stats.csv',WAV_Extract); %Export data -%% RewP_WAV_Freqs -[~,WAV_Extract_freq_time] = max(WAV_diff(:,1:2,:),[],3); %Determine max peak times for delta activity -[~,WAV_Extract_freq_time(:,3:7)] = min(WAV_diff(:,3:7,:),[],3); %Determine min peak times for theta activity -WAV_Time = -200:2:998; %Create a variable of time in milliseconds -csvwrite('RewP_WAV_Freqs.csv',[nanmean(WAV_Both(:,1:7,:),3),WAV_Time(WAV_Extract_freq_time)]); %Export data diff --git a/src/article/peprocess50.m b/src/article/peprocess50.m deleted file mode 100644 index b4295d8..0000000 --- a/src/article/peprocess50.m +++ /dev/null @@ -1,100 +0,0 @@ -%%%% script démarrage pre process copimod2411 -restoredefaultpath; -clear all; -% close all; -clc; %First, clean the environment -addpath(genpath('/home/nforde/Documents/StageEEGpre/src')) - -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/article/MATLAB-EEG-fileIO-master')) -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/article/MATLAB-EEG-icaTools-master')) -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/article/MATLAB-EEG-preProcessing-master')) -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/article/MATLAB-EEG-timeFrequencyAnalysis-master')) - - - -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/eeglab_current/eeglab2021.1')) - - -%% EXTRACTION OF DATA %% -%% Aggregate data across participants -clear all; close all; clc; %First, clean the environment -cd('/home/nforde/Documents/StageEEGpre/data/Stage 2.3. Part 1.1'); %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 - - %FFT frequency resolution is 0.67, so must extract proper frequencies - All_FFT(:,1,:,participant) = mean(EEG.FFT.data(:,1:2,:),2); %Store all the FFT data into a single variable (1 Hz) - All_FFT(:,2,:,participant) = EEG.FFT.data(:,3,:); %Store all the FFT data into a single variable (2 Hz) - All_FFT(:,3,:,participant) = mean(EEG.FFT.data(:,4:5,:),2); %Store all the FFT data into a single variable (3 Hz) - All_FFT(:,4,:,participant) = EEG.FFT.data(:,6,:); %Store all the FFT data into a single variable (4 Hz) - All_FFT(:,5,:,participant) = mean(EEG.FFT.data(:,7:8,:),2); %Store all the FFT data into a single variable (5 Hz) - All_FFT(:,6,:,participant) = EEG.FFT.data(:,9,:); %Store all the FFT data into a single variable (6 Hz) - All_FFT(:,7,:,participant) = mean(EEG.FFT.data(:,10:11,:),2); %Store all the FFT data into a single variable (7 Hz) - All_FFT(:,8,:,participant) = EEG.FFT.data(:,12,:); %Store all the FFT data into a single variable (8 Hz) - All_FFT(:,9,:,participant) = mean(EEG.FFT.data(:,13:14,:),2); %Store all the FFT data into a single variable (9 Hz) - All_FFT(:,10,:,participant) = EEG.FFT.data(:,15,:); %Store all the FFT data into a single variable (10 Hz) - - All_WAV(:,:,:,:,participant) = EEG.WAV.data; %%Store all the WAV data into a single variable -end - -save('All_ERP', 'All_ERP'); %Save ERP Data -save('All_FFT', 'All_FFT'); %Save FFT Data -%Unfortunately, the WAV file is unable to save as it is too large, thus when using the All_WAV variable, this section must be run first -%% Load Data -load('All_ERP.mat'); %Load saved ERP data -load('All_FFT.mat'); %Load saved FFT data -%% RewP_Waveforms -csvwrite('RewP_Waveforms.csv',[(-200:2:998)',mean(squeeze(All_ERP(26,:,1,:)),2),mean(squeeze(All_ERP(26,:,2,:)),2),mean(squeeze(All_ERP(26,:,1,:)),2)-mean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_Waveforms_AllPs -csvwrite('RewP_Waveforms_AllPs.csv',[squeeze(All_ERP(26,:,1,:)),squeeze(All_ERP(26,:,2,:))]'); %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 -csvwrite('RewP_Latency.csv',peak_loc'); %Export data -%% RewP_FFT -csvwrite('RewP_FFT.csv',[(1:1:10)',squeeze(mean(All_FFT(26,:,1,:),4))',squeeze(mean(All_FFT(26,:,2,:),4))',(squeeze(mean(All_FFT(26,:,1,:),4))-squeeze(mean(All_FFT(26,:,2,:),4)))']); %Export data. Conditions: Frequency, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_FFT_AllPs -csvwrite('RewP_FFT_AllPs.csv',[squeeze(All_FFT(26,:,1,:)),squeeze(All_FFT(26,:,2,:))]'); %Export data. Conditions: Loss, Win. -%% RewP_WAV_Stats -delta_extract=zeros(30,600); %Create empty matrix for delta extraction -delta_extract(1:2,:)=(squeeze(mean(All_WAV(26,1:2,151:750,:,:),[4,5])))>5.5; %Determine delta effect via the collapsed localizer method with a power threshhold of 5.5 -theta_extract=zeros(30,600); %Create empty matrix for theta extraction -theta_extract(3:7,:)=(squeeze(mean(All_WAV(26,3:7,151:750,:,:),[4,5])))>5.5; %Determine theta effect via the collapsed localizer method with a power threshhold of 5.5 -both_extract = (delta_extract+theta_extract>0); %Determine all effects via the collapsed localizer method - -WAV_data1 = permute(squeeze(All_WAV(26,:,151:750,1,:)),[3,1,2]); %Extract participants time-frequency condition 1 -WAV_data2 = permute(squeeze(All_WAV(26,:,151:750,2,:)),[3,1,2]); %Extract participants time-frequency condition 2 -WAV_diff = WAV_data1-WAV_data2; %Create difference wave -nb=100; -for participant = 1:nb %Cycle through participants - WAV_Delta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*delta_extract; %Confine data to significant delta activity for difference WAV - WAV_Theta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*theta_extract; %Confine data to significant theta activity for difference WAV - gain_WAV_Delta(participant,:,:) = squeeze(WAV_data1(participant,:,:)).*delta_extract; %Confine data to significant delta activity for condition 1 - loss_WAV_Theta(participant,:,:) = squeeze(WAV_data2(participant,:,:)).*theta_extract; %Confine data to significant theta activity for condition 2 - WAV_Both(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*both_extract; %Determine all frequencies -end - -WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -gain_WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -loss_WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -WAV_Both(WAV_Both==0) = nan; %Remove non-significant data - -WAV_Extract(:,1) = nanmean(WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of difference WAV -WAV_Extract(:,2) = nanmean(WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of difference WAV -WAV_Extract(:,3) = nanmean(loss_WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of condition 2 -WAV_Extract(:,4) = nanmean(gain_WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of condition 1 -csvwrite('RewP_WAV_Stats.csv',WAV_Extract); %Export data -%% RewP_WAV_Freqs -[~,WAV_Extract_freq_time] = max(WAV_diff(:,1:2,:),[],3); %Determine max peak times for delta activity -[~,WAV_Extract_freq_time(:,3:7)] = min(WAV_diff(:,3:7,:),[],3); %Determine min peak times for theta activity -WAV_Time = -200:2:998; %Create a variable of time in milliseconds -csvwrite('RewP_WAV_Freqs.csv',[nanmean(WAV_Both(:,1:7,:),3),WAV_Time(WAV_Extract_freq_time)]); %Export data diff --git a/src/article/sansICA75.m b/src/article/sansICA75.m deleted file mode 100644 index 7a8a077..0000000 --- a/src/article/sansICA75.m +++ /dev/null @@ -1,102 +0,0 @@ -%%%% script démarrage pre process copimod2411 -restoredefaultpath; -clear all; -% close all; -clc; %First, clean the environment -addpath(genpath('/home/StageEEGpre/src')) - -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-fileIO-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-icaTools-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-preProcessing-master')) -addpath(genpath('/home/StageEEGpre/dependencies/article/MATLAB-EEG-timeFrequencyAnalysis-master')) - - -% addpath('/home/StageEEGpre/dependencies/eeglab_current/eeglab2021.1') - -addpath(genpath('/home/dependencies/eeglab_current/eeglab2021.1')) - - - -%% EXTRACTION OF DATA %% -%% Aggregate data across participants -clear all; close all; clc; %First, clean the environment -cd('/home/StageEEGpre/data/Rawdata10'); %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 - - %FFT frequency resolution is 0.67, so must extract proper frequencies - All_FFT(:,1,:,participant) = mean(EEG.FFT.data(:,1:2,:),2); %Store all the FFT data into a single variable (1 Hz) - All_FFT(:,2,:,participant) = EEG.FFT.data(:,3,:); %Store all the FFT data into a single variable (2 Hz) - All_FFT(:,3,:,participant) = mean(EEG.FFT.data(:,4:5,:),2); %Store all the FFT data into a single variable (3 Hz) - All_FFT(:,4,:,participant) = EEG.FFT.data(:,6,:); %Store all the FFT data into a single variable (4 Hz) - All_FFT(:,5,:,participant) = mean(EEG.FFT.data(:,7:8,:),2); %Store all the FFT data into a single variable (5 Hz) - All_FFT(:,6,:,participant) = EEG.FFT.data(:,9,:); %Store all the FFT data into a single variable (6 Hz) - All_FFT(:,7,:,participant) = mean(EEG.FFT.data(:,10:11,:),2); %Store all the FFT data into a single variable (7 Hz) - All_FFT(:,8,:,participant) = EEG.FFT.data(:,12,:); %Store all the FFT data into a single variable (8 Hz) - All_FFT(:,9,:,participant) = mean(EEG.FFT.data(:,13:14,:),2); %Store all the FFT data into a single variable (9 Hz) - All_FFT(:,10,:,participant) = EEG.FFT.data(:,15,:); %Store all the FFT data into a single variable (10 Hz) - - All_WAV(:,:,:,:,participant) = EEG.WAV.data; %%Store all the WAV data into a single variable -end - -save('All_ERP', 'All_ERP'); %Save ERP Data -save('All_FFT', 'All_FFT'); %Save FFT Data -%Unfortunately, the WAV file is unable to save as it is too large, thus when using the All_WAV variable, this section must be run first -%% Load Data -load('All_ERP.mat'); %Load saved ERP data -load('All_FFT.mat'); %Load saved FFT data -%% RewP_Waveforms -csvwrite('RewP_Waveforms.csv',[(-200:2:998)',mean(squeeze(All_ERP(26,:,1,:)),2),mean(squeeze(All_ERP(26,:,2,:)),2),mean(squeeze(All_ERP(26,:,1,:)),2)-mean(squeeze(All_ERP(26,:,2,:)),2)]); %Export data. Conditions: Time, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_Waveforms_AllPs -csvwrite('RewP_Waveforms_AllPs.csv',[squeeze(All_ERP(26,:,1,:)),squeeze(All_ERP(26,:,2,:))]'); %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 -csvwrite('RewP_Latency.csv',peak_loc'); %Export data -%% RewP_FFT -csvwrite('RewP_FFT.csv',[(1:1:10)',squeeze(mean(All_FFT(26,:,1,:),4))',squeeze(mean(All_FFT(26,:,2,:),4))',(squeeze(mean(All_FFT(26,:,1,:),4))-squeeze(mean(All_FFT(26,:,2,:),4)))']); %Export data. Conditions: Frequency, Loss, Win, Difference. Electrode 26 is FCz. -%% RewP_FFT_AllPs -csvwrite('RewP_FFT_AllPs.csv',[squeeze(All_FFT(26,:,1,:)),squeeze(All_FFT(26,:,2,:))]'); %Export data. Conditions: Loss, Win. -%% RewP_WAV_Stats -delta_extract=zeros(30,600); %Create empty matrix for delta extraction -delta_extract(1:2,:)=(squeeze(mean(All_WAV(26,1:2,151:750,:,:),[4,5])))>5.5; %Determine delta effect via the collapsed localizer method with a power threshhold of 5.5 -theta_extract=zeros(30,600); %Create empty matrix for theta extraction -theta_extract(3:7,:)=(squeeze(mean(All_WAV(26,3:7,151:750,:,:),[4,5])))>5.5; %Determine theta effect via the collapsed localizer method with a power threshhold of 5.5 -both_extract = (delta_extract+theta_extract>0); %Determine all effects via the collapsed localizer method - -WAV_data1 = permute(squeeze(All_WAV(26,:,151:750,1,:)),[3,1,2]); %Extract participants time-frequency condition 1 -WAV_data2 = permute(squeeze(All_WAV(26,:,151:750,2,:)),[3,1,2]); %Extract participants time-frequency condition 2 -WAV_diff = WAV_data1-WAV_data2; %Create difference wave -nb=75; -for participant = 1:nb %Cycle through participants - WAV_Delta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*delta_extract; %Confine data to significant delta activity for difference WAV - WAV_Theta(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*theta_extract; %Confine data to significant theta activity for difference WAV - gain_WAV_Delta(participant,:,:) = squeeze(WAV_data1(participant,:,:)).*delta_extract; %Confine data to significant delta activity for condition 1 - loss_WAV_Theta(participant,:,:) = squeeze(WAV_data2(participant,:,:)).*theta_extract; %Confine data to significant theta activity for condition 2 - WAV_Both(participant,:,:) = squeeze(WAV_diff(participant,:,:)).*both_extract; %Determine all frequencies -end - -WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -gain_WAV_Delta(WAV_Delta==0) = nan; %Remove non-significant data -loss_WAV_Theta(WAV_Theta==0) = nan; %Remove non-significant data -WAV_Both(WAV_Both==0) = nan; %Remove non-significant data - -WAV_Extract(:,1) = nanmean(WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of difference WAV -WAV_Extract(:,2) = nanmean(WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of difference WAV -WAV_Extract(:,3) = nanmean(loss_WAV_Theta,[2,3]); %Average data across resulting time and frequencies for theta of condition 2 -WAV_Extract(:,4) = nanmean(gain_WAV_Delta,[2,3]); %Average data across resulting time and frequencies for delta of condition 1 -csvwrite('RewP_WAV_Stats.csv',WAV_Extract); %Export data -%% RewP_WAV_Freqs -[~,WAV_Extract_freq_time] = max(WAV_diff(:,1:2,:),[],3); %Determine max peak times for delta activity -[~,WAV_Extract_freq_time(:,3:7)] = min(WAV_diff(:,3:7,:),[],3); %Determine min peak times for theta activity -WAV_Time = -200:2:998; %Create a variable of time in milliseconds -csvwrite('RewP_WAV_Freqs.csv',[nanmean(WAV_Both(:,1:7,:),3),WAV_Time(WAV_Extract_freq_time)]); %Export data diff --git a/src/automagic/ConvertToAutomagic.m b/src/automagic/ConvertToAutomagic.m deleted file mode 100644 index cf138e5..0000000 --- a/src/automagic/ConvertToAutomagic.m +++ /dev/null @@ -1,11 +0,0 @@ -cd('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set'); %Find and change working folder to raw EEG data - -for participant = 1:500 %Cycle through participants - - disp(['Participant: ', num2str(participant)]) %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - SubjectName = ['Subject' participant_number{2}]; - mkdir(SubjectName); - filename=['set_' participant_number{2} '.set']; - movefile(filename,SubjectName); -end \ No newline at end of file diff --git a/src/automagic/autohbn.m b/src/automagic/autohbn.m deleted file mode 100644 index 2223e4e..0000000 --- a/src/automagic/autohbn.m +++ /dev/null @@ -1,64 +0,0 @@ - -%Parameter=load('/udd/nforde/Nina/StageEEGpre/data/RAWformat_Test4_results/project_state.mat') - - - -addAutomagicPaths(); - -name = 'commandline_project2'; -dataFolder = '/udd/nforde/Nina/StageEEGpre/data/EEG/raw/'; -resultsFolder = '/tmp/data/resultsautohbn/'; -ext = '.mat'; -Params= struct('FilterParams', struct('notch', struct('freq', 50), ... - 'high', struct('freq', 0.5, 'order', []),... % Default order for filtering - 'low', struct('freq', 30, 'order', []),... - 'zapline', struct('freq', 50, 'ncomps', 3)),... - 'CRDParams', struct('ChannelCriterion', 0.85,... - 'LineNoiseCriterion', 4,... - 'BurstCriterion', 5,... - 'WindowCriterion', 0.25, ... - 'Highpass', [0.25 0.75]), ... - 'PrepParams', struct(), ... - 'HighvarParams', struct('sd', 25), ... - 'InterpolationParams', struct('method', 'spherical'), ... - 'RPCAParams', struct('lambda', [], ... % Default lambda by alm_rpca - 'tol', 1e-7, ... - 'maxIter', 1000), ... - 'MARAParams', struct('chanlocMap', containers.Map, ... - 'largeMap', 0, ... - 'high', struct('freq', 1.0, 'order', [])), ... - 'ICLabelParams', struct('brainTher', [], ... - 'muscleTher', [0.8 1], ... - 'eyeTher', [0.8 1], ... - 'heartTher', [0.8 1], ... - 'lineNoiseTher', [0.8 1], ... - 'channelNoiseTher', [0.8 1], ... - 'otherTher', [], ... - 'includeSelected', 0, ... - 'high', struct('freq', 2.0, 'order', [])),... - 'EOGRegressionParams', struct(), ... - 'ChannelReductionParams', struct('tobeExcludedChans', []), ... - 'EEGSystem', struct('name', 'EGI',... - 'sys10_20', 0, ... - 'locFile', '', ... - 'refChan', struct('idx', []), ... - 'fileLocType', '',... - 'eogChans', [1 8 14 17 21 25 32 125 126 127 128], ... - 'powerLineFreq', []), ... - 'Settings', struct('trackAllSteps', 0, ... - 'pathToSteps', '/allSteps.mat',... - 'colormap','Default',... - 'sortChans', 0)... - ); - - -VisualisationParams = struct(); -samplingrate=[]; - -%save(Params, 'config_20211007.mat') - - - -project = Project(name, dataFolder, resultsFolder, ext, Params, VisualisationParams, samplingrate); -project.preprocessAll(); -project.interpolateSelected(); \ No newline at end of file diff --git a/src/automagic/compa_automagic.m b/src/automagic/compa_automagic.m deleted file mode 100644 index ac5b835..0000000 --- a/src/automagic/compa_automagic.m +++ /dev/null @@ -1,62 +0,0 @@ -%compa -%init -% ordre generation fichiers dans allstep: -% 1/ EEGorig -% 2/ EEGPrep -% 3/ EEGFiltered -% 4/ EEG Regressed -% 5/ highvaried -% 6/min varied -% 7/ EEG Final -% -% 8/Regressed - - -%A80 -allgui=load('/tmp/DataResults/dataorigin3_results/A80/allSteps_A80 20170712 0853003.mat'); -allscript= load('/tmp/DataResults/commandline_results/A80/allSteps_A80 20170712 0853003.mat'); -DataField = fieldnames(allgui); - -%guifinal=allgui.(DataField{1}).data; - -for k=1:7 - -%tmp2=allgui.(DataField{k}).data; -allgui2=allgui.(DataField{k}).data; -allscript2=allscript.(DataField{k}).data; - - -diff=allgui2-allscript2; -diff = diff(~isnan(diff)); -mat0=zeros(size(diff)); - -if diff==mat0 - disp 'vrai' - -else disp 'false' - [I_row, I_col] = ind2sub(size(diff),I) -end - -end - -disp 'reduced data now : ' - -reducedgui=load('/tmp/DataResults/dataorigin3_results/A80/reduced2_A80 20170712 0853003.mat'); -reducedscript= load('/tmp/DataResults/commandline_results/A80/reduced2_A80 20170712 0853003.mat'); - -gui=reducedgui.reduced.data; -script=reducedscript.reduced.data; - -diff=gui-script; -diff = diff(~isnan(diff)); -mat0=zeros(size(diff)); - - -if diff==mat0 - disp 'vrai' - -else disp 'false' -end - - - \ No newline at end of file diff --git a/src/automagic/erp2511.m b/src/automagic/erp2511.m deleted file mode 100644 index 1b94806..0000000 --- a/src/automagic/erp2511.m +++ /dev/null @@ -1,47 +0,0 @@ - - -addAutomagicPaths(); - -name = 'script_projecterp'; -dataFolder = '/home/StageEEGpre/data/Raw Data Part 1/'; -resultsFolder = '/home/StageEEGpre/data/Raw Data Part 1_erp2511_results/'; -ext = '.dat'; - -Params= struct('FilterParams', struct('high', struct('freq', 0.1000, 'order', []),... % Default order for filtering - 'low', struct('freq', 35, 'order', [])),... - 'CRDParams', struct([]),... - 'PrepParams', struct('lineFrequencies', 50), ... - 'HighvarParams', struct('sd', 25,... - 'cutoff', 100, ... - 'rejRatio', 0.5000), ... - 'MinvarParams', struct('sd', 1),... - 'InterpolationParams', struct('method', 'spherical'), ... - 'ICLabelParams', struct([]),... - 'EEGSystem', struct('name', 'EGI',... - 'sys10_20', 0, ... - 'locFile', '', ... - 'refChan', struct('idx', []), ... - 'fileLocType', '',... - 'eogChans', [], ... - 'tobeExcludedChans', []), ... - 'Settings', struct('trackAllSteps', 1, ... - 'pathToSteps', '/allSteps.mat',... - 'colormap','Default',... - 'sortChans', 0)... - ); - - -sRate= []; -VisualisationParams = struct(); - - - - - -% Template for all the parameters - - - -project = Project(name, dataFolder, resultsFolder, ext, Params, VisualisationParams,sRate); -project.preprocessAll(); -%project.interpolateSelected(); \ No newline at end of file diff --git a/src/automagic/janvier/pbiclabel.m b/src/automagic/janvier/pbiclabel.m deleted file mode 100644 index 021404e..0000000 --- a/src/automagic/janvier/pbiclabel.m +++ /dev/null @@ -1,54 +0,0 @@ -%Parameter=load('/udd/nforde/Nina/StageEEGpre/data/RAWformat_Test4_results/project_state.mat') -restoredefaultpath; -clear all; -clc; %First, clean the environment - -addpath( '/home/nforde/Documents/StageEEGpre/dependencies/automagic-master'); - - -addAutomagicPaths(); - -name = 'testsansICla'; -dataFolder = '/home/nforde/Documents/StageEEGpre/data/formatBIDS'; -resultsFolder = '/home/nforde/Documents/StageEEGpre/data/formatBIDS_testsansICla_results'; -ext = '.dat'; - -% Params= struct('FilterParams', struct('high', struct('freq', 0.1000, 'order', []),... % Default order for filtering -% 'low', struct('freq', 35, 'order', [])),... -% 'CRDParams', struct([]),... -% 'PrepParams', struct('lineFrequencies', 50), ... -% 'HighvarParams', struct('sd', 25,... -% 'cutoff', 100, ... -% 'rejRatio', 0.5000), ... -% 'MinvarParams', struct('sd', 1),... -% 'InterpolationParams', struct('method', 'spherical'), ... -% 'ICLabelParams', struct([]),... -% 'EEGSystem', struct('name', 'EGI',... -% 'sys10_20', 0, ... -% 'locFile', '', ... -% 'refChan', struct('idx', []), ... -% 'fileLocType', '',... -% 'eogChans', [], ... -% 'tobeExcludedChans', []), ... -% 'Settings', struct('trackAllSteps', 1, ... -% 'pathToSteps', '/allSteps.mat',... -% 'colormap','Default',... -% 'sortChans', 0)... -% ); - -Params=load('/home/nforde/Documents/StageEEGpre/data/formatBIDS_testsansICla_results/project_state.mat'); -Params=Params.self.params; -sRate= []; -VisualisationParams = struct(); - - - - - -% Template for all the parameters - - - -project = Project(name, dataFolder, resultsFolder, ext, Params, VisualisationParams,sRate); -project.preprocessAll(); -%project.interpolateSelected(); \ No newline at end of file diff --git a/src/automagic/janvier/script25janv.m b/src/automagic/janvier/script25janv.m deleted file mode 100644 index 1f0ca52..0000000 --- a/src/automagic/janvier/script25janv.m +++ /dev/null @@ -1,54 +0,0 @@ -%Parameter=load('/udd/nforde/Nina/StageEEGpre/data/RAWformat_Test4_results/project_state.mat') -restoredefaultpath; -clear all; -clc; %First, clean the environment - -addpath( '/home/nforde/Documents/StageEEGpre/dependencies/automagic-master'); - - -addAutomagicPaths(); - -name = '24janv2'; -dataFolder = '/home/nforde/Documents/StageEEGpre/data/formatBIDS'; -resultsFolder = '/home/nforde/Documents/StageEEGpre/data/formatBIDS_24janv_results2'; -ext = '.dat'; - -% Params= struct('FilterParams', struct('high', struct('freq', 0.1000, 'order', []),... % Default order for filtering -% 'low', struct('freq', 35, 'order', [])),... -% 'CRDParams', struct([]),... -% 'PrepParams', struct('lineFrequencies', 50), ... -% 'HighvarParams', struct('sd', 25,... -% 'cutoff', 100, ... -% 'rejRatio', 0.5000), ... -% 'MinvarParams', struct('sd', 1),... -% 'InterpolationParams', struct('method', 'spherical'), ... -% 'ICLabelParams', struct([]),... -% 'EEGSystem', struct('name', 'EGI',... -% 'sys10_20', 0, ... -% 'locFile', '', ... -% 'refChan', struct('idx', []), ... -% 'fileLocType', '',... -% 'eogChans', [], ... -% 'tobeExcludedChans', []), ... -% 'Settings', struct('trackAllSteps', 1, ... -% 'pathToSteps', '/allSteps.mat',... -% 'colormap','Default',... -% 'sortChans', 0)... -% ); - -Params=load('/home/nforde/Documents/StageEEGpre/data/formatBIDS_24janv_results2/project_state.mat'); -Params=Params.self.params; -sRate= []; -VisualisationParams = struct(); - - - - - -% Template for all the parameters - - - -project = Project(name, dataFolder, resultsFolder, ext, Params, VisualisationParams,sRate); -project.preprocessAll(); -%project.interpolateSelected(); \ No newline at end of file diff --git a/src/automagic/janvier/segmentation.m b/src/automagic/janvier/segmentation.m deleted file mode 100644 index a660606..0000000 --- a/src/automagic/janvier/segmentation.m +++ /dev/null @@ -1,37 +0,0 @@ -%dd -restoredefaultpath; -clear all; -clc; - -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/article/MATLAB-EEG-preProcessing-master'))% pour doSegment -addpath(genpath('/home/nforde/Documents/StageEEGpre/dependencies/eeglab_current/eeglab2021.1'))% pour pop rmbase - - -dataFolder = '/home/nforde/Documents/StageEEGpre/data/formatBIDS_24janv_results2/sujets2'; -cd(dataFolder) -filenames = dir('allSteps_RewardProcessing*') %Compile list of all data - - -for participant = 1:2%Cycle through participants - - %Get participant name information - disp(['Participant: ', num2str(participant)]) %Display current participant being processed - participant_number = strsplit(filenames(participant).name(1:end-4),'_') %Split filename into components - participant_varname = ['segmented',participant_number{3}] %Create new file name - EEG = []; %Clear past data - file=load(filenames(participant).name); %Load participant output - EEG=file.EEGFinal - - - - - %Determine markers of interest for WAV later - markers = {'S110','S111'}; %Loss, win - - - [EEG] = doSegmentData(EEG,markers,[-500 1498]); %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 correction - EEG = pop_rmbase( EEG, [-200/1000,0]); - save(participant_varname,'EEG') %Save the current output - -end diff --git a/src/automagic/testauto.m b/src/automagic/testauto.m deleted file mode 100644 index 8d80273..0000000 --- a/src/automagic/testauto.m +++ /dev/null @@ -1,44 +0,0 @@ -addAutomagicPaths(); - -name = 'commandline_project'; -dataFolder = '/tmp/RAWformat/'; -resultsFolder = '/tmp/DataResults/commandline_results/'; -ext = '.RAW'; -Params= struct('FilterParams', struct('high', struct('freq', 0.1000, 'order', []),... % Default order for filtering - 'low', struct('freq', 35, 'order', [])),... - 'CRDParams', struct([]),... - 'PrepParams', struct('lineFrequencies', 50), ... - 'HighvarParams', struct('sd', 25,... - 'cutoff', 100, ... - 'rejRatio', 0.5000), ... - 'MinvarParams', struct('sd', 1),... - 'InterpolationParams', struct('method', 'spherical'), ... - 'ICLabelParams', struct([]),... - 'EEGSystem', struct('name', 'EGI',... - 'sys10_20', 0, ... - 'locFile', '', ... - 'refChan', struct('idx', []), ... - 'fileLocType', '',... - 'eogChans', [], ... - 'tobeExcludedChans', []), ... - 'Settings', struct('trackAllSteps', 1, ... - 'pathToSteps', '/allSteps.mat',... - 'colormap','Default',... - 'sortChans', 0)... - ); - - -sRate= []; -VisualisationParams = struct(); - - - - - -% Template for all the parameters - - - -project = Project(name, dataFolder, resultsFolder, ext, Params, VisualisationParams,sRate); -project.preprocessAll(); -%project.interpolateSelected(); \ No newline at end of file diff --git a/src/fieldtrip/ftread.m b/src/fieldtrip/ftpreprocessing.m similarity index 72% rename from src/fieldtrip/ftread.m rename to src/fieldtrip/ftpreprocessing.m index 5c96488..7ece5b7 100644 --- a/src/fieldtrip/ftread.m +++ b/src/fieldtrip/ftpreprocessing.m @@ -1,19 +1,23 @@ -% % load channel list -load('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/Subject001/channelsTokeep.mat'); +% % --- This is the preprocessing workflow reproduced by FieldTrip functions--- +% % For contact: aya.kabbara7@gmail.com + +% % Regulate some parameters that will be used in interpolation +load('src/channelsTokeep.mat'); + % % prepare neighborhood template for bad channels interpolation reduced_subjects=ss; cfg = []; -cfg.dataset = '/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/Subject001/set_001.set'; +cfg.dataset = '/Raw Data Part 1/set/Subject001/set_001.set'; cfg.continuous = 'yes'; % force it to be continuous -cfg.channel = ChannelsTokeep; +cfg.channel = ChannelsTokeep; data_prepare = ft_preprocessing(cfg); cfg=[]; cfg.method ='distance'; cfg.neighbourdist=0.7; [neighbours, cfg] = ft_prepare_neighbours(cfg, data_prepare) -% % go to folder -cd('/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1'); %Find and change working folder to raw EEG data +% % go to data folder +cd('/Raw Data Part 1'); %Find and change working folder to raw EEG data filenames = dir('*.vhdr') nb=500; trials_loss=[]; @@ -25,7 +29,7 @@ disp(['Participant: ', num2str(participant)]) %Display current participant being processed participant_number = strsplit(filenames(participant).name(1:end-5),'_'); %Split filename into components - RawFile = ['/Users/ayakabbara/Desktop/projects/EEG_PreProcessing/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; + RawFile = ['/Raw Data Part 1/set/Subject' participant_number{2} '/set_' participant_number{2} '.set']; SubjectName = ['participant_' participant_number{2}]; % % Check if the folder contains the required files @@ -35,15 +39,15 @@ % % read data cfg = []; cfg.dataset = RawFile; - cfg.continuous = 'yes'; + cfg.continuous = 'yes'; % force it to be continuous if(length(find(reduced_subjects==participant))==0) - cfg.channel = ChannelsTokeep; + cfg.channel = ChannelsTokeep; end % % re-referencing cfg.reref = 'yes'; - cfg.refchannel = {'TP9', 'TP10'}; + cfg.refchannel = {'TP9', 'TP10'}; data_eeg = ft_preprocessing(cfg); @@ -91,7 +95,7 @@ try [data_interp] = ft_channelrepair(cfg, data_eeg_filtered) - % % define trials, segment + % % define trials, segment cfg = []; cfg.dataset= RawFile; cfg.trialdef.eventtype = 'Stimulus'; @@ -151,18 +155,34 @@ try [timelock] = ft_timelockanalysis(cfg, data_loss_final); trials_loss(participant)=length(data_loss_final.trial); - All_ERP_ft(1,:,:,participant) = timelock.avg; + All_ERP_ft(1,:,:,participant) = timelock.avg; catch end try [timelock] = ft_timelockanalysis(cfg, data_win_final); trials_win(participant)=length(data_win_final.trial); - All_ERP_ft(2,:,:,participant) = timelock.avg; + All_ERP_ft(2,:,:,participant) = timelock.avg; catch end catch end end - save('All_ERP_ft_final.mat','All_ERP_ft'); \ No newline at end of file + 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,:,:)); + +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('ft_RewP_Waveforms_AllPs_final22.csv',[tt1,tt2]'); %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_final22.csv',peak_loc'); %Export data +