1 % Run independent components analysis (ICA) on MEG data,
   2 % classify components related to blinks, eye-movements, 
   3 % and pulse artefacts, then remove artefact component(s) 
   4 % related to, e.g., blinks.
   5 %
   6 % This script assumes you've run the session 2 script
   7 %   cbu_meeg_masterclass2
   8 % or at least that you have a directory for the project
   9 % directory in your imaging space.
  10 %
  11 
  12 %--Set Parameters------------------------------------------
  13 
  14 % Working directory 
  15 wd = '/home/ms02/MEGMasterclass/Marie/'; %% !!!CHANGE TO YOUR OWN !!!
  16 cd(wd)
  17 
  18 % Extra bit for MEG masterclass - copy ICAed files from here:
  19 %ica_sourcedir = '/imaging/jt03/public/demo/';
  20 ica_sourcedir = '/imaging/jt03/public/demo/icadata/';
  21 
  22 % Add relevant directories to your path:
  23 addpath /imaging/local/spm/spm5/cbu_updates/
  24 addpath /imaging/local/meg_misc/  
  25 addpath /imaging/local/spm_eeglab/  % ica + other EEGLAB functions
  26 addpath /imaging/jt03/public/demo/  % functions not yet uploaded to meg_misc!
  27 
  28 % Processing flags & variables for ICA
  29 do_runica    = 0; % do ICA (takes time!)?
  30 do_copyica   = 1; % just copy ICAed files instead!
  31 mem_flag     = 0; % flag if run into memory problems!
  32 
  33 icaflag      = [1 1 1 0];  % don't run on grms
  34 ica_Npc      = 32;  % Number of ICs to extract (PCA run first)
  35 ica_samppct  = 1;  % Proportion of data to give to ICA
  36 ica_class    = {'blink','horiz','pulse'}; % Artefact(s) to classify
  37 ica_chans    = {'veog','heog','ecg'};     % Reference channels for correlations
  38 ica_rem      = {'blink'};                 % Artefact(s) to remove 
  39 
  40 %%%%%%%%%% Different types of sensor
  41 typ          = {'mags' 'grds' 'eeg' 'grms'};
  42 Ntyp         = length(typ);
  43 eegflag      = [0 0 1 0];  % Binary flag for which typ is EEG
  44 
  45 %%%%%%%%%%% subject-specific parameters
  46 dosubs       = 1:8;
  47 Nsub         = length(dosubs);
  48 
  49 
  50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  51 %--COPY ICAed DATA-----------------------------------------
  52 % Extra bit for MEG masterclass. Instead of running ICA on 
  53 % all 8 subjects, just copy ICAed files:
  54 
  55 if do_copyica
  56 	for sub=1:Nsub
  57 		swd = sprintf('%ssub%d',wd,sub);
  58 		try   cd(swd)
  59     catch mkdir(swd)
  60           cd(swd)
  61     end
  62     
  63 		ses=1;
  64 	
  65 		% Assign filenames:
  66 		basefile=sprintf('dsub%d_ses%d_trans',sub,ses);
  67 		icafile=sprintf('ica_%s',basefile);
  68 		[p fstem ext]=fileparts(basefile);
  69 		
  70 		for typn = 1:Ntyp
  71       if icaflag(typn)
  72         % Copy ica*.mat file to local directory:
  73         fn=sprintf('%s-%s.mat',icafile,typ{typn});
  74         if exist(fn)~=2
  75           disp(sprintf('Copying file %s...',fn))
  76           [cpfail tmp]=unix(sprintf('! cp %s/%s ./',ica_sourcedir,fn));
  77           if cpfail
  78             error('Copy failed!')
  79           end
  80           disp(sprintf('Done.'))
  81         end
  82         % Ensure .dat file is present (copy over if not):
  83         D=spm_eeg_ldata(fn);
  84         D.path=swd;
  85         D.fnamedat=sprintf('%s-%s.dat',fstem,typ{typn});
  86         if exist(D.fnamedat)~=2
  87           disp(sprintf('Copying .dat file ...'))
  88           [cpfail tmp]=unix(sprintf('! cp %s/%s ./',ica_sourcedir,D.fnamedat));
  89           if cpfail
  90             error('Copy failed!')
  91           end
  92           disp(sprintf('Done.'))
  93         end
  94         % Remove any previous classifications:
  95         D.ica.class=[];
  96         save(D.fname,'D');
  97 
  98         clear D
  99       end % icaflag
 100 		end % typn
 101 	end % sub
 102 end % do_copyica
 103 
 104 
 105 % End of extra bit
 106 %----------------------------------------------------------
 107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 108 
 109 
 110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 111 %--Processing begins here!---------------------------------
 112 
 113 for sub = 1:Nsub;
 114 	subnum    = dosubs(sub);
 115 	Nses(sub)=1;
 116 
 117 	cd(wd)
 118 	swd = sprintf('sub%d',subnum);
 119 	try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
 120 
 121 	for ses = 1:Nses(sub)
 122 
 123 		% We are using the files that have been imported to SPM,
 124 		% downsampled, and split (separated into mags and grds)
 125 		% Thus, basefile for subject 1 is 'dsub1_ses1_trans'
 126 		basefile=sprintf('dsub%d_ses%d_trans',sub,ses);
 127 
 128 		for typn = 1:Ntyp
 129 
 130 
 131 %%%%% Artefact compensation using ICA
 132 
 133 			if icaflag(typn)
 134 
 135 %--Run ICA--------------------------------------------------
 136 %
 137 % ICA takes somewhere between 5-10 minutes on each file,
 138 %  longer if pca=65 or if using non-split mags+grds file.
 139 %  I would recommend running this bit of code in a batch over
 140 %  subjects, say overnight, then running the next bit 
 141 %  (classify and remove artefacts) in the morning!
 142 
 143 			if do_runica % If...else...end is Extra bit for MEG masterclass
 144 	
 145 				clear S; S.D = sprintf('%s-%s.mat',basefile,typ{typn});
 146 				S.samppct = ica_samppct;
 147 				S.args    = {'extended',1,...
 148 										 'maxsteps',800,...
 149 										 'pca',ica_Npc};
 150 				if ~mem_flag
 151 					D = spm_eeglab_runica(S);
 152 				elseif mem_flag
 153 					save('tmp.mat','S');
 154 					! matlab -nodesktop -nosplash -r "spm eeg; load tmp.mat; spm_eeglab_runica(S); exit"
 155 					try delete('tmp.mat'); end
 156 					D.fname=['ica_' S.D];
 157 				end %mem_flag
 158 				
 159 			else % Extra bit for MEG masterclass
 160 					
 161 					D.fname=sprintf('ica_%s-%s.mat',basefile,typ{typn});
 162 					
 163 			end % Extra bit for MEG masterclass
 164 
 165 					
 166 %--Classify and Remove Artefact ICs-------------------------
 167 %
 168 % View raw data on blink- and pulse-sensitive channels,
 169 %  compute IC activations, correlate with physiological 
 170 %  signals, plot corr Z-scores, plot topographies,
 171 %  classify ICs as blink, horizontal eye mov'ts, pulse,
 172 %  and remove selected artefact components.
 173 
 174 				clear S; S.D = D.fname;
 175 				S.mode     = {'both'};  % classify + remove
 176 				S.artlabs  = ica_class;
 177 				S.chanlabs = ica_chans;
 178 				S.remlabs  = ica_rem;
 179 				S.remfname = ['rem_' S.D];
 180 
 181 				if ~mem_flag
 182 					D = meg_ica_artefact(S);
 183 				elseif mem_flag
 184 					save('tmp.mat','S');
 185 					! matlab -nodesktop -nosplash -r "spm eeg; load tmp.mat; meg_ica_artefact(S); exit"
 186 					try delete('tmp.mat'); end
 187 					D.fname=S.remfname;
 188 				end %mem_flag
 189 
 190 %%%%% Done!
 191 
 192 			end % if icaflag
 193 			
 194 			% The rest of pre-processing filter, epoch, average, etc.
 195 			%  would go here.
 196 	
 197 			
 198 		end % typn
 199 	end % ses
 200 end % sub
 201 		
 202