Contents

%%%Om #!/usr/local/bin/matlab%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to run ICA and prune Data using SASICA
% 1. RUN ICA (binica) for each participant and save dataset
% 2. Load the ICAed dataset, run SASICA on it to prune it from artefactual
%    components. This can be done manually, but SASICA does a fairly good
%    job automatically, provided the parameters are selected properly.
%    % We've chosen auto-correlation and ADJUST algorithm based rejections
%    % For reasonably good data (namely, adult EEG data), this may an
%    % overkill; use EOG channel correlations in that case (we didn't have
%    % that choice here!)
% 3. NOW, and here's the crux of the thing, load THE_MAIN_PREICA_DATASET,
%    which is the one with the bad channels removed, but without the
%    pre-ICA manual rejections (from the last step).
% 4. Copy the ICA related things (wins, spheres, weights..) from the ICAed
%    and SASICAed dataset to THE_MAIN_PREICA_DATASET. But not the
%    activations.  These should be calculated (see below).
% 5. Copy the bad components list generated by SASICA/ADJUST from the ICAed
%    and SASICAed dataset to THE_MAIN_PREICA_DATASET.
% 6. Now, forget the ICAed and SASICAed dataset.
% 7. Prune data: Reject bad components from THE_MAIN_PREICA_DATASET.
% 8. SAVE this dataset separately, and check (with channel scroll) whether the blink
%    artefacts you had in THE_MAIN_PREICA_DATASET are gone!  They must have!
% 9. This is THE_ICA_PRUNED_DATASET.
%
% With valuable advice from: Alessandro Tavano
%
% Author: R. Muralikrishnan -- 24.02.2017; Last modified: 07.03.2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%m
clear all;
close all;

eeglab;

V_Folder_EEGLabData = '/home/muralir/Documents/AB/Experiments/NeuroFemi/EEGLab-Data/';

V_Participant = 'xx';
% CA_x => CellArray
CA_Participants = {'NF01','NF02','NF03','NF05','NF06','NF07','NF09','NF10','NF11','NF12','NF13','NF14','NF15','NF16','NF20','NF22','NF23','NF24','NF25','NF27','NF28','NF30','NF32','NF33','NF34','NF37','NF38','NF39','NF42','NF43','NF45','NF46','NF47','NF48','NF49','NF51','NF52','NF53','NF54','NF55','NF56'};

Participants who are ruled out a priori, and for whom NO ICA needs to be done

% NF41, NF36 : Incomplete/non-existent data
% NF29, NF44, NF50 : Left-handed/handedness unclear
% NF04, NF08, NF17, NF18, NF19, NF26, NF31, NF35, NF40 : Too few correct responses. See Avr-Rej-Details-Filetered.txt (in EEG-Data folder).
% NF21      ADJUST SASICA didn't work for some reason!

% Remaining participants 56-15 = 41.

RUN ICA

CA_Participants = {'NF56'};
for V_Counter = 1:length(CA_Participants)
    V_Participant = CA_Participants{V_Counter};

    V_DataFile = fullfile(V_Folder_EEGLabData, V_Participant, [V_Participant, '_bp03_40r_ica_in.set']);
    EEG = pop_loadset('filename', V_DataFile);


    % Change to the participant's folder, so the binica.* files will be stored there
    V_Command = sprintf('../%s', V_Participant);
    cd(V_Command);

    % ICA

    %N_of_Components = rank(EEG.data(:,:));
    N_of_Components = 64;

    % Run binICA -- This uses ica_linux through binica.m
    [EEG.icaweights, EEG.icasphere] = binica(EEG.data, 'extended', 1, 'pca', N_of_Components, 'filenum', 0);
    % Without the binary binica installed, the same can be done as in Katrin's
    % script as follows:
    % pop_runica(EEG,'icatype','binica','options',{'extended', 1, 'pca', N_of_Channels_Remaining, 'filenum' 1});


    % Calculate ICA Activations -- http://sccn.ucsd.edu/wiki/Binica
    % After this, 'EEGLab redraw' should show ICA Weights: Yes for the dataset
    EEG.icaact = EEG.icaweights*EEG.icasphere*EEG.data;  % Matrix multiplication

