From 1ab0863a84b433e9aa9ac3d83f1530ca1475e22e Mon Sep 17 00:00:00 2001 From: Felix Zakirov <fzaki001@fiu.edu> Date: Sun, 8 Dec 2024 23:39:45 -0500 Subject: [PATCH] Updated scripts --- code/preprocessing-eeg/create_mat.m | 63 ++++++------- code/preprocessing-eeg/create_mat_csd.m | 91 +++++++------------ .../tf/MainScript_Calculate_TF_ITPS_ICPS.m | 38 +++++--- .../MainScript_Calculate_TF_ITPS_ICPS_para.m | 35 ++++--- code/preprocessing-eeg/tf/timefreq_script.m | 5 +- code/preprocessing-eeg/transform_csd.m | 70 ++++++++------ 6 files changed, 152 insertions(+), 150 deletions(-) diff --git a/code/preprocessing-eeg/create_mat.m b/code/preprocessing-eeg/create_mat.m index 660f442..42766d1 100644 --- a/code/preprocessing-eeg/create_mat.m +++ b/code/preprocessing-eeg/create_mat.m @@ -14,6 +14,8 @@ % addpath(genpath([main_dir filesep 'MADE-EEG-preprocessing-pipeline']));% enter the path of the EEGLAB folder in this line addpath(genpath('/home/data/NDClab/tools/lab-devOps/scripts/MADE_pipeline_standard/eeg_preprocessing'));% enter the path of the folder in this line +addpath(genpath('/home/data/NDClab/analyses/thrive-theta-ddm/code/matlab/')) + %Location of "EEG % addpath(genpath([main_dir filesep 'eeglab13_4_4b']));% enter the path of the EEGLAB folder in this line addpath(genpath('/home/data/NDClab/tools/lab-devOps/scripts/MADE_pipeline_standard/eeglab13_4_4b'));% enter the path of the EEGLAB folder in this line @@ -30,7 +32,7 @@ %location of dataset folder % dataset_dir = '/Users/fzaki001/thrive-theta-ddm'; -dataset_dir = '/home/data/NDClab/datasets/thrive-dataset/derivatives/preprocessed/'; +dataset_dir = '/home/data/NDClab/datasets/thrive-dataset'; % summary_csv_path = '/Users/fzaki001/thrive-theta-ddm/derivatives/behavior/summary.csv'; summary_csv_path = '/home/data/NDClab/analyses/thrive-theta-ddm/derivatives/behavior/summary.csv'; @@ -81,7 +83,10 @@ % outputHeader = {'id, s_resp_incon_error, s_resp_incon_corr, ns_resp_incon_error, ns_resp_incon_corr'}; % dlmwrite(strcat('thrive_trialCounts_respOnly', date, '.csv'), outputHeader, 'delimiter', '', '-append'); +diary(sprintf('erp_log_%s.log', datestr(now, 'mm_dd_yyyy_HH_MM_SS'))) + for subject=1:length(datafile_names) + EEG=[]; fprintf('\n\n\n*** Processing subject %d (%s) ***\n\n\n', subject, datafile_names{subject}); @@ -138,14 +143,14 @@ %count how many of each event type (combination of event types) of %interest are present - s_resp_incon_error = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - s_resp_incon_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - ns_resp_incon_error = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - ns_resp_incon_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - s_stim_incon_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - s_stim_con_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "c")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - ns_stim_incon_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - ns_stim_con_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "c")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); + s_resp_incon_error = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + s_resp_incon_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + ns_resp_incon_error = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + ns_resp_incon_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + s_stim_incon_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + s_stim_con_corr = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "c")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + ns_stim_incon_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + ns_stim_con_corr = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "stim")) & (strcmp({EEG.event.congruency}, "c")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); %Create the trial counts table for trial counts counts_table=table({datafile_names{subject}}, {s_resp_incon_error}, {s_resp_incon_corr}, {ns_resp_incon_error}, {ns_resp_incon_corr}, {s_stim_incon_corr}, {s_stim_con_corr}, {ns_stim_incon_corr}, {ns_stim_con_corr}); @@ -165,7 +170,7 @@ %specify min number of trials per condition (if file contains less than %this number for ANY condition, then they will be skipped for ALL conditions -minTrials = 8; +minTrials = 6; %specify min accuracy per condition (if file contains less than %this number for ANY condition, then they will be skipped for ALL conditions @@ -208,32 +213,18 @@ EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on'); EEG = eeg_checkset( EEG ); - % NOTE % - %the logic of checking conditions and then looping over conditions - %below is fairly hard-coded and could be be substantially improved to - %allow for easier reuse when number of conditions or number of - %variables per condition changes. - % - %before pulling trials of interest, for any conditions, check to make - %sure this file/participant has more than minTrials for EACH condition. - %If the file/participant is below minTrials for even one of the - %conditions that will be pulled, then the file/participant is skipped - %entirely and no condition data at all will be pulled for this specific - %file (but the participant can still have data pulled for another one - %of their files from another visit). - % %count trials for each condition of interest and store in numTrials vector - numTrials(1) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(2) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(3) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(4) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); + numTrials(1) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(2) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(3) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(4) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); %logical test if the number of trials for each condition (numTrials vector) %are NOTE all >= minTrials. If statement is true, then participant/file %is skipped and for loop over files continues to next file - % if ~(sum(numTrials >= minTrials) == length(numTrials)) - % continue - % end + if ~(sum(numTrials >= minTrials) == length(numTrials)) + continue + end % loop through conditions of interest for this file (combo of event types) % @@ -251,6 +242,7 @@ accuracy = 0; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==2) % social correct observation = 's'; eventType = 'resp'; @@ -258,6 +250,7 @@ accuracy = 1; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==3) % nonsocial error observation = 'ns'; eventType = 'resp'; @@ -265,6 +258,7 @@ accuracy = 0; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==4) % nonsocial correct observation = 'ns'; eventType = 'resp'; @@ -272,11 +266,12 @@ accuracy = 1; responded = 1; validRt = 1; + extraResponse = 0; end %select combintion of event types of interest based on vars above - EEG1 = pop_selectevent( EEG, 'latency','-1<=1','observation',observation,'eventType',eventType,'congruency',congruency,'accuracy',accuracy,'responded',responded,'validRt',validRt,'deleteevents','on','deleteepochs','on','invertepochs','off'); - EEG1 = eeg_checkset( EEG1 ); + EEG1 = pop_selectevent( EEG, 'latency', '-1<=1', 'observation', observation, 'eventType', eventType, 'congruency', congruency, 'accuracy', accuracy, 'responded', responded, 'validRt', validRt, 'extraResponse', extraResponse, 'deleteevents', 'on', 'deleteepochs', 'on', 'invertepochs', 'off'); + EEG1 = eeg_checkset(EEG1); % Average across epoch dimension % this all Channel ERP only needs to be computed once @@ -300,4 +295,4 @@ end %save the erps and subject list -save('thrive_Resp_erps_min_8t_60acc.mat','erpDat_data', 'erpDat_subIds') \ No newline at end of file +save(sprintf('thrive_Resp_erps_min_6t_60acc_%s.mat', datestr(now, 'mm_dd_yyyy_HH_MM_SS')), 'erpDat_data', 'erpDat_subIds') diff --git a/code/preprocessing-eeg/create_mat_csd.m b/code/preprocessing-eeg/create_mat_csd.m index c82eab6..b4f9c88 100644 --- a/code/preprocessing-eeg/create_mat_csd.m +++ b/code/preprocessing-eeg/create_mat_csd.m @@ -31,16 +31,16 @@ analysis_dir = '/home/data/NDClab/analyses/thrive-theta-ddm'; %location of dataset folder -% dataset_dir = '/Users/fzaki001/thrive-theta-ddm'; -dataset_dir = '/home/data/NDClab/datasets/thrive-dataset'; +dataset_dir = '/home/data/NDClab/analyses/thrive-theta-ddm'; +%dataset_dir = '/home/data/NDClab/datasets/thrive-dataset'; % summary_csv_path = '/Users/fzaki001/thrive-theta-ddm/derivatives/behavior/summary.csv'; summary_csv_path = '/home/data/NDClab/analyses/thrive-theta-ddm/derivatives/behavior/summary.csv'; % Setting up other things % 1. Enter the path of the folder that has the data to be analyzed -data_location = [dataset_dir filesep 'derivatives' filesep 'preprocessed']; - +%data_location = [dataset_dir filesep 'derivatives' filesep 'preprocessed']; +data_location = [dataset_dir filesep 'derivatives' filesep 'preprocessed' filesep 'csd_data']; % 2. Enter the path of the folder where you want to save the postprocessing outputs output_location = [analysis_dir filesep 'derivatives' filesep 'preprocessed/erp_check']; @@ -58,7 +58,8 @@ visitFileName = 's1_r1_e1'; %file names include "e1" designation % Read files to analyses -datafile_info=dir([data_location filesep 'sub-*' filesep visitDirName filesep 'eeg' filesep 'sub-*_' task '_eeg_*' procStage '_' visitFileName '.set']); +%datafile_info=dir([data_location filesep 'sub-*' filesep visitDirName filesep 'eeg' filesep 'sub-*_' task '_eeg_*' procStage '_' visitFileName '.set']); +datafile_info=dir([data_location filesep 'sub-*_' task '_eeg_*' procStage '_' visitFileName '.set']); datafile_info=datafile_info(~ismember({datafile_info.name},{'.', '..', '.DS_Store'})); datafile_names={datafile_info.name}; datafile_paths={datafile_info.folder}; @@ -72,15 +73,15 @@ % Create output folders to save data if exist(output_location, 'dir') == 0 - mkdir(output_location) + mkdir(output_location); end -for site = 1:64 - trodes{site} = num2str(site) -end -Montage_64=ExtractMontage('/home/data/NDClab/analyses/thrive-theta-ddm/code/preprocessing-eeg/64ch_bv_montage_csd.csd', trodes'); +%for site = 1:64 +% trodes{site} = num2str(site); +%end +%Montage_64=ExtractMontage('/home/data/NDClab/analyses/thrive-theta-ddm/code/preprocessing-eeg/64ch_bv_montage_csd.csd', trodes'); % MapMontage(Montage_64); -[G, H] = GetGH(Montage_64); +%[G, H] = GetGH(Montage_64); %% Count trials % switch to output directory @@ -90,7 +91,7 @@ % outputHeader = {'id, s_resp_incon_error, s_resp_incon_corr, ns_resp_incon_error, ns_resp_incon_corr'}; % dlmwrite(strcat('thrive_trialCounts_respOnly', date, '.csv'), outputHeader, 'delimiter', '', '-append'); -diary(sprintf('erp_log_%s.log', datestr(now, 'mm-dd-yyyy_HH_MM_SS'))) +diary(sprintf('erp_log_%s.log', datestr(now, 'mm_dd_yyyy_HH_MM_SS'))) %% pull resp-locked erp mat file @@ -125,16 +126,6 @@ %find row in behavior file corresponding to this participant behavior_id_match_idxs = find(behavior_info{:,'sub'} == str2num(subNumText)); - %if participant has low accuracy in either condition, skip that - %participant for ALL conditions - % if (behavior_info{behavior_id_match_idxs,'acc_nonsoc'} < acc_cutoff || behavior_info{behavior_id_match_idxs,'acc_soc'} < acc_cutoff) - % continue - % end - % - % if (behavior_info{behavior_id_match_idxs,'x6_or_more_err_nonsoc'} < 6 || behavior_info{behavior_id_match_idxs,'x6_or_more_err_soc'} < 6) - % continue - % end - %load the original data set EEG = pop_loadset( 'filename', datafile_names{subject}, 'filepath', datafile_paths{subject}); EEG = eeg_checkset( EEG ); @@ -142,34 +133,20 @@ EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on'); EEG = eeg_checkset( EEG ); - for ne = 1:length(EEG.epoch) - myEEG = single(EEG.data(:, :, ne)); - MyResults = CSD(myEEG, G, H); % compute CSD for <channels-by-samples> 2-D epoch - data(:, :, ne) = MyResults; - end - EEG.data = data - - data(:,:,:) = NaN; - - % NOTE % - %the logic of checking conditions and then looping over conditions - %below is fairly hard-coded and could be be substantially improved to - %allow for easier reuse when number of conditions or number of - %variables per condition changes. - % - %before pulling trials of interest, for any conditions, check to make - %sure this file/participant has more than minTrials for EACH condition. - %If the file/participant is below minTrials for even one of the - %conditions that will be pulled, then the file/participant is skipped - %entirely and no condition data at all will be pulled for this specific - %file (but the participant can still have data pulled for another one - %of their files from another visit). - % + % for ne = 1:length(EEG.epoch) + % myEEG = single(EEG.data(:, :, ne)); + % MyResults = CSD(myEEG, G, H); % compute CSD for <channels-by-samples> 2-D epoch + % data(:, :, ne) = MyResults; + % end + % EEG.data = data + + % data(:,:,:) = NaN; + %count trials for each condition of interest and store in numTrials vector - numTrials(1) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(2) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(3) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); - numTrials(4) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) )); + numTrials(1) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(2) = length(find( (strcmp({EEG.event.observation}, "s")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(3) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 0) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); + numTrials(4) = length(find( (strcmp({EEG.event.observation}, "ns")) & (strcmp({EEG.event.eventType}, "resp")) & (strcmp({EEG.event.congruency}, "i")) & ([EEG.event.accuracy] == 1) & ([EEG.event.responded] == 1) & ([EEG.event.validRt] == 1) & ([EEG.event.extraResponse] == 0) )); %logical test if the number of trials for each condition (numTrials vector) %are NOTE all >= minTrials. If statement is true, then participant/file @@ -178,13 +155,7 @@ continue end - % loop through conditions of interest for this file (combo of event types) - % - % specify number of conditions using a seperate conditionNums var, so - % that it can be referenced below when iterating idx counters (to only - %iterate when c == length(conditionNums); conditionNums = 1:4; - % for c = conditionNums if (c==1) % social error @@ -194,6 +165,7 @@ accuracy = 0; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==2) % social correct observation = 's'; eventType = 'resp'; @@ -201,6 +173,7 @@ accuracy = 1; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==3) % nonsocial error observation = 'ns'; eventType = 'resp'; @@ -208,6 +181,7 @@ accuracy = 0; responded = 1; validRt = 1; + extraResponse = 0; elseif (c==4) % nonsocial correct observation = 'ns'; eventType = 'resp'; @@ -215,10 +189,11 @@ accuracy = 1; responded = 1; validRt = 1; + extraResponse = 0; end - %select combintion of event types of interest based on vars above - EEG1 = pop_selectevent( EEG, 'latency','-1<=1','observation',observation,'eventType',eventType,'congruency',congruency,'accuracy',accuracy,'responded',responded,'validRt',validRt,'deleteevents','on','deleteepochs','on','invertepochs','off'); + %select combination of event types of interest based on vars above + EEG1 = pop_selectevent( EEG, 'latency','-1<=1','observation', observation, 'eventType', eventType, 'congruency', congruency, 'accuracy', accuracy, 'responded', responded,'validRt', validRt, 'extraResponse', extraResponse, 'deleteevents','on','deleteepochs','on','invertepochs','off'); EEG1 = eeg_checkset( EEG1 ); % Average across epoch dimension @@ -243,4 +218,4 @@ end %save the erps and subject list -save('thrive_Resp_erps_csd_min_6t_60acc.mat','erpDat_data', 'erpDat_subIds') +save(sprintf('thrive_Resp_erps_csd_min_6t_60acc_%s.mat', datestr(now, 'mm_dd_yyyy_HH_MM_SS')), 'erpDat_data', 'erpDat_subIds') diff --git a/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS.m b/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS.m index 5dffa4b..ffb6413 100644 --- a/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS.m +++ b/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS.m @@ -1,25 +1,30 @@ %%to Run on FIU HPC% % create a local cluster object -cluster = parcluster('local'); +%cluster = parcluster('local'); % start matlabpool with max workers set in the slurm file -parpool(cluster, str2num(getenv('SLURM_CPUS_ON_NODE'))) +%parpool(cluster, str2num(getenv('SLURM_CPUS_ON_NODE'))) %Compute TF, ITPS, ICPS, and wPLI measures for EEG data %Maureen Bowers 6/29/2021 based on scripts by Ranjan Debnath -% Kianoosh Hosseini has edited this script for mfe_b study on 07/18/2023! + clear all; clc; %% %%%%% Setting paths %%%%% -main_dir = '/home/data/NDClab/datasets/thrive-theta-ddm'; %directory on the HPC +main_dir = '/home/data/NDClab/analyses/thrive-theta-ddm'; %directory on the HPC %1. Data Location data_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'csd_data']; %2. Save Data Location -save_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'eeg' filesep 'TF_outputs' filesep 'main' filesep ]; +save_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'TF_outputs' filesep 'test_fix_seed']; +%disp(save_location) +% Create output folders to save data +if exist(save_location, 'dir') == 0 + mkdir(save_location); +end %3. Scripts Location scripts_location = [main_dir filesep 'code' filesep 'preprocessing-eeg' filesep 'tf']; @@ -122,10 +127,15 @@ eeglab % Loading EEGLAB +%rng('default'); % Reset random number generator +%s = RandStream('mlfg6331_64'); % Create a random number stream + +rng(2,"twister") + %%%%%%%%%%%%%%%%%%%% COMPUTATIONS BEGIN BELOW HERE %%%%%%%%%%%%%%%% %% loop through all subject -for sub=1:5 - +for sub=1:2 + % Use the same stream for all workers % Initialize objects for this participant: timefreqs_data = []; phase_data=[]; @@ -141,8 +151,8 @@ EEG=pop_loadset('filename', [subject], 'filepath', data_location); EEG = pop_selectevent( EEG, 'latency','-.1 <= .1','deleteevents','on'); - EEG = pop_editeventfield( EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Condition','NaN'); - EEG = eeg_checkset( EEG ); + EEG = pop_editeventfield(EEG, 'indices', strcat('1:', int2str(length(EEG.event))), 'Condition','NaN'); + EEG = eeg_checkset(EEG); for t=1:length(EEG.event) @@ -181,8 +191,8 @@ % Trials that are labeled as responded and validRt (trials that have larger than 150 ms RT) % will only be included. try - EEG = pop_selectevent( EEG, 'validRt', 1, 'extraResponse', 0, deleteevents','on','deleteepochs','on','invertepochs','off'); - EEG = eeg_checkset( EEG ); + EEG = pop_selectevent(EEG, 'validRt', 1, 'extraResponse', 0, 'deleteevents', 'on', 'deleteepochs', 'on', 'invertepochs', 'off'); + EEG = eeg_checkset(EEG); catch ME warning('Error occurred when selecting only valid trials for subject %d: %s', subject, ME.message); @@ -190,8 +200,8 @@ end try - EEG = pop_selectevent( EEG, 'responded',1, 'deleteevents','on','deleteepochs','on','invertepochs','off'); - EEG = eeg_checkset( EEG ); + EEG = pop_selectevent(EEG, 'responded', 1, 'deleteevents', 'on', 'deleteepochs', 'on', 'invertepochs', 'off'); + EEG = eeg_checkset(EEG); catch ME warning('Error occurred when selecting only valid trials (responded) for subject %d: %s', subject, ME.message); @@ -214,4 +224,4 @@ %Save out table with trial numbers TrialNums_table = struct2table(TrialNums); -writetable(TrialNums_table,[save_location 'TrialNums.csv']); \ No newline at end of file +writetable(TrialNums_table,[save_location 'TrialNums.csv']); diff --git a/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS_para.m b/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS_para.m index 85dd898..e402950 100644 --- a/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS_para.m +++ b/code/preprocessing-eeg/tf/MainScript_Calculate_TF_ITPS_ICPS_para.m @@ -19,7 +19,7 @@ data_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'csd_data']; %2. Save Data Location -save_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'TF_outputs' filesep 'main' filesep 'resp' filesep]; +save_location = [main_dir filesep 'derivatives' filesep 'preprocessed' filesep 'TF_outputs' filesep 'main' filesep 'stim' filesep]; %disp(save_location) % Create output folders to save data if exist(save_location, 'dir') == 0 @@ -51,20 +51,33 @@ %7. What are your conditions of interest if using Event-Related Data? This %naming convention should come from the Edit_events.m script provided. %Note: not needed for resting state data. - +%Super_Conds = { +%{'resp_s_i_0'},... +%{'resp_s_i_1'},... +%{'resp_s_c_1'},... +%{'resp_ns_i_0'},... +%{'resp_ns_i_1'},... +%{'resp_ns_c_1'},... +%} % RESPONSE CONDITIONS -Conds = { +%Conds = { %'resp_s_i_0',... %'resp_s_i_1',... %'resp_s_c_1',... %'resp_ns_i_0',... %'resp_ns_i_1',... -'resp_ns_c_1',... -}; % resp_i_0 = incong error response; resp_i_1 = incong correct response +%'resp_ns_c_1',... +%}; % resp_i_0 = incong error response; resp_i_1 = incong correct response % STIMULUS CONDITIONS -%Conds = {'stim_s_i_0','stim_s_i_1', 'stim_s_c_1', ... -% 'stim_ns_i_0','stim_ns_i_1', 'stim_ns_c_1',}; % resp_i_0 = incong error stim; resp_i_1 = incong correct response +Conds = { +%'stim_s_i_0',... +%'stim_s_i_1',... +%'stim_s_c_1',... +%'stim_ns_i_0',... +%'stim_ns_i_1',... +'stim_ns_c_1',... +}; % resp_i_0 = incong error stim; resp_i_1 = incong correct response %8. Minimum number of trials to analyze mintrialnum = 6; %If the participant does not have enough trials in a condition based on this cutoff, a "notenoughdata.mat" file will be saved into save_location. @@ -97,7 +110,7 @@ %%%%% Questions about phase-based measures %%%%% %14. Would you like to calculate inter-trial phase synchrony (ITPS) in addition to TF? -ITPS_calc = 0; %(1=yes,0=no) +ITPS_calc = 1; %(1=yes,0=no) %15. Would you like to subsample trials? This is recommended for event-related paradigms, % especially when there are uneven numbers of trials in conditions. @@ -114,7 +127,7 @@ ICPS_calc = 1; %(1=yes,0=no) %17. Would you like to calculate coherence or weighted phaselagidx? -ICPS_or_wPLI = 0; %(1=coherence, 0=wPLI) +ICPS_or_wPLI = 1; %(1=coherence, 0=wPLI) %18. Inter-channel phase synchrony over trials or connectivity over time? % Kia: over trials results in a frequency by time matrix and over time @@ -146,12 +159,8 @@ TrialNums = zeros(length(subject_list), length(Conds)); %rng(2, 'twister'); % to fix seed -%D = parallel.pool.DataQueue; -%waitbar = waitbar(0, 'Processing subjects...'); -%afterEach(D, @(x) waitbar(x/length(subject_list), waitbar)); %%%%%%%%%%%%%%%%%%%% COMPUTATIONS BEGIN BELOW HERE %%%%%%%%%%%%%%%% %% loop through all subject - parfor sub=1:length(subject_list) %parfor sub=1:2 % rng(2, 'twister'); % if you want to fix the seed put it inside the loop as well diff --git a/code/preprocessing-eeg/tf/timefreq_script.m b/code/preprocessing-eeg/tf/timefreq_script.m index 249d36a..601f1b0 100644 --- a/code/preprocessing-eeg/tf/timefreq_script.m +++ b/code/preprocessing-eeg/tf/timefreq_script.m @@ -190,8 +190,9 @@ %pull out condition data cond_data = []; cond_EEG=[]; - cond_EEG = pop_selectevent( EEG, 'Condition',(Conds{cond}), 'deleteevents','on','deleteepochs','on','invertepochs','off'); - cond_data = eeg_checkset( cond_EEG.data );%add pulling out condition specific data + cond_EEG = pop_selectevent(EEG, 'Condition',(Conds{cond}), 'deleteevents','on','deleteepochs','on','invertepochs','off'); + cond_data = eeg_checkset(cond_EEG.data);%add pulling out condition specific data + fprintf('\n* Condition %s data shape: %d\n *\n', Conds{cond}, size(cond_data, 3)); %Create TrialNums structure to be saved out TrialNums( sub+(cond-1)+((sub-1)*(length(Conds)-1)) ).subject = subject; diff --git a/code/preprocessing-eeg/transform_csd.m b/code/preprocessing-eeg/transform_csd.m index 648643a..c040c3f 100644 --- a/code/preprocessing-eeg/transform_csd.m +++ b/code/preprocessing-eeg/transform_csd.m @@ -5,11 +5,17 @@ % Workshop on 02/22. This script uses parts of the "set up" structure from % the MADE preprocessing pipeline (Debnath, Buzzell, et. al., 2020) -clear % clear matlab workspace -clc % clear matlab command window - +%clear % clear matlab workspace +%clc % clear matlab command window +cluster = parcluster('local'); %% Setting up other things - +parpool(cluster, str2num(getenv('SLURM_CPUS_PER_TASK'))) % this should be same as --cpus-per-task +pool = gcp('nocreate'); % Get the current parallel pool without creating a new one +if isempty(pool) + disp('No parallel pool is currently running.'); +else + disp(['Parallel pool with ', num2str(pool.NumWorkers), ' workers is running.']); +end %Location of MADE and ADJUSTED-ADJUST scripts % addpath(genpath([main_dir filesep 'MADE-EEG-preprocessing-pipeline']));% enter the path of the EEGLAB folder in this line addpath(genpath('/home/data/NDClab/tools/lab-devOps/scripts/MADE_pipeline_standard/eeg_preprocessing'));% enter the path of the folder in this line @@ -71,35 +77,41 @@ trodes{site} = num2str(site); end Montage_64=ExtractMontage('/home/data/NDClab/analyses/thrive-theta-ddm/code/preprocessing-eeg/64ch_bv_montage_csd.csd', trodes'); -% MapMontage(Montage_64); [G, H] = GetGH(Montage_64); diary(sprintf('csd_log_%s.log', datestr(now, 'mm-dd-yyyy_HH_MM_SS'))) % loop through each participant in the study -for subject = 1:length(datafile_names) - - % extract participant number - subNumText = datafile_names{subject}(5:11); - - %load the original data set - EEG = pop_loadset('filename', datafile_names{subject}, 'filepath', datafile_paths{subject}); - EEG = eeg_checkset(EEG); - disp(datafile_names{subject}) - disp(save_location) - %remove all the non-stim-locking markers (should have done already...) - EEG = pop_selectevent(EEG, 'latency','-.1 <= .1','deleteevents','on'); - EEG = eeg_checkset(EEG); - - for ne = 1:length(EEG.epoch) - myEEG = single(EEG.data(:, :, ne)); - MyResults = CSD(myEEG, G, H); % compute CSD for <channels-by-samples> 2-D epoch - data(:, :, ne) = MyResults; +parfor subject = 1:length(datafile_names) + try + % extract participant number + subNumText = datafile_names{subject}(5:11); + + %load the original data set + EEG = pop_loadset('filename', datafile_names{subject}, 'filepath', datafile_paths{subject}); + EEG = eeg_checkset(EEG); + + %remove all the non-stim-locking markers (should have done already...) + EEG = pop_selectevent(EEG, 'latency','-.1 <= .1','deleteevents','on'); + EEG = eeg_checkset(EEG); + fprintf('Subject %s: Processing %d events\n', subNumText, length(EEG.event)); + data = zeros(size(EEG.data), 'single'); + for ne = 1:length(EEG.epoch) + myEEG = single(EEG.data(:, :, ne)); + MyResults = CSD(myEEG, G, H); % compute CSD for <channels-by-samples> 2-D epoch + data(:, :, ne) = MyResults; + end + + EEG.data = data; + EEG = eeg_checkset(EEG); + disp(size(EEG.event)) + + EEG = pop_editset(EEG, 'setname', datafile_names{subject}); + EEG = pop_saveset(EEG, 'filename', datafile_names{subject}, 'filepath', save_location); + fprintf('Subject %s: Successfully processed and saved\n', subNumText); +% clear data; + catch ME + fprintf('Error processing subject %s: %s\n', subNumText, ME.message); + continue; end - EEG.data = data; - - data(:,:,:) = NaN; - - EEG = pop_editset(EEG, 'setname', datafile_names{subject}); - EEG = pop_saveset(EEG, 'filename', datafile_names{subject}, 'filepath', save_location); end