Attachment 'Session4-script.m'
Download 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
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.