Change back to the scripts folder

    V_Command = sprintf('../AB-Scripts');
    cd(V_Command);

    % Save the resultant file as an EEGLab .set

    EEG.filename = [V_Participant '_bp03_40r_ica_out'];
    EEG.filepath = [V_Folder_EEGLabData, V_Participant];
    EEG.comments = char([V_Participant, '_bp03_40r_ica_out'], ['; Sampling Rate: 500 Hz; BPF: 0.3-40 Hz; ReRef: Linked Mastoids; With ICA Weights']);
    pop_saveset(EEG, 'filename', EEG.filename, 'filepath', EEG.filepath);

    fprintf('########################################################\n\n');
    fprintf('Saved ICA data for ... %s.\n\n', V_Participant);

  %Om
end;

PRUNE DATA

Since non-repetitive and out of the ordinary artefacts have been rejected before ICA as part of the ICA Clean-up, the latencies for events in the resultant ICAed data above have changed! We actually want the triggers to be in place. So removing components from this data won't do for further analysis. Instead, we do the following: -

% and run SASICA on the unrejected file that was used as input for
% ICA clean-up.  So we just copy the ICA weights etc. to that file from the
% ICA output file.

Run SASICA, based on which remove components from the main pre-ICA clean-up file

%CA_Participants = {'NF22','NF23','NF24','NF25','NF27','NF28','NF30','NF32','NF33','NF34','NF37','NF38','NF39','NF42','NF43','NF45','NF46','NF47','NF48','NF49','NF51','NF52','NF53','NF54','NF55','NF56'};

