script_09_ica_clustering.m

% Wakeman & Henson Data analysis: ICA clustering.
%
% Authors: Arnaud Delorme, Ramon Martinez-Cancino, Johanna Wagner, Romain Grandchamp

% Clearing all is recommended to avoid variable not being erased between calls
clear;

% Path to data below. Using relative paths so no need to update.
% This script loads the epoched data generated by Session 2
path2data = fullfile(pwd,'ds000117_pruned', 'derivatives', 'meg_derivatives', 'sub-01', 'ses-meg/', 'meg/'); % Path to data
filename_famous = 'wh_S01_run_01_ERP_Analysis_Session_2_famous_out.set';
filename_unfamiliar = 'wh_S01_run_01_ERP_Analysis_Session_2_unfamiliar_out.set';
filename_scrambled = 'wh_S01_run_01_ERP_Analysis_Session_2_scrambled_out.set';

% Start EEGLAB
[ALLEEG, EEG, CURRENTSET] = eeglab;

%% Loading epoched data from Session 2
% Load the preprocessed data first (needed for ICA component info)
EEG = pop_loadset('filename', 'wh_S01_run_01_preprocessing_data_session_1_out.set','filepath',path2data);
[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1);

% Load the three epoched datasets from Session 2
ALLEEG(2) = pop_loadset('filename', filename_famous,'filepath',path2data);
ALLEEG(3) = pop_loadset('filename', filename_unfamiliar,'filepath',path2data);
ALLEEG(4) = pop_loadset('filename', filename_scrambled,'filepath',path2data);

% Copy ICA information from the continuous dataset to epoched datasets if not present
if ~isfield(ALLEEG(2), 'icaweights') && isfield(ALLEEG(1), 'icaweights')
    for i = 2:4
        ALLEEG(i).icasphere = ALLEEG(1).icasphere;
        ALLEEG(i).icaweights = ALLEEG(1).icaweights;
        ALLEEG(i).icawinv = ALLEEG(1).icawinv;
        ALLEEG(i).icachansind = ALLEEG(1).icachansind;
        if isfield(ALLEEG(1), 'etc') && isfield(ALLEEG(1).etc, 'ic_classification')
            ALLEEG(i).etc.ic_classification = ALLEEG(1).etc.ic_classification;
        end
        % Recalculate ICA activations for epoched data
        ALLEEG(i) = eeg_checkset(ALLEEG(i), 'ica');
    end
end

%% ----------------------
%% BELOW IS PLOTTING ONLY
%% ----------------------

%% plot ERP scalp distribution 
figure; pop_timtopo(ALLEEG(2), [-100  600], [NaN], 'ERP data and scalp maps of Famous Epoched');
figure; pop_timtopo(ALLEEG(3), [-100  600], [NaN], 'ERP data and scalp maps of Unfamiliar Epoched');
figure; pop_timtopo(ALLEEG(4), [-100  600], [NaN], 'ERP data and scalp maps of Scrambled Epoched');

%% plot Component 5 contributions to ERP scalp envelope with channel projection 
pop_envtopo(ALLEEG(2), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'compnums',[5],'title', 'Famous','electrodes','off', 'plotproj', 1);

%% plot Component contributions to ERP scalp envelope 
pop_envtopo(ALLEEG(2), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'title', 'Largest ERP components of Famous Epoched','electrodes','off');
pop_envtopo(ALLEEG(3), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'title', 'Largest ERP components of Unfamiliar Epoched','electrodes','off');
pop_envtopo(ALLEEG(4), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'title', 'Largest ERP components of Scrambled Epoched','electrodes','off');