%CA_Participants = {'NF21'}; % The Adjust SASICA doesn't work for NF21!
V_Participants_Ignored = 0;
for V_Counter = 1:length(CA_Participants)

    V_Participant = CA_Participants{V_Counter};

    V_DataFile = fullfile(V_Folder_EEGLabData, V_Participant, [V_Participant, '_bp03_40r_ica_out.set']);
    EEG_ICA = pop_loadset('filename', V_DataFile);

    % Run SASICA and identify ICA components to remove.
    % (We have 64 ICA components in all...see N_of_Components above!).
    % This step will generate EEG.reject.gcompreject, which would be a vector
    % containing 1's or 0's for each component depending upon whether the
    % component is marked for removal by SASICA or not!
    % This is the vector used in pop_subcomp (when it is used with GUI)
    % to prompt for whether we want these to be removed.


    %     % If one or more of the VEOG / HEOG channels are not found in the
    %     % dataset of the current participant, abort running SASICA on this
    %     % dataset.
    %
    %     % http://stackoverflow.com/questions/14480876/using-find-with-a-struct
    %     % VEOGTopPresent = ~isempty( find( cellfun(@(x)strcmp(x,'E124'),{EEG_ICA.chanlocs.labels}) ) );
    %     VEOGT_Present = ~isempty( find( cellfun(@(x)strcmp(x,'E9'),{EEG_ICA.chanlocs.labels}), 1 ) );
    %     VEOGB_Present = ~isempty( find( cellfun(@(x)strcmp(x,'E127'),{EEG_ICA.chanlocs.labels}), 1 ) );
    %     HEOGL_Present = ~isempty( find( cellfun(@(x)strcmp(x,'E22'),{EEG_ICA.chanlocs.labels}), 1 ) );
    %     HEOGR_Present = ~isempty( find( cellfun(@(x)strcmp(x,'E9'),{EEG_ICA.chanlocs.labels}), 1 ) );
    %
    %
    %
    %
    %     if ~(VEOGT_Present && VEOGB_Present && HEOGL_Present && HEOGR_Present)
    %         fprintf('One or more of the EOG Channels are not present in the dataset!\n\n');
    %         fprintf('Aborting SASICA for Participant %s. \n\n', V_Participant);
    %         V_Participants_Ignored = V_Participants_Ignored + 1;
    %
    %         continue;
    %         % Don't do SASICA...go to the next for-loop iteration...next participant!!!
    %         %http://stackoverflow.com/questions/1082605/jump-command-in-matlab
    %
    %     end;
    %
    % SASICA with VEOG E14, E126, and HEOG E21 and E14...the thing will run if at least one of the specified EOG electrodes is found in the dataset.
    % EEG_ICA = eeg_SASICA(EEG_ICA,'MARA_enable',0,'FASTER_enable',0,'FASTER_blinkchanname','No channel','ADJUST_enable',0,'chancorr_enable',0,'chancorr_channames','No channel','chancorr_corthresh','auto 4','EOGcorr_enable',1,'EOGcorr_Heogchannames','E14 E21','EOGcorr_corthreshH','auto 4','EOGcorr_Veogchannames','E14 E126','EOGcorr_corthreshV','auto 4','resvar_enable',0,'resvar_thresh',15,'SNR_enable',0,'SNR_snrcut',1,'SNR_snrBL',[-Inf 0] ,'SNR_snrPOI',[0 Inf] ,'trialfoc_enable',0,'trialfoc_focaltrialout','auto','focalcomp_enable',1,'focalcomp_focalICAout','auto','autocorr_enable',1,'autocorr_autocorrint',20,'autocorr_dropautocorr','auto','opts_noplot',1,'opts_nocompute',0,'opts_FontSize',14);

    % SASICA with frontal channel correlations...the thing will run if at...the thing will run if at least one of the specified frontal electrodes is found in the dataset.
    %EEG_ICA = eeg_SASICA(EEG_ICA,'MARA_enable',0,'FASTER_enable',0,'FASTER_blinkchanname','No channel','ADJUST_enable',0,'chancorr_enable',1,'chancorr_channames','E21 E127 E14 E126 E22 E15 E9','chancorr_corthresh','auto 4','EOGcorr_enable',0,'EOGcorr_Heogchannames','No channel','EOGcorr_corthreshH','auto 4','EOGcorr_Veogchannames','No channel','EOGcorr_corthreshV','auto 4','resvar_enable',0,'resvar_thresh',15,'SNR_enable',0,'SNR_snrcut',1,'SNR_snrBL',[-Inf 0] ,'SNR_snrPOI',[0 Inf] ,'trialfoc_enable',0,'trialfoc_focaltrialout','auto','focalcomp_enable',1,'focalcomp_focalICAout','auto','autocorr_enable',1,'autocorr_autocorrint',20,'autocorr_dropautocorr','auto','opts_noplot',1,'opts_nocompute',0,'opts_FontSize',14);

    % SASICA with ADJUST enabled; no EOG or other channels specified
    EEG_ICA = eeg_SASICA(EEG_ICA,'MARA_enable',0,'FASTER_enable',0,'FASTER_blinkchanname','No channel','ADJUST_enable',1,'chancorr_enable',0,'chancorr_channames','No channel','chancorr_corthresh','auto 4','EOGcorr_enable',0,'EOGcorr_Heogchannames','No channel','EOGcorr_corthreshH','auto 4','EOGcorr_Veogchannames','No channel','EOGcorr_corthreshV','auto 4','resvar_enable',0,'resvar_thresh',15,'SNR_enable',0,'SNR_snrcut',1,'SNR_snrBL',[-Inf 0] ,'SNR_snrPOI',[0 Inf] ,'trialfoc_enable',0,'trialfoc_focaltrialout','auto','focalcomp_enable',1,'focalcomp_focalICAout','auto','autocorr_enable',1,'autocorr_autocorrint',20,'autocorr_dropautocorr','auto','opts_noplot',1,'opts_nocompute',0,'opts_FontSize',14);


    %EEG_ICA = SASICA();
    %waitfor(gcf);

    % Now, instead of removing the components from the ICAed data,
    % we need to remove them from the main data (i.e., the data that was
    % not used for pre-ICA clean-up.  This is because, due to the clean-up,
    % we have lost considerable amount of the original data (including trials).  Further, and as a consequence, the
    % events latencies in the resultant data have changed.  That's of no use for us.

    % So, we just copy the ICA weights and spheres from the ICAed data,
    % as well as the SASICA generated EEG.reject.gcompreject from above
    % to the main data (i.e., data used for pre-ICA clean-up).

    % Load the main data to be used for further analysis
    % This is the data with the bad channels removed, but in a state prior
    % to the pre-ICA clean-up.

    V_DataFile = fullfile(V_Folder_EEGLabData, V_Participant, [V_Participant, '_bp03_40r_ica_bc.set']);
    EEG = pop_loadset('filename', V_DataFile);

    % Copy the ICA weights and spheres and SASICA rejections to the main data set

    % EEG.icaact should not be copied from the ICA dataset to the main one,
    % for the number of data points is different (because of pre-ICA
    % cleanup). Instead, copy the weights and spheres, and calculate the
    % activations for the main data.

    EEG.icawinv = EEG_ICA.icawinv;
    EEG.icasphere = EEG_ICA.icasphere;
    EEG.icaweights = EEG_ICA.icaweights;
    EEG.icachansind = EEG_ICA.icachansind;

    EEG.icaact = EEG.icaweights*EEG.icasphere*EEG.data;  % Matrix multiplication

    EEG.reject.gcompreject = EEG_ICA.reject.gcompreject;
    EEG.reject.SASICA = EEG_ICA.reject.SASICA;

    % Extract the Component numbers to remove from the output of
    % SASICA in EEG.reject.gcompreject.  Wherever this is 1, that
    % particular ICA component needs to be removed from the data.

    V_Components_to_remove = find(EEG.reject.gcompreject == 1);

    % Remove the components from the main data without further ado!

    EEG = pop_subcomp(EEG, V_Components_to_remove);
    % EEG = pop_subcomp(EEG, V_Components_to_remove, 1); % To ask for
    % confirmation, and visualise the difference that the removals would
    % make to the data.

    % Save the resultant file as an EEGLab .set

    EEG.filename = [V_Participant '_bp03_40r_ica_pruned'];
    EEG.filepath = [V_Folder_EEGLabData, V_Participant];
    EEG.comments = char([V_Participant, '_bp03_40r_ica_pruned'], ['; Sampling Rate: 500 Hz; BPF: 0.3-40 Hz; ReRef: Linked Mastoids; SASICA based ICA Pruned Original Data']);
    pop_saveset(EEG, 'filename', EEG.filename, 'filepath', EEG.filepath);

    fprintf('^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n\n');
    fprintf('Saved ICA Pruned data for ... %s.\n\n', V_Participant);
    fprintf('Number of Removed Components ... %d. \n\n', length(V_Components_to_remove));
    fprintf('########################################################\n\n');

%  Om



end;

fprintf('Number of participants for which SASICA was not performed: %d. \n\n', V_Participants_Ignored);

The dataset ABxx_bp03_40r_ica_pruned.set is the one that will be used in the next step!

%# See EEGLab_WriteCntData_Interpolate_PostICA.m

%  Om
























%
%
% for i = 1%:length(A_Participants)
%   V_Participant = CA_Participants{i};
%   V_FileIn = fullfile(V_FolderIn, [V_Participant '_bp03_100r.set']);
%   fprintf('########################################################\n\n');
%   fprintf('Reading %s.\n\n', V_FileIn);
%   EEG = pop_loadset(V_FileIn);
%   %EEG = pop_repchan(EEG, [1,4]);
%   EEG = pop_runica(EEG,'runica','pca', 32);
% EEG.filepath = [V_FolderOut, V_Participant, '/R'];
%
% %ncomp = rank(EEG.data(:,:));
% %EEG = pop_runica(EEG,'icatype','binica','options',{'extended' 1 'pca' ncomp, 'filenum' 1});
%
% %EEG = binica(EEG, 'options',{'extended' 1, 'pca' ncomp, 'filenum' 1});
%
%   pop_saveset(EEG, 'filename', fullfile([V_Participant '__bp03_100r_ICA']));
% end