%% plot Component contributions to ERP scalp envelope - removing artefact ICs
pop_envtopo(ALLEEG(2), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Largest ERP components of Famous Epoched','electrodes','off');
pop_envtopo(ALLEEG(3), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Largest ERP components of Unfamiliar Epoched','electrodes','off');
pop_envtopo(ALLEEG(4), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Largest ERP components of Scrambled Epoched','electrodes','off');

%% plot Component contributions to ERP scalp envelope - removing artefact ICs - plot differences between conditions
pop_envtopo(([ALLEEG(2)  ALLEEG(4)]), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Famous - Scrambled','electrodes','off');
pop_envtopo(([ALLEEG(2)  ALLEEG(3)]), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Famous - Unfamiliar','electrodes','off');
pop_envtopo(([ALLEEG(3)  ALLEEG(4)]), [-100  600] ,'limcontrib',[0 400],'compsplot',[7],'subcomps',[1 6],'title', 'Unfamiliar - Scrambled','electrodes','off');

%% Identify Brain ICs using IC Label classification results
if isfield(ALLEEG(1).etc, 'ic_classification')
    [M,I] = max(ALLEEG(1).etc.ic_classification.ICLabel.classifications,[],2);                       % Use max prob for classification
    Brain_comps = find(I == find(strcmp(ALLEEG(1).etc.ic_classification.ICLabel.classes, 'Brain')));
    
    %% Subtract artefactual components from the EEG
    ALLEEG(2) = pop_subcomp( ALLEEG(2), Brain_comps, 0, 1);
    ALLEEG(3) = pop_subcomp( ALLEEG(3), Brain_comps, 0, 1);
    ALLEEG(4) = pop_subcomp( ALLEEG(4), Brain_comps, 0, 1);
end

%% Rename datasets
ALLEEG(2) = pop_editset(ALLEEG(2), 'setname', 'Famous', 'run', []);
ALLEEG(3) = pop_editset(ALLEEG(3), 'setname', 'Unfamiliar', 'run', []);
ALLEEG(4) = pop_editset(ALLEEG(4), 'setname', 'Scrambled', 'run', []);

%% plot ERP scalp distribution 
figure; pop_timtopo(ALLEEG(2), [-100  600], [NaN], 'Famous');
figure; pop_timtopo(ALLEEG(3), [-100  600], [NaN], 'Unfamiliar');
figure; pop_timtopo(ALLEEG(4), [-100  600], [NaN], 'Scrambled');

%% plot ERP scalp distribution at each ERP peak
figure; pop_timtopo(ALLEEG(2), [-100  600], [120  170  250], 'Famous');
figure; pop_timtopo(ALLEEG(3), [-100  600], [120  170  250], 'Unfamiliar');
figure; pop_timtopo(ALLEEG(4), [-100  600], [120  170  250], 'Scrambled');

%% plot 3 largest contributing ICs to ERP
pop_envtopo(ALLEEG(2), [-100  600] ,'limcontrib',[0 400],'compsplot',[3],'title', 'Famous','electrodes','off');
pop_envtopo(ALLEEG(3), [-100  600] ,'limcontrib',[0 400],'compsplot',[3],'title', 'Unfamiliar','electrodes','off');
pop_envtopo(ALLEEG(4), [-100  600] ,'limcontrib',[0 400],'compsplot',[3],'title', 'Scrambled','electrodes','off');

%% Visualize channel ERPs in 2D
pop_topoplot(ALLEEG(2), 1, [25:25:300] ,'Famous',[3 4] ,0,'electrodes','on');
pop_topoplot(ALLEEG(3), 1, [25:25:300] ,'Unfamiliar',[3 4] ,0,'electrodes','on');
pop_topoplot(ALLEEG(4), 1, [25:25:300] ,'Scrambled',[3 4] ,0,'electrodes','on');

%% Plot channel ERPs in topographic array
figure; pop_plottopo(ALLEEG(2), [1:EEG.nbchan] , 'Famous', 0, 'ydir',1);
figure; pop_plottopo(ALLEEG(3), [1:EEG.nbchan] , 'Unfamiliar', 0, 'ydir',1);
figure; pop_plottopo(ALLEEG(4), [1:EEG.nbchan] , 'Scrambled', 0, 'ydir',1);


%% plot average ERPs for each condition with standard deviation

% find channel index of eeg065
Chanind = find(strcmp({ALLEEG(2).chanlocs.labels},'EEG065'));
if isempty(Chanind)
    Chanind = 1;
end

% create timevector for plotting 
[val, indL] = min(abs(ALLEEG(2).times+200)); %get timepoints for -200 and 800 Latencies
[val, indU] = min(abs(ALLEEG(2).times-800));
timevec = ALLEEG(2).times(indL:indU); % create timevector

% create datavectors for plotting each condition
av_datavecF = mean(ALLEEG(2).data(Chanind,indL:indU,:),3); % average
std_datavecF = std(ALLEEG(2).data(Chanind,indL:indU,:),1,3); % standard deviation

figure;
X2 = [[timevec],fliplr([timevec])];                %#create continuous x value array for plotting
Y2 = [av_datavecF-std_datavecF,fliplr(av_datavecF+std_datavecF)];              %#create y values for out and then back
fill(X2,Y2,[153/255 204/255 255/255]);
hold on
plot(timevec,av_datavecF, 'b', 'LineWidth',2)
xline(0, 'LineWidth',2)
yline(0, 'LineWidth',2)
xlabel('Latency ms')
ylabel('mu Volt')
title([ 'famous faces channel ' EEG.chanlocs(Chanind).labels ]);
set(gca, 'FontSize', 15)

av_datavecU = mean(ALLEEG(3).data(Chanind,indL:indU,:),3); % average
std_datavecU = std(ALLEEG(3).data(Chanind,indL:indU,:),1,3); % standard deviation

figure;
X2 = [[timevec],fliplr([timevec])];                %#create continuous x value array for plotting
Y2 = [av_datavecU-std_datavecU,fliplr(av_datavecU+std_datavecU)];              %#create y values for out and then back
fill(X2,Y2,[153/255 204/255 255/255]);
hold on
plot(timevec,av_datavecU, 'b', 'LineWidth',2)
xline(0, 'LineWidth',2)
yline(0, 'LineWidth',2)
xlabel('Latency ms')
ylabel('mu Volt')
title([ 'unfamiliar faces channel ' EEG.chanlocs(Chanind).labels ]);
set(gca, 'FontSize', 15)

av_datavecS = mean(ALLEEG(4).data(Chanind,indL:indU,:),3); % average
std_datavecS = std(ALLEEG(4).data(Chanind,indL:indU,:),1,3); % standard deviation

figure;
X2 = [[timevec],fliplr([timevec])];                %#create continuous x value array for plotting
Y2 = [av_datavecS-std_datavecS,fliplr(av_datavecS+std_datavecS)];              %#create y values for out and then back
fill(X2,Y2,[153/255 204/255 255/255]);
hold on
plot(timevec,av_datavecS, 'b', 'LineWidth',2)
xline(0, 'LineWidth',2)
yline(0, 'LineWidth',2)
xlabel('Latency ms')
ylabel('mu Volt')
title([ 'scrambled faces channel ' EEG.chanlocs(Chanind).labels ]);
set(gca, 'FontSize', 15)

%% plot superimposed ERPs

figure;plot(timevec,av_datavecF, 'LineWidth',2, 'color', 'r'); hold on
plot(timevec,av_datavecU, 'LineWidth',2, 'color', 'b')
plot(timevec,av_datavecS, 'LineWidth',2, 'color', 'g')
fillcurves(timevec,av_datavecF-std_datavecF,av_datavecF+std_datavecF, 'r', 0.2);
fillcurves(timevec,av_datavecU-std_datavecU,av_datavecU+std_datavecU, 'b', 0.2);
fillcurves(timevec,av_datavecS-std_datavecS,av_datavecS+std_datavecS, 'g', 0.2);
xline(0, 'LineWidth',2)
yline(0, 'LineWidth',2)
xlabel('Latency ms')
ylabel('mu Volt')
legend('famous', 'unfamiliar', 'scrambled')
set(gca, 'FontSize', 15)
title([ 'face types for channel ' EEG.chanlocs(Chanind).labels ]);

%% ERPimage
figure; pop_erpimage(ALLEEG(2),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{},[],'' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );
figure; pop_erpimage(ALLEEG(3),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{},[],'' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );
figure; pop_erpimage(ALLEEG(4),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{},[],'' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );

% sort by event latency (when there are enough of left_nonsym and right_sym events
% figure; pop_erpimage(ALLEEG(2),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{ 'left_nonsym' 'right_sym'},[],'latency' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );
% figure; pop_erpimage(ALLEEG(3),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{ 'left_nonsym' 'right_sym'},[],'latency' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );
% figure; pop_erpimage(ALLEEG(4),1, [Chanind],[[]],EEG.chanlocs(Chanind).labels,3,1,{ 'left_nonsym' 'right_sym'},[],'latency' ,'yerplabel','\muV','erp','on','limits',[-100 1200 NaN NaN NaN NaN NaN NaN] ,'cbar','on','topo', { [Chanind] EEG.chanlocs EEG.chaninfo } );