Attachment 'cbu_meeg_spm5_pipeline.m'
Download 1 % Batch script (linear pipeline) for MEG analysis in SPM5
2 % Rik Henson MRC CBU Sep 08
3 % Latest Revision July 09
4 %
5 % Note that would benefit from modularisation (eg for stop/start
6 % running), but wanted all steps in one file to ease exposition
7 %
8 % First start SPM with (at the CBU) "spm 5 eeg" from Linux prompt
9
10 clear
11
12 wd = '/imaging/rh01/Methods/CBUMEGDemo/Group'; % !!!REPLACE WITH YOUR DIRECTORY!!!
13 cd(wd)
14
15 addpath /imaging/local/spm/spm5/cbu_updates/
16 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/spm5_cbu_updates/devel/
17 addpath /imaging/local/meg_misc/
18 % Or update direct from SVN http://imaging.mrc-cbu.cam.ac.uk/svn/meg_misc/devel/
19
20 %%%%%%%%%% Different types of sensor
21 sennam{1} = 'mags'; sennam{2} = 'grds'; sennam{3} = 'eeg';
22 Nsen = length(sennam);
23 eegflag = [0 0 1]; % Binary flag for which sensor is EEG
24
25 % Num channels per sensor type (if differs across subjects, can re-specify below)
26 Nchan(1) = 102; Nchan(2) = 204; Nchan(3) = 70;
27
28 % Any user-thresholds for rejection for each of above (Inf=none)
29 thr(1) = Inf; thr(2) = Inf; thr(3) = 120;
30 EOGthr = 120; % EOG threshold (EOG is duplicated in each file)
31
32 %%%%%%%%%% Experiment-specific parameters
33 expwin = [-100 400]; % Epoch window (in ms)
34 expcodes = [1:8]; % Codes in FIF file to extract
35 offset = [ones(1,length(expcodes))*34]; % t=0 offset (eg visual delay)
36 newcodes = [1 1 1 1 2 2 2 2]; % Any new (recoding) of above
37 expcons = [1 0; 0 1; 1 -1; 1/2 1/2]; % Contrasts (on NEW codes)
38 Ncons = size(expcons,1);
39
40 %%%%%%%%%% Raw FIF files on network
41 % (cell of cells for each session of each sub (only one session here))
42 % Eg if instead two sessions per subject:
43 % Fifdata = { { '/megdata/cbu/lfp/meg08_0014/080117/block1_raw.fif';
44 % '/megdata/cbu/lfp/meg08_0014/080117/block2_raw.fif'} };
45 fifdata = {
46 {'/megdata/cbu/lfp/meg08_0014/080117/block4_raw.fif'};
47 {'/megdata/cbu/lfp/meg08_0015/080118/block4_raw.fif'};
48 {'/megdata/cbu/lfp/meg08_0016/080118/block4_raw.fif'};
49 {'/megdata/cbu/lfp/meg08_0038/080131/block4_raw.fif'};
50 {'/megdata/cbu/lfp/meg08_0064/080212/block4_raw.fif'};
51 {'/megdata/cbu/lfp/meg08_0111/080303/block4_raw.fif'};
52 {'/megdata/cbu/lfp/meg08_0112/080303/block4_raw.fif'};
53 {'/megdata/cbu/lfp/meg08_0113/080303/block4_raw.fif'};
54 {'/megdata/cbu/lfp/meg08_0127/080310/block4_raw.fif'};
55 }
56
57
58 %%%%%%%%%% Raw DICOM files on network (one per subject)
59 mridata = {
60 '/mridata/cbu/CBU080769_CBU080769/20080916_095501/Series_002_CBU_MPRAGE';
61 '/mridata/cbu/CBU071255_CBU071255/20071218_150851/Series_002_CBU_MPRAGE';
62 '/mridata/cbu/CBU071151_CBU071151/20071116_140044/Series_002_CBU_MPRAGE';
63 '/mridata/cbu/CBU071229_CBU071229/20071211_112249/Series_002_CBU_MPRAGE';
64 '/mridata/cbu/CBU080034_RIS13667/20080111_102845/Series_002_CBU_MPRAGE';
65 '/mridata/cbu/CBU080292_CBU080292/20080418_161843/Series_002_CBU_MPRAGE';
66 '/mridata/cbu/CBU080358_CBU080358/20080506_151233/Series_003_CBU_MPRAGE';
67 '/mridata/cbu/CBU080413_CBU080413/20080522_124313/Series_002_CBU_MPRAGE';
68 '/mridata/cbu/CBU080802_CBU080802/20080930_090538/Series_002_CBU_MPRAGE';
69 };
70
71
72 %%%%%%%%%% Bad channel names noticed by Operator or post hoc (eeg/meg,sub,ses)
73 badchans{1,1,1} = [1143];
74 badchans{1,2,1} = [1143 0913 1013 0313 0933 2222 0642 1232 1913];
75 badchans{1,3,1} = [1143];
76 badchans{1,4,1} = [1143 2223 0432];
77 badchans{1,5,1} = [1143];
78 badchans{1,6,1} = [1143];
79 badchans{1,7,1} = [1143 2223];
80 badchans{1,8,1} = [1143 1412];
81 badchans{1,9,1} = [1143 2223];
82
83 %%%%%%%%%% Bad EEG channels noticed by Operator or post hoc (in terms of NAME)
84 badchans{2,1,1} = [48 50 73]; % Poor contact
85 badchans{2,2,1} = [7 13 16 24 46 48 57]; % Drifting
86 badchans{2,3,1} = [48]; % Drifting
87 badchans{2,4,1} = [48]; % Drifting
88 badchans{2,5,1} = [15 40 49]; % Poor contact
89 badchans{2,6,1} = [48]; % Drifting
90 badchans{2,7,1} = [48 72]; % No deflection (odd N170 topo)?
91 badchans{2,8,1} = [65 48]; % odd in final topo of N170?
92 badchans{2,9,1} = [20 31 45]; % 45 poor contact; 20+31 highfreq noise
93
94 % (Just happens that any EEG Channels 61 onwards are numbered 65 onwards in our lab!)
95 for sub=1:size(badchans,2)
96 for ses=1:size(badchans,3)
97 f = find(badchans{2,sub,ses}>60);
98 badchans{2,sub,ses}(f) = badchans{2,sub,ses}(f)-4;
99 end
100 end
101
102 %%%%%%%%%%% subject-specific parameters
103 %dosubs = [1:9]; % Subjects to run, eg for subset: dosubs = [5:8];
104 dosubs = [1:8]; % Exclude sub9 because MRI segmentation requires manual re-
105 % positioning first (can try yourself; see MRI section below!)
106 % (in fact for many, eg sub5 too, normalisation is much
107 % better after manual repositioning)
108 Nsub = length(dosubs);
109
110 % Below assumes that two EOG channels (VEOG and HEOG), which will be
111 % appended at end of every data-file. Note that this script does not yet
112 % handle concurrent ECG (see Jason Taylor for how to do this)
113 for s = 1:Nsub
114 subnum = dosubs(s);
115 for m = 1:Nsen;
116 Nsubchan(s,m) = Nchan(m); % In case different number of channels (eg EEG) for some subjects
117 chan_thr{s,m} = [ones(1,Nsubchan(s,m))*thr(m) EOGthr EOGthr];
118 end
119 end
120
121 VitECapsules = ones(1,max(dosubs)); % Whether MRIs have Vitamin E capsule (all here)
122
123 %%%%%%%%%%% control flags (1=do step; 0=omit step)
124 maxfiltflag = 1; maxmvcompflag = 1; maxtransflag = 1;
125 usemaxmvcompflag = 0; usemaxtransflag = 1;
126 chancheckflag = 0; respmflag = 1; reprocessflag = 1;
127 getmriflag = 1;
128 write2Dflag = [1 0 1]; write3Dflag = 1;
129 invertflag = [1 1 1]; % Invert all modalities
130 fuseinvertflag = 1; % Whether want to simultaneously invert (fuse) modalities (sennam)
131 % (Henson et al, 2009b)
132 grpinvertflag = 1; % Whether invert using group priors (after individual inversions)
133 % (Litvak & Friston, 2008, Neuroimage)
134 fmriinvertflag = 1; % Whether want to use fMRI priors on inversion
135 % (Henson et al, submitted)
136
137 %%%%%%%%%% General parameters
138 % Maxfilter
139 downsample_maxf = '-ds 4'; % maxfilter downsampling factor (eg 4 to save space)
140 % though beware of effect of downsampling on triggers
141 % (http://imaging.mrc-cbu.cam.ac.uk/meg/CbuSpmParameters)
142 %downsample_maxf = ''; % (if none)
143 format_maxf = '-format float'; % maxfilter data format (native = float32)
144 %format_maxf = '-format short'; % (if to save space)
145 trans_offset_maxf = [0 -13 +6]; % New centre for trans default
146 % (see our maxfilter WIKI)
147
148 % SPM
149 downsample_spm = 100; % spm new sample rate (Hz; should be at
150 % least twice filt_cut_spm below)
151 filt_cut_spm = 40; % Lowpass filter cutoff for SPM (Hz)
152 dformat_spm = 'int16'; % data format for spm (eg to save space)
153 sensor_smooth_spm = [5 5 20]; % Smoothing in SensorSpace-Time [mm mm ms];
154 source_smooth_spm = [10 10 10]; % Smoothing in Source space [mm mm mm]
155
156 %%%%%%%%%%% Initialisation
157 subsphere = [];
158 allsssbad = cell(Nsub,1);
159 nbadchan = cell(Nsen,Nsub);
160 nrejects = cell(Nsen,Nsub);
161 nevents = cell(Nsen,Nsub);
162
163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164 %%%%%%%%%%% Preprocess MEG/EEG data
165 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
166
167 for sub = 1:Nsub
168
169 subnum = dosubs(sub);
170 Nses(sub) = size(fifdata{subnum},1);
171
172 cd(wd)
173 swd = sprintf('sub%d',subnum);
174 try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end
175
176 for ses = 1:Nses(sub)
177
178 allsssbad{sub} = cell(1,Nses(sub));
179
180 fnm_stem = sprintf('sub%d_ses%d',subnum,ses);
181
182 if maxfiltflag
183
184 %%%%% Use autobad to detect bad channels from first 20s, and select from output
185 rawfile = fifdata{subnum}{ses};
186 sssfile = sprintf('%s_sss.fif',fnm_stem)
187 logfile = sprintf('%s_autobad.log',fnm_stem)
188 headptsfile = sprintf('%s_headpoints.txt',fnm_stem)
189
190 if exist(sssfile), delete(sssfile), end
191 if exist(logfile), delete(logfile), end
192 [Centre,Radius] = meg_fit_sphere(rawfile,pwd,headptsfile);
193 disp(sprintf('Subject %d (%d): Session %d: Sphere centre [%d %d %d] and radius %d mm',subnum,sub,ses,Centre(1),Centre(2),Centre(3),Radius))
194 subsphere(sub,:) = [Centre Radius];
195
196 disp('Checking first 20s for bad channels...')
197 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d -autobad 20 -skip 21 9999999 -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),logfile))
198
199 badfile = sprintf('%s_badchan.txt',fnm_stem)
200 if exist(badfile), delete(badfile), end
201 eval(sprintf('!cat %s | sed -n ''/Static/p'' | cut -f 5- -d '' '' > %s',logfile,badfile));
202 sssbadchans{sub}{ses} = unique([textread(badfile,'%d')' badchans{1,subnum,ses}]);
203 for bc=1:length(sssbadchans{sub}{ses}); allsssbad{sub}{ses} = [allsssbad{sub}{ses} sprintf(' %04d',sssbadchans{sub}{ses}(bc))]; end
204
205 %%%%% Now mark bad channels and run ST
206 logfile = sprintf('%s_ssst.log',fnm_stem)
207 if exist(sssfile), delete(sssfile), end % Don't need previous file from autobad
208 sssfile = sprintf('%s_ssst.fif',fnm_stem)
209 if exist(logfile), delete(logfile), end
210 posfile = sprintf('%s_mvpos.txt',fnm_stem)
211
212 disp(['Now maxfiltering with bad channels (and tSSS): ' allsssbad{sub}{ses}])
213 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -headpos -hpistep 1000 -hpisubt amp -hp %s %s -v | tee %s',rawfile,sssfile,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},posfile,downsample_maxf,logfile))
214
215 S.logfile = logfile;
216 [max_d(sub),max_r(sub)] = meg_plot_maxfilter_move(S);
217 % Do check the movement parameters - particularly goodness of fit -
218 % because movement compensation below will be invalidated if poor fit
219 % (also, the 'inter' MaxFilter option is used below in case the HPI
220 % temporarily fails (which it does for sub5), but if your subject's
221 % HPI failed permanently from the start, you might not know this otherwise,
222 % because 'inter' would handle it fine!).
223
224 infile = sssfile;
225
226 %%%%% Now mark bad channels and run MVCOMP
227 if maxmvcompflag
228
229 sss2file = sprintf('%s_mvcomp.fif',fnm_stem)
230 logfile = sprintf('%s_mvcomp.log',fnm_stem)
231 if exist(sss2file), delete(sss2file), end
232 if exist(logfile), delete(logfile), end
233
234 disp(['Now compensating movement... (every 1000ms)'])
235 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -movecomp inter -hpistep 200 -hpisubt amp %s -v | tee %s',rawfile,sss2file,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},downsample_maxf,logfile))
236 % eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -frame head -origin %d %d %d %s -autobad off -bad %s -st 4 -movecomp -hpistep 1000 -hpisubt amp %s -v | tee %s',rawfile,sss2file,Centre(1),Centre(2),Centre(3),format_maxf,allsssbad{sub}{ses},downsample_maxf,logfile))
237
238 infile = sss2file;
239 end
240
241 %%%%% Now Trans Default (if requested)
242 if maxtransflag
243 trans_origin = Centre(1:3) + trans_offset_maxf;
244 logfile = sprintf('%s_trans.log',fnm_stem)
245 sss3file = sprintf('%s_trans.fif',fnm_stem)
246 if exist(sss3file), delete(sss3file), end
247 if exist(logfile), delete(logfile), end
248
249 disp(['Now trans-ing to DEFAULT ' num2str(round(trans_origin))])
250 eval(sprintf('!/neuro/bin/util/maxfilter -f %s -o %s -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat -autobad off -trans default -frame head -origin %d %d %d %s -v -force | tee %s',infile,sss3file,trans_origin(1),trans_origin(2),trans_origin(3),format_maxf,logfile))
251 end
252 end
253
254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255 %%%%% Olaf's channel checking (have not used for a while; should put in SVN!)
256 if chancheckflag
257 addpath /imaging/olaf/MEG/artefact_scan_tool/
258 global fiffs tasks values criteria sensors datapara figs;
259 ini_check4artefacts; % Initialisation of variables
260 epochlength = 5;
261 amplitude_threshold_mag = 5000;
262 amplitude_threshold_gra = 10000;
263 out_dir = fullfile(wd,'ChanCheck');
264 addfiff(infile,fnm_stem);
265 addtask('maxminamp'); % maximum minus minimum amplitude within epoch
266 addtask('maxmingrad'); % maximum signal gradient (amplitude difference between successive data points in time) with epoch
267 addtask('threshold'); % number of epochs with max-min amplitudes above specified threshold
268 addtask('crosscorr'); % intercorrelation matrix for magnetometers and gradiometers separately
269 addtask('trendamp'); % linear trend in max-min differences across all epochs
270 addtask('trendgrad'); % linear trend in maximum signal gradients across all epochs
271 check4artefacts; % That's where it all happens... lean back and enjoy!
272 end
273 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
274
275 %%%%% Read into SPM (writes separate files for MEG and EEG)
276
277 clear S
278 if usemaxtransflag
279 S.Fdata = sprintf('%s_trans.fif',fnm_stem);
280 if maxmvcompflag
281 S.HPIfile = sprintf('%s_mvcomp.fif',fnm_stem); % If trans'ed (hpi lost)
282 else
283 S.HPIfile = sprintf('%s_sss.fif',fnm_stem); % If trans'ed (hpi lost)
284 end
285 elseif usemaxmvcompflag
286 S.Fdata = sprintf('%s_mvcomp.fif',fnm_stem);
287 S.HPIfile = [];
288 else
289 S.Fdata = sprintf('%s_sss.fif',fnm_stem);
290 S.HPIfile = [];
291 end
292 % meg_plot_fifpoints(S); % If want to check points
293 S.Fchannels = 'FIF306_setup.mat';
294 S.veogchan = [62 63]; S.veog_unipolar = 0;
295 S.heogchan = [61 64]; S.heog_unipolar = 0;
296 S.trig_chan = 'STI101';
297 S.twin = [0 Inf];
298 S.eeg_ref = 'average';
299 S.trig_type = 'positive_from_zero';
300 % S.Dtype = 'float32'; % make int16 if files too big
301 S.Fchannels_eeg = fullfile(spm('dir'),'EEGtemplates', sprintf('cbu_meg_%deeg_montage_old.mat',Nsubchan(sub,find(eegflag)))); % use default montage (which simply explains how to display in 2D), prior to new CBU caps in mid-2008
302 % S.Fchannels_eeg = fullfile(wd,swd,sprintf('eeg_montage_sub%02d.mat',subnum)); % use subject-specific (care labels may not be standard)
303
304 if respmflag
305 D = spm_eeg_rdata_FIF(S);
306 basefile = D.fname(1:(end-4));
307 else
308 basefile = S.Fdata(1:(end-4));
309 end
310
311 %%%%% Downsample, then split Mags and Grads (file too big otherwise)
312
313 if reprocessflag
314 clear S; S.D = sprintf('%s-eeg.mat',basefile); disp(S.D)
315
316 %%%%% Take opportunity to relabel EEG channels according to approx 10-20 system
317 %%%%% (see http://imaging.mrc-cbu.cam.ac.uk/meg/EEGChannelNames)
318 D = spm_eeg_ldata(S.D);
319 D.channels.ctf = fullfile(spm('dir'),'EEGtemplates', sprintf('cbu_meg_%deeg_montage_named.mat',Nsubchan(sub,find(eegflag))));
320 tmp = load(D.channels.ctf);
321 D.channels.name = tmp.Cnames; % NB: This assumes channel order
322 % has not been changed from 1-70!
323 spm_eeg_inv_save(D);
324
325
326 S.Radc_new = downsample_spm;
327 D = spm_eeg_downsample(S);
328
329 clear S; S.D = basefile;
330 S.Radc_new = downsample_spm;
331 D = spm_eeg_downsample(S);
332 basefile = D.fname(1:(end-4));
333
334 clear S; S.D = D.fname;
335 D = spm_eeg_splitFIF(S);
336 end
337
338 %%%%% If want to view rawdata (with Danny's tool):
339 % meg_viewdata(D.fname);
340
341 %%%%% Preprocess each sensor type
342
343 for sen = 1:Nsen
344
345 if reprocessflag
346 clear S
347 S.D = sprintf('%s-%s.mat',basefile,sennam{sen}); disp(S.D)
348 S.filter.type = 'butterworth';
349 S.filter.band = 'low'; S.filter.PHz = filt_cut_spm;
350 % S.filter.band = 'stop'; S.filter.PHz = [49 51];
351 S.Dtype = dformat_spm;
352 D = spm_eeg_filter(S);
353
354 %%%%%% !!!Insert ICA here if want to remove EOG/ECG artifacts!!!
355 %%%%%% (See Jason Taylor; assume filtered, downsampled data appropriate)
356
357 clear S; S.D = D.fname;
358 S.events.start = expwin(1);
359 S.events.stop = expwin(2);
360 S.events.types = expcodes;
361 S.events.Inewlist = 0;
362 S.events.resynch = offset;
363 D = spm_eeg_epochs(S);
364
365 % re-classify event-codes
366 for e = 1:length(expcodes);
367 fi = find(D.events.code==expcodes(e));
368 if ~isempty(fi)
369 D.events.code(fi) = newcodes(e);
370 end
371 end
372 D.events.types = unique(D.events.code);
373 D.events.Ntypes = length(D.events.types);
374 save(D.fname,'D')
375 clear S; S.D = D.fname;
376
377 else
378 S.D = sprintf('e_fd%s-%s.mat',basefile,sennam{sen}); disp(S.D)
379 D = spm_eeg_ldata(S.D);
380 end
381
382 %%%%% Threshold channels (eg EOG)
383
384 nc = D.Nchannels;
385 S.thresholds.External_list = 0;
386 S.thresholds.Check_Threshold = 1;
387 S.thresholds.threshold = chan_thr{sub,sen};
388 S.BadThr = 0.2;
389 S.artefact.weighted = 0;
390 if eegflag(sen)==0
391 % S.channels.Bad = unique([D.channels.Bad strcmpi(D.channels.name....sssbadchans{sub}{ses}]); % If dont trust MaxFilter channel reconstruction
392 else
393 S.channels.Bad = unique([D.channels.Bad badchans{2,subnum,ses}]);
394 end
395 D = spm_eeg_artefact(S);
396
397 nbadchan{sen,sub} = [nbadchan{sen,sub} length(D.channels.Bad)];
398 nrejects{sen,sub} = [nrejects{sen,sub} length(find(D.events.reject))];
399
400 %%%%% Re-reference EEG if removed channels
401 if eegflag(sen) & ~isempty(D.channels.Bad)
402 D = spm_eeg_ldata(D.fname);
403 meg_reavg_ref(D);
404 save(D.fname,'D')
405 end
406 sesnam{sen,ses} = D.fname;
407
408 %%%%% Uncomment if want to average and contrast each session separately
409 %%%%% (as well as after merging below)
410 %
411 % clear S; S.D = D.fname;
412 % D = spm_eeg_average(S);
413 %
414 % clear S; S.D = D.fname;
415 % S.c = expcons;
416 % S.WeightAve = 1;
417 % D = spm_eeg_weight_epochs(S);
418 % nevents{sen,sub} = D.events.repl;
419
420 end % end of sen loop
421 end % end of ses loop
422
423 %%%%% Merge sessions, average, contrast
424
425 for sen = 1:Nsen
426 if Nses(sub)>1
427 if ~usemaxtransflag
428 disp(['Concatenating sessions... (NB: merged files will not have correct channel Loc/Orient info unless trans-ed to first or to default!!)']);
429 end
430 clear S; S.D = strvcat(sesnam{sen,:});
431 S.recode = cell(1,Nses(sub));
432 D = spm_eeg_merge(S);
433 clear S; S.D = D.fname;
434 else
435 clear S; S.D = sesnam{sen,1};
436 D = spm_eeg_ldata(S.D);
437 end
438 preavgnam{sen,sub} = fullfile(D.path,S.D);
439 Nsubchan(sub,sen) = length(D.channels.eeg) - length(D.channels.Bad);
440 end
441
442
443 %%%% If want to take RMS of grds before averaging (for Sensor-Time image files)
444 %%%% Noise will not cancel as well, but results less biased by number of trials
445 %%%% (for more info, see end of http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm)
446 %
447 % for sen = 1:Nsen
448 % if strcmp(sennam{sen},'grds')
449 % clear S; S.D = preavgnam{sen,sub};
450 % S.bc = 1;
451 % S.do_write = 1;
452 % D = meg_grds2grms(S); % RMS of grads; just for sensor-level analyses
453 % sennam{Nsen+1} = 'grms';
454 % preavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
455 % nevents{Nsen+1,sub} = nevents{sen,sub};
456 % nrejects{Nsen+1,sub} = nrejects{sen,sub};
457 % nbadchan{Nsen+1,sub} = nbadchan{sen,sub};
458 % Nsubchan(sub,Nsen+1) = Nsubchan(sub,sen);
459 % write2Dflag(Nsen+1) = 1;
460 % invertflag(Nsen+1) = 0;
461 % end
462 % end
463
464 for sen = 1:length(sennam)
465
466 clear S; S.D = preavgnam{sen,sub};
467 D = spm_eeg_average(S);
468
469 clear S; S.D = D.fname;
470 S.c = expcons;
471 S.WeightAve = 0;
472 D = spm_eeg_weight_epochs(S);
473
474 nevents{sen,sub} = D.events.repl;
475 avgnam{sen,sub} = fullfile(D.path,D.fname);
476 end
477 end
478
479 cd(wd)
480
481 for sen = 1:length(sennam)
482 disp(sprintf('%s...',sennam{sen}))
483 for sub = 1:Nsub
484 subnum = dosubs(sub);
485 disp(sprintf('\tSub %d (%d)...',subnum,sub))
486 disp(sprintf('\t\tNevents (per condition): %s',mat2str(nevents{sen,sub})))
487 disp(sprintf('\t\tRejects (per session): %s',mat2str(nrejects{sen,sub})))
488 disp(sprintf('\t\tBadchans (per session): %s',mat2str(nbadchan{sen,sub})))
489 end
490 end
491
492 %%%% If want to take RMS after averaging (cf. commented section above)
493 %%%% In theory more sensitive, but biased by trial numbers and less Gaussian
494 %%%% (See end of http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm page)
495
496 for sen = 1:Nsen
497 if strcmp(sennam{sen},'grds')
498 for sub = 1:Nsub
499 cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
500 clear S; S.D = avgnam{sen,sub};
501 S.bc = 1;
502 S.do_write = 1;
503 S.log = 1; % If your RMS are skewed (see above WIKI page)
504 D = meg_grds2grms(S);
505 avgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
506 nevents{Nsen+1,sub} = nevents{sen,sub};
507 end
508 Nsen=Nsen+1; sennam{Nsen} = 'grms';
509 write2Dflag(Nsen) = 1;
510 end
511 end
512
513 %%%% Write out Sensor-Time image files for each subject
514
515 for sen = find(write2Dflag)
516 for sub = 1:Nsub
517 cd(wd); swd = sprintf('sub%d',dosubs(sub)); cd(swd)
518 disp(avgnam{sen,sub})
519 [pth,nam,ext] = fileparts(avgnam{sen,sub});
520 clear S; S.Fname = [nam ext];
521 S.interpolate_bad = 0;
522 S.n = 32;
523 S.pixsize = 3;
524 S.trialtypes = [1:Ncons];
525 P = spm_eeg_convertmat2ana3D(S);
526 if any(sensor_smooth_spm)
527 if length(sensor_smooth_spm)==1; sensor_smooth_spm = ones(1,3)*sensor_smooth_spm; end
528 for con = 1:size(P,1);
529 [pth,nam,ext] = fileparts(P(con,:));
530 Pout = fullfile(pth,['s' nam ext]);
531 spm_smooth(spm_vol(P(con,:)),Pout,sensor_smooth_spm);
532 Pin = strvcat(P(con,:),Pout);
533 spm_imcalc_ui(Pin,Pout,'((i1+eps).*i2)./(i1+eps)',{[],[],'float32',0}); %Reinsert NaNs
534 SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,Pout);
535 end
536 else
537 for con = 1:size(P,1);
538 SensorTimeImgs{sen}{sub}(con,:) = fullfile(wd,swd,P(con,:));
539 end
540 end
541 end
542 end
543
544 %%%% Create Grand Mean across subjects
545 cd(wd)
546 for sen = 1:length(sennam)
547 clear S; S.P = strvcat(avgnam{sen,:});
548 S.Pout = sprintf('G%d-%s.mat',Nsub,sennam{sen})
549 D = spm_eeg_grandmean(S);
550 end
551
552 eval(sprintf('save preproc%d.mat subsphere nrejects nsubchan nbadchan nevents allsssbad dosubs preavgnam avgnam',length(dosubs)))
553
554 disp('!!!NOW USE SPM-Display-M/EEG button to examine results, eg G8-*.mat files !!!')
555
556 %return % uncommment if want script to finish here
557
558
559 % Once displaying one of the modalities, you can hold shift and select both
560 % conditions 1+2, and click on a posterior right channel to see the M170
561 % difference between faces (1) and scrambled faces (2). You can also select
562 % condition 3, which is the differential ERF/ERP between faces and scrambled,
563 % then press "topography", enter -100ms for initial time, then press and hold
564 % the time slider to see topography of difference at every time point (displayed
565 % in Matlab window) - random noise until ~100ms, when strong differences emerge
566
567
568 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
569 %%%%%%%%%%%%%%%%% 2D sensor x 1D Time SPMs
570 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
571
572 stdir = fullfile(wd,'SensorTimeSPMs');
573 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
574
575 clear grpsubs grpcons imgfiles
576 grpsubs{1} = [1:Nsub];
577 % If want to split subjects, eg into two groups, eg:
578 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
579 grpcons{1} = [1:2]; % Include each condition for mags
580 grpcons{3} = [1:2]; % Include each condition for eeg
581 grpcons{4} = [1:2]; % Diff of RMS (be careful if diff number of trials for diff conditions)
582 %grpcons{4} = [3]; % Or RMS of diff (not equivalent to above)
583
584 Ngrp = length(grpsubs);
585 for sen = find(write2Dflag)
586 outdir = fullfile(stdir,sennam{sen})
587 for grp = 1:Ngrp
588 for sub = 1:length(grpsubs{grp});
589 subnum = grpsubs{grp}(sub);
590 imgfiles{grp}{sub} = SensorTimeImgs{sen}{subnum}(grpcons{sen},:);
591 end
592 end
593 meg_batch_anova(imgfiles,outdir);
594 end
595
596 disp('!!!NOW press Results button in SPM!!!')
597
598 %return
599
600 % You need to know a bit about the SPM Results interface, but press the Results
601 % button (make sure SPM is in "EEG" mode) (or better still, use the "ns" toolbox
602 % under "Toolbox" button - see WIKI page below) and select the SPM.mat file for
603 % one of the modalities, select the default effects of interest F-contrast
604 % (because we don't care about polarity, at least for Mags and EEG), choose
605 % 0.001 uncorrected. Then press "Volume" to get table of stats (ignore brain image).
606 % You should see, at least in Mags and EEG, frontal and posterior clusters
607 % (simply of different polarity) maximal around 170ms. Few individual
608 % "voxels" survive whole-image correction because only 8 subjects (so low
609 % df's, ie low power, and RFT conservative) - though the space-time extent
610 % of clusters around 170ms do survive correction at the cluster level
611 % for Mags and EEG (using "ns" toolbox); For GRMS, you can evaluate a
612 % T-contrast [1 -1] for Faces>Scrambled, threshold at p<001 uncorrected
613 % (using ns toolbox) and find some corrected clusters around 170ms,
614 % though the results are not as convincing as Mags or EEG (probably
615 % related to RMS issues: see end of WIKI page below)
616 %
617 % Correction for height would also be significant if you used SVC because
618 % of an a priori timewindow of interest around 170ms (again see WIKI page below).
619 %
620 % For more info on interpreting these SPMs, see
621 % http://imaging.mrc-cbu.cam.ac.uk/meg/SensorSpm
622 %
623 % Also, for an example use with EEG data, see Henson et al (2008), Neuroimage.
624 %
625 % Note that you can do SVC for a priori timewindows of interest, but the main
626 % use of such SPMs is to define timewindows of interest (in the absence of
627 % a priori ones) for the subsequent source localisation below...
628 %
629 % (Note that these SPMs can also handle different numbers and locations of channels
630 % across subjects, provided their locations are digitised in same space, because the
631 % data are interpolated to same grid. For MEG data, particularly gradiometers, the
632 % results benefit from using 'trans -default' above, to realign across different
633 % headpositions for different subjects - see Smith et al (2009), HBM Abstract.
634
635
636 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
637 %%%%%%%%%%%%%%%%% Get MRIS
638 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
639
640 for sub=1:Nsub
641 subnum = dosubs(sub);
642 swd = sprintf('sub%d',subnum);
643 mwd{subnum} = fullfile(wd,swd,'MRI');
644 try cd(mwd{subnum}); catch eval(sprintf('!mkdir %s',mwd{subnum})); cd(mwd{subnum}); end
645
646 if getmriflag
647 hdr = spm_select('FPList',mridata{subnum},'^*.dcm');
648 hdr = spm_dicom_headers(hdr);
649
650 patsex{subnum} = hdr{1}.PatientsSex;
651 patage{subnum} = hdr{1}.PatientsAge; % NB At time of scanning, not time of MEG!
652
653 spm_dicom_convert(hdr,'all','flat','nii');
654
655 mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
656 % (filter setting just to ensure picks up raw structural - assuming CBU onvention)
657 mrifile{sub} = strtrim(mrifile{sub}(1,:)); % Take first if multiple previous
658
659 disp('!!!GOOD IDEA TO MANUALLY REPOSITION STRUCTURAL VIA DISPLAY BUTTON!!!')
660 % pause
661
662 if VitECapsules(subnum)
663 meg_headmask(mrifile{sub},'display'); % Remove any VitE capsules
664 disp('!!!Check difference image!!!');
665 % pause
666 eval(sprintf('!mv %s_masked.nii %s',mrifile{sub}(1:end-4),mrifile{sub}))
667 end
668 patage
669 patsex
670
671 else
672 mrifile{sub} = spm_select('FPList',mwd{subnum},'^s.*-0[0-9].nii$');
673 mrifile{sub} = strtrim(mrifile{sub}(1,:)); % Take first if multiple previous
674 wmrifile{sub} = spm_select('FPList',mwd{subnum},'^wms.*-0[0-9].nii$');
675 wmrifile{sub} = strtrim(wmrifile{sub}(1,:)); % Take first if multiple previous
676 end
677 end
678 spm_check_registration(strvcat(mrifile)); % Only feasible for <~20 images
679 try, spm_check_registration(strvcat(wmrifile)); end % (sub5 bit too warped at back)
680
681 disp('!!!NOW DISPLAY MRI IMAGE in SPM AND DETERMINE FIDUCIALS!!!')
682
683 %%%%%%%%%% Coordinates of fiducials from above MRIs (nasion, left, right)
684 % Below are from using MRI as is (without manually reorienting) - very approximate
685 % For information about manual repositioning, see
686 % http://imaging.mrc-cbu.cam.ac.uk/meg/RepositioningMRIs
687 smri_fids{1} = [ 3 111 22; -78 18 -31; 77 11 -33];
688 smri_fids{2} = [ 6 107 -28; -69 26 -41; 68 14 -45];
689 smri_fids{3} = [-2 124 -23; -78 22 -52; 70 21 -58];
690 smri_fids{4} = [ 3 109 -8; -72 18 -41; 72 18 -47];
691 smri_fids{5} = [-15 103 -33; -79 3 -38; 62 13 -46];
692 smri_fids{6} = [-4 109 -20; -71 15 -32; 63 11 -35];
693 smri_fids{7} = [-1 93 -5; -72 14 -44; 70 10 -44];
694 smri_fids{8} = [ 1 99 -22; -73 8 -27; 69 10 -42];
695 smri_fids{9} = [-3 125 -24; -76 22 -55; 69 22 -58];
696 % Below are after my manual reorienting (just for your reference)
697 %smri_fids{1} = [ 0.7 80.0 -31.2; -76.9 -14.3 -45.9; 79.9 -19.2 -41.0];
698 %smri_fids{2} = [ 1.3 78.7 -30.9; -71.8 -2.5 -58.0; 71.5 1.8 -58.0];
699 %smri_fids{3} = [ 1.7 88.8 -24.6; -74.0 -1.4 -66.1; 71.7 -2.8 -66.1];
700 %smri_fids{4} = [-1.6 78.9 -22.4; -74.8 -13.3 -57.7; 71.6 -7.2 -54.6];
701 %smri_fids{5} = {TBA};
702 %smri_fids{6} = [-2.4 79.1 -27.1; -72.1 -9.8 -44.6; 64.3 -9.8 -51.9];
703 %smri_fids{7} = [-2.3 73.2 -7.0; -71.8 -2.0 -46.9; 70.2 -4.9 -49.8];
704 %smri_fids{8} = [ 0.0 78.1 -30.3; -71.3 -1.5 -48.5; 72.8 6.0 -51.5];
705 %smri_fids{9} = [ 0.0 88.9 -31.1; -76.8 -15.5 -56.6; 78.0 -12.7 -55.0];
706
707
708 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
709 %%%%%%%%%%%%%%%%% Localise on Mesh
710 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
711
712 addpath /imaging/local/Brainstorm/Toolbox
713
714 % !!! There are several options for source localisation, eg choice of
715 % forward model (eg Spheres or BEMs), and choice of inversion method (eg
716 % MSP, COH, MNM). Below illustrates just two possible inversions (each
717 % inversion is indexed by "val"), the second of which has been commented
718 % out because BEMs take a VERY LONG time (and don't work well for EEG yet).
719 % Of course many more options can be tried (and compared using the model-evidence,
720 % recorded in mod_evd below, provided that you don't also change the
721 % data, eg epoch inverted).
722 %
723 % For more information about forward models, see
724 % http://imaging.mrc-cbu.cam.ac.uk/meg/SpmForwardModels
725 % For more information about inversion parameters, see
726 % Friston et al (2008) Neuroimage
727 % (PDF available here http://imaging.mrc-cbu.cam.ac.uk/meg/SpmAnalysis)
728
729 %%%%%%%%%%%%%%%%%%%%
730 val = 1; % Spherical forward model using canonical cortical mesh but individual head meshes
731 comment{val} = 'Concentric Spheres';
732 cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
733 CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
734 iskull{val} = 'Ind'; % Individual (SPM-created) inner skull mesh
735 oskull{val} = ''; % (Outer skull irrelevant for spherical models)
736 scalp{val} = 'Ind'; % Individual (SPM-created) scalp mesh
737 formod{val} = 'Sph'; % Concentric spheres
738 surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all sensor-types
739
740 invcons{val} = [1 2]; % Conditions to invert
741 invtype{val} = 'GS'; % Type of inversion (This is MSP, with a Greedy Search (GS))
742 invtwin{val} = expwin; % Invert whole epoch
743
744 %%%%% Options for subsequent time-freq contrast
745 contwin{val}{1} = [140 190]; % Time window for contrast
746 contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
747 conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
748
749 %%%%%%%%%%%%%%%%%%%%
750 val = 2; % As val=1 except IID (ie standard minimum norm) used to invert
751 comment{val} = 'Concentric Spheres';
752 cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
753 CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
754 iskull{val} = 'Ind'; % Individual (SPM-created) inner skull mesh
755 oskull{val} = ''; % (Outer skull irrelevant for spherical models)
756 scalp{val} = 'Ind'; % Individual (SPM-created) scalp mesh
757 formod{val} = 'Sph'; % Concentric spheres
758 surfit(val,:)= [2 2 2]; % Use 2nd surface to fit sphere (here, innerskull) for all data-types
759
760 invcons{val} = [1 2]; % Conditions to invert
761 invtype{val} = 'IID'; % Type of inversion
762 invtwin{val} = expwin; % Invert whole epoch
763
764 %%%%% Options for subsequent time-freq contrast
765 contwin{val}{1} = [140 190]; % Time window for contrast
766 contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
767 conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
768
769 %%%%%%%%%%%%%%%%%%%%
770 %val = 2; % BEM using canonical meshes (Warning: BEMs take a long time to compute!)
771 %
772 %comment{val} = 'Boundary Element Model';
773 %cortex{val} = 'Can'; % Canonical (inverse-normalised template) cortical mesh
774 %CtxSize(val) = 8196; % One of 5124, 8196 or 20484 dipoles in canonical cortical mesh
775 %iskull{val} = 'Can'; % Canonical (inverse-normalised template) inner skull mesh
776 %oskull{val} = 'Can'; % Canonical (inverse-normalised template) outer skull mesh
777 %scalp{val} = 'Can'; % Canonical (inverse-normalised template) scalp mesh
778 %formod{val} = 'BEM'; % Boundary Element Model
779 %surfit(val,:)= [2 2 4]; % Surfaces for BEM: up to 2nd for MEG (iskull) and 4th for EEG (+oskull+scalp)
780 %invcons{val} = [1 2]; % Conditions to invert
781 %invtype{val} = 'GS'; % Type of inversion (This is MSP, with a Greedy Search (GS))
782 %invtwin{val} = expwin; % Invert whole epoch
783
784 %%%%% Options for subsequent time-freq contrast
785 %contwin{val}{1} = [140 190]; % Time window for contrast
786 %contwin{val}{2} = []; % Frequency window for contrast ([]=all freqs)
787 %conengy{val} = 'evoked'; % or 'induced', but spm_eeg_invert_fuse doesnt handle yet
788
789 % Other options you can consider are "overlapping spheres" for the
790 % forward model (see spm_eeg_inv_BSTfwdsol.m), minimumum norm (IID) for
791 % the inversion, etc. See Henson et al (2009), Neuroimage, for more options.
792
793 Nval = length(comment);
794 mod_evd=[];
795
796 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
797
798 % If want to pass prior estimates of sensor-level error (noise), eg
799 % estimates from empty-room data in case of MEG(see Henson et al, 2009b)...
800 tmp = load('/imaging/rh01/Methods/MEG_Group/MeanCn.mat');
801 Cn{1} = tmp.Cn{1};
802 Cn{2} = tmp.Cn{2};
803 %Cn{1} = tmp.Cn{1}-eye(size(tmp.Cn{1}))*mean(diag(tmp.Cn{1})); % Orthogonalise
804 %Cn{2} = tmp.Cn{2}-eye(size(tmp.Cn{2}))*mean(diag(tmp.Cn{2})); % Orthogonalise
805 % !!!Care: Cn only relevant to SSSt and 40Hz lowpass (see Henson et al, 2009b)!!!
806
807 for sub = 1:Nsub
808
809 cd(wd)
810 subnum = dosubs(sub);
811 swd = sprintf('sub%d',subnum);
812 cd(swd)
813
814 % If want to pass prior estimates of sensor-level error (noise)...
815 % !!!Care - sensor order in sennam must match order in Cn, ie Cn{1}=mags, Cn{2}=grads!!!
816 Ce{1} = {speye(Nsubchan(sub,1)) Cn{1}}; % Mags (Cn assumes no bad channels!!!!)
817 Ce{2} = {speye(Nsubchan(sub,2)) Cn{2}}; % Grds (Cn assumes no bad channels!!!!)
818 Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
819 % Otherwise, do not pass anything, corresponding to just IID ('white') noise only:
820 % Ce = cell(1,Nsen);
821
822 for val = 1:Nval
823
824 for sen = find(invertflag) % Nsen
825
826 clear S;
827 S.D = char(preavgnam{sen,sub});
828 % This is using all trials to estimate sources, ie total (induced+evoked) power.
829 % This gives more sample points in total, and the pure evoked part can always be
830 % extracted when estimating time-freq contrasts afterwards (see below). If you want
831 % to localise purely evoked power, then use instead:
832 % S.D = char(avgnam{sen,sub});
833 % Also may prefer "mvcomp" rather than "trans" data, if don't trust "trans default"
834
835 D = spm_eeg_ldata(S.D); disp(S.D)
836
837 %If want to clear previous localisations (safer?)
838 if val==1
839 if isfield(D,'inv'); D = rmfield(D,'inv'); save(D.fname,'D'); end
840 if isfield(D,'con'); D = rmfield(D,'con'); save(D.fname,'D'); end
841 end
842
843 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initialise
844
845 disp(comment{val})
846 D.val = val; D.fname
847 D.inv{val}.method = 'imaging';
848 load('defaults_eeg_inv.mat')
849 D.inv{val} = invdefts.imag;
850 D.inv{val}.date = mat2str(fix(clock));
851 D.inv{val}.comment = comment{val};
852
853 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MRI
854
855 %%%%%%%%%%%%%%%%%%% Meshing 1: Normalising native MRI
856
857 D.inv{val}.mesh.sMRI = mrifile{sub}; % In case many from previous
858
859 [pth,nam,ext] = spm_fileparts(D.inv{val}.mesh.sMRI);
860 mri_stem = nam;
861 wmrifile{sub} = fullfile(pth,['wm' nam ext]);
862 if ~exist(wmrifile{sub})
863 D = spm_eeg_inv_spatnorm(D);
864 else
865 def_name = [nam '_sn_1.mat'];
866 isndef_name = [nam '_inv_sn_1.mat'];
867 D.inv{val}.mesh.def = fullfile(pth,def_name);
868 D.inv{val}.mesh.invdef = fullfile(pth,isndef_name);
869 D.inv{val}.mesh.nobias = fullfile(pth,['m' nam ext]);
870 D.inv{val}.mesh.wmMRI = fullfile(pth,['wm' nam ext]);
871 end
872
873 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating binary images for scalp (and inner skull)
874
875 D.inv{val}.mesh.msk_scalp = spm_select('List',mwd{subnum},'^*scalp\.nii$');
876
877 if isempty(D.inv{val}.mesh.msk_scalp)
878 D = spm_eeg_inv_getmasks(D);
879 else % assume other masks exist too!
880 D.inv{val}.mesh.msk_scalp = fullfile(mwd{subnum},D.inv{val}.mesh.msk_scalp);
881 D.inv{val}.mesh.msk_iskull = spm_select('FPList',mwd{subnum},'^*iskull\.nii$');
882 D.inv{val}.mesh.msk_cortex = spm_select('FPList',mwd{subnum},'^*cortex\.nii$');
883 D.inv{val}.mesh.msk_flags = struct('img_norm',0,'ne',1,'ng',2,'thr_im',[.5 .05]);
884 end
885 spm_check_registration(char(D.inv{val}.mesh.msk_scalp,D.inv{val}.mesh.msk_iskull,D.inv{val}.mesh.msk_cortex,D.inv{val}.mesh.sMRI));
886
887 %%%%%%%%%%%%%%%%%%% Meshing 2: Creating meshes
888
889 D.inv{val}.mesh.template = 0;
890 D.inv{val}.mesh.canonical = 1;
891 D.inv{val}.mesh.Msize = CtxSize(val);
892
893 % Prompt for canonical or user-specified cortical meshs (SPM does not currently
894 % produce individual ones, because automatic meshing of cortex is difficult!)
895
896 if strcmp(cortex{val},'Can')
897 Tmesh = load(fullfile(spm('dir'),'EEGtemplates',sprintf('mni_mesh_cortex_%d.mat',CtxSize(val))));
898 D.inv{val}.mesh.tess_mni.vert = Tmesh.vert;
899 D.inv{val}.mesh.tess_mni.face = uint16(Tmesh.face);
900 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
901 D.inv{val}.mesh.tess_ctx.vert = Tmesh.vert;
902 D.inv{val}.mesh.tess_ctx.face = uint16(Tmesh.face);
903 elseif strcmp(cortex{val},'Load')
904 Tmesh = spm_select(1,'mat','Select cortex mat file containing vert and face',[],pwd,'.*mat');
905 D.inv{val}.mesh.tess_ctx.vert = Tmesh.vert;
906 D.inv{val}.mesh.tess_ctx.face = uint16(Tmesh.face);
907 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.invdef);
908 D.inv{val}.mesh.tess_mni.vert = Tmesh.vert;
909 D.inv{val}.mesh.tess_mni.face = uint16(Tmesh.face);
910 else
911 error('No cortical mesh specified????')
912 end
913
914 iskull_mesh_name = fullfile('MRI',[mri_stem '_iskull_mesh.mat']);
915 scalp_mesh_name = fullfile('MRI',[mri_stem '_scalp_mesh.mat']);
916
917 % Prompt for canonical, individual (from SPM) or user-specified inner skull
918 % (Use canonical if worried about canonical cortex piercing innner skull)
919
920 if strcmp(iskull{val},'Can')
921 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_iskull_2562.mat'));
922 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
923 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
924 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
925 elseif strcmp(iskull{val},'Ind')
926 if ~exist(iskull_mesh_name)
927 D = spm_eeg_inv_getmeshes(D,1);
928 vert = D.inv{val}.mesh.tess_iskull.vert;
929 face = D.inv{val}.mesh.tess_iskull.face;
930 save(iskull_mesh_name,'vert','face');
931 else
932 Tmesh = load(iskull_mesh_name);
933 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
934 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
935 end
936 elseif strcmp(iskull{val},'Load')
937 Tmesh = spm_select(1,'mat','Select iskull mat file containing vert and face',[],pwd,'.*mat');
938 D.inv{val}.mesh.tess_iskull.vert = Tmesh.vert;
939 D.inv{val}.mesh.tess_iskull.face = uint16(Tmesh.face);
940 else
941 warning('No inner skull mesh specified')
942 end
943
944 % Prompt for canonical or user-specified outer skull (SPM does not currently
945 % produce individual ones, because segmenting outer skull from scalp difficult)
946 % (Note that outer skull only really necessary for EEG BEMs)
947
948 if strcmp(oskull{val},'Can')
949 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_oskull_2562.mat'));
950 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
951 D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
952 D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
953 elseif strcmp(oskull{val},'Load')
954 Tmesh = spm_select(1,'mat','Selectoskull mat file containing vert and face',[],pwd,'.*mat');
955 D.inv{val}.mesh.tess_oskull.vert = Tmesh.vert;
956 D.inv{val}.mesh.tess_oskull.face = uint16(Tmesh.face);
957 else
958 warning('No outer skull mesh specified')
959 end
960
961 % Prompt for canonical, individual (from SPM) or user-specified scalp
962 % (Use canonical if worried about canonical cortex or inner skull piercing scalp)
963
964 if strcmp(scalp{val},'Can')
965 Tmesh = load(fullfile(spm('dir'),'EEGtemplates','mni_mesh_scalp_2562.mat'));
966 Tmesh.vert = spm_get_orig_coord(Tmesh.vert,D.inv{val}.mesh.def);
967 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
968 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
969 elseif strcmp(scalp{val},'Ind')
970 if ~exist(scalp_mesh_name)
971 D = spm_eeg_inv_getmeshes(D,2);
972 vert = D.inv{val}.mesh.tess_scalp.vert;
973 face = D.inv{val}.mesh.tess_scalp.face;
974 save(scalp_mesh_name,'vert','face');
975 else
976 Tmesh = load(scalp_mesh_name);
977 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
978 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
979 end
980 elseif strcmp(scalp{val},'Load')
981 Tmesh = spm_select(1,'mat','Select scalp mat file containing vert and face',[],pwd,'.*mat');
982 D.inv{val}.mesh.tess_scalp.vert = Tmesh.vert;
983 D.inv{val}.mesh.tess_scalp.face = uint16(Tmesh.face);
984 else
985 warning('No scalp mesh specified')
986 end
987
988 % Always individual scalp for datareg
989
990 if ~exist(scalp_mesh_name)
991 S = spm_eeg_inv_getmeshes(D,2); % lastarg==1 returns iskull (not scalp)
992 vert = S.inv{val}.mesh.tess_scalp.vert;
993 face = S.inv{val}.mesh.tess_scalp.face;
994 save(scalp_mesh_name,'vert','face');
995 D.inv{val}.mesh.tess_scalp_datareg.vert = vert;
996 D.inv{val}.mesh.tess_scalp_datareg.face = uint16(face);
997 else
998 Tmesh = load(scalp_mesh_name);
999 D.inv{val}.mesh.tess_scalp_datareg.vert = Tmesh.vert;
1000 D.inv{val}.mesh.tess_scalp_datareg.face = uint16(Tmesh.face);
1001 end
1002
1003 spm_eeg_inv_checkmeshes(D);
1004 save(D.fname,'D');
1005
1006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DataReg
1007
1008 if isempty(D.inv{val}.datareg.fid_coreg)
1009
1010 D.inv{val}.datareg.sensors = D.channels.Loc';
1011 D.inv{val}.datareg.fid_eeg = D.channels.fid_eeg;
1012 D.inv{val}.datareg.headshape = D.channels.headshape;
1013 if ~eegflag(sen)
1014 D.inv{val}.datareg.megorient = D.channels.Orient';
1015 end
1016 D.inv{val}.datareg.scalpvert = D.inv{val}.mesh.tess_scalp_datareg.vert;
1017 D.inv{val}.datareg.fid_mri = smri_fids{subnum};
1018 D = spm_eeg_inv_datareg(D);
1019 fid_err(sub) = sqrt(mean(mean((D.inv{val}.datareg.fid_coreg - D.inv{val}.datareg.fid_mri).^2)))
1020 end
1021
1022 spm_eeg_inv_checkdatareg(D);
1023 save(D.fname,'D');
1024 clear S; S.D=D.fname;
1025 % meg_pretty_datareg(S); % Prettier rendering of head and sensors
1026
1027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Forward Model
1028
1029 rotate3d off
1030 D.inv{val}.method = 'Imaging';
1031 if ~eegflag(sen)
1032 if strcmp(formod{val},'Sph'), meth = 'meg_sphere';
1033 elseif strcmp(formod{val},'BEM'), meth = 'meg_bem';
1034 elseif strcmp(formod{val},'OS'), meth = 'meg_os';
1035 else, error('Unknown forward model!'); end
1036 else
1037 if strcmp(formod{val},'Sph'), meth = 'eeg_3sphereBerg';
1038 elseif strcmp(formod{val},'BEM'), meth = 'eeg_bem';
1039 elseif strcmp(formod{val},'OS'), meth = 'eeg_os';
1040 else, error('Unknown forward model!'); end
1041 end
1042 D.inv{val}.forward.method = meth;
1043 D.inv{val}.forward.surface2fit = surfit(val,sen);
1044 % D.inv{val}.forward.NVertMax = 1000; % If want to save some time!
1045 gain_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatrix_' upper(sennam{sen}) ...
1046 '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
1047 gxyz_name = fullfile(mwd{subnum},[mri_stem '_SPMgainmatxyz_' upper(sennam{sen}) ...
1048 '_' formod{val} '_' num2str(surfit(val,sen)) '.mat'])
1049
1050 % delete(gain_name); % if want to be safe (uncomment if want to save time)
1051 if ~exist(gain_name) %| ~exist(gxyz_name)
1052 [D,OPTIONS] = spm_eeg_inv_BSTfwdsol(D);
1053 eval(sprintf('!mv %s %s',D.inv{val}.forward.gainmat,gain_name))
1054 eval(sprintf('!mv %s %s',D.inv{val}.forward.gainxyz,gxyz_name))
1055 end
1056
1057 D.inv{val}.forward.gainmat = gain_name;
1058 D.inv{val}.forward.gainxyz = gxyz_name;
1059 save(D.fname,'D');
1060
1061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Invert each modality separately
1062
1063 D.inv{val}.inverse = [];
1064 D.inv{val}.inverse.trials = invcons{val};
1065 D.inv{val}.inverse.type = invtype{val};
1066 D.inv{val}.inverse.woi = invtwin{val};
1067 D.inv{val}.inverse.lpf = 0; % Confusing, but hpf refers to cutoff
1068 D.inv{val}.inverse.hpf = 40; % of lowpass (and lpf to cutoff of highpass)
1069 % Some extra parameters that you could play with...(but defaults seem ok!)
1070 % D.inv{val}.inverse.smooth = 0.8;
1071 % D.inv{val}.inverse.Np = 256;
1072 % D.inv{val}.inverse.Han = 0;
1073 % D.inv{val}.inverse.Nm = [];
1074 % D.inv{val}.inverse.NmTOL = -Inf;
1075 % D.inv{val}.inverse.NrTOL = -1;
1076 D.inv{val}.inverse.pQe = Ce{sen};
1077 D.inv{val}.inverse.pQp = [];
1078
1079 D = spm_eeg_invert(D);
1080 % D = spm_eeg_invert_fuse(D); % If want to compare exactly with fused below
1081 save(D.fname,'D');
1082 % mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end); % if want to include model complexity
1083 mod_evd(sub,val,sen) = D.inv{val}.inverse.F(end);
1084
1085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1086
1087 D.inv{val}.contrast.woi = contwin{val}{1};
1088 D.inv{val}.contrast.fboi = contwin{val}{2};
1089 D.inv{val}.contrast.type = conengy{val};
1090 D.val = val;
1091 D = spm_eeg_inv_results(D);
1092 spm_eeg_inv_results_display(D);
1093
1094 if write3Dflag
1095 D.inv{val}.contrast.smooth = source_smooth_spm;
1096 D.inv{val}.contrast.rms = 1;
1097 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1098 D = spm_eeg_inv_Mesh2Voxels(D);
1099 SourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1100 end
1101
1102 % If want to save on filesize
1103 % if val>2
1104 % canvert = D.inv{val-1}.mesh.tess_mni.vert;
1105 % canface = D.inv{val-1}.mesh.tess_mni.face;
1106 % D.inv{val-1}.mesh = [];
1107 % D.inv{val-1}.mesh.tess_mni.vert = canvert;
1108 % D.inv{val-1}.mesh.tess_mni.face = canface;
1109 % D.inv{val-1}.datareg = [];
1110 % clear canvert canface
1111 % end
1112
1113 disp(sprintf('Sub=%d, %s, val=%d: %s',sub,sennam{sen},val,D.inv{val}.comment))
1114 save(D.fname,'D');
1115
1116 invertedflag(sen) = 1;
1117
1118 end % end sen
1119
1120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1121 % Fusion
1122
1123 if fuseinvertflag
1124 DD = {};
1125 for sen = find(invertflag)
1126 DD{sen} = spm_eeg_ldata(preavgnam{sen,sub});
1127 DD{sen}.val = val;
1128 DD{sen}.inv{val}.inverse.pQe = Ce{sen};
1129 DD{sen}.inv{val}.inverse.pQp = [];
1130 end
1131 DD{1}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1132 DD{1}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1133 % DD{1}.inv{val}.inverse.NrTOL = -1;
1134
1135 DD{1}.inv{val}.inverse.OutExt = '_fused';
1136 D = spm_eeg_invert_fuse(DD);
1137 % mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
1138 mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F(end);
1139
1140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1141
1142 D.inv{val}.contrast.woi = contwin{val}{1};
1143 D.inv{val}.contrast.fboi = contwin{val}{2};
1144 D.inv{val}.contrast.type = conengy{val};
1145 D.val = val;
1146 D = spm_eeg_inv_results(D);
1147 spm_eeg_inv_results_display(D);
1148
1149 if write3Dflag
1150 D.inv{val}.contrast.smooth = source_smooth_spm;
1151 D.inv{val}.contrast.rms = 1;
1152 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1153 D = spm_eeg_inv_Mesh2Voxels(D);
1154 SourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1155 end
1156
1157 sennam{Nsen+1} = 'fused';
1158 preavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1159 invertedflag(Nsen+1) = 1;
1160 end
1161
1162 end % end val
1163
1164 end % end subs
1165
1166 % Review model-evidences (note cannot compare fused to individual modality
1167 % inversions, since data differs)
1168 for sub = 1:Nsub
1169 subnum = dosubs(sub);
1170 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1171 disp(sprintf('\t\tFiducial error (RMS mm): %3.2f',fid_err(sub)))
1172 for val=1:Nval
1173 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(mod_evd(sub,val,:))))))
1174 end
1175 end
1176
1177 eval(sprintf('save preproc%d.mat subsphere nrejects Nsubchan nbadchan nevents allsssbad dosubs preavgnam sennam invertflag SourceImgs Nval mod_evd fid_err Ce',length(dosubs)))
1178
1179 % Fid errors may be quite large (~8mm) for some subjects, simply because
1180 % I did not bother to mark them carefully on the MRI. Note that the
1181 % actual coregistration is determined primarily by the digitised headshape;
1182 % the fiducials are only really used to check the quality of registration.
1183
1184 %spm_check_registration(strvcat(wmrifile)) % double-check normaliation worked
1185
1186 disp('!!!Review single-subject localisations if you wish!!!')
1187
1188 % At this point, you can review a single-subject's localisation by
1189 % pressing the "3D source reconstruction" button, and then loading one
1190 % of the ae* files (for one modality, or the "*_fused.mat" for the
1191 % fusion of the modalities). You can then press the "mip" button to see
1192 % the overall reconstruction, or the "display" button under "window" to
1193 % see the power in the time-freq contrast (for the selected
1194 % condition). You can also display a movie and change the start-stop
1195 % times, etc. Don't press the other buttons for inversion, fusion, group
1196 % inversion etc because they will be covered below...
1197
1198 %%%%%%%%%%%%%%%%% Or a quick look at mean over subjects...
1199
1200 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1201 S.fig = 2;
1202 for val=1:Nval
1203 S.val = val;
1204 for sen=find(invertedflag)
1205 S.fig = S.fig+1;
1206 S.P = strvcat(preavgnam{sen,:});
1207 meg_inv_display_con(S);
1208 end
1209 end
1210 lastfig = S.fig;
1211
1212 % You should see bilateral fusiform, maximal anteriorly in MEG but more
1213 % posteriorly in EEG (more occipital), plus a mixture after fusion
1214 % Note that this is just the mean across subjects, and stats (below) differ
1215
1216 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1217 %%%%%%%%%%%%%%%%% 3D SPMs
1218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1219
1220 stdir = fullfile(wd,'SourceSPMs');
1221 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1222
1223 clear grpsubs imgfiles grpcons
1224 grpsubs{1} = [1:Nsub];
1225 % If want to split subjects, eg into two groups, eg:
1226 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1227
1228 grpcons = [1:2];
1229 Ngrp = length(grpsubs);
1230 for val=1:Nval
1231 valdir = fullfile(stdir,sprintf('Inv%d',val));
1232 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1233 for sen = find(invertedflag)
1234 outdir = fullfile(valdir,sennam{sen})
1235 for grp = 1:Ngrp
1236 for sub = 1:length(grpsubs{grp});
1237 subnum = grpsubs{grp}(sub);
1238 imgfiles{grp}{sub} = SourceImgs{sen}{val}{subnum}(grpcons,:);
1239 end
1240 end
1241 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1242 end
1243 end
1244 cd(stdir)
1245 disp('!!!Review the group 3D stats if you wish!!!')
1246
1247 %return
1248
1249 % Again you need to know a bit about the SPM Results interface, but press
1250 % Results and select the SPM.mat file for the "mags" directory in the
1251 % "Inv1" (MSP(GS)) directory, and enter a new T-contrast [1 -1] to find
1252 % voxels in which greater source amplitude for faces than scrambled faces
1253 % Choose 0.001 uncorrected and then press "Volume" to get table of stats.
1254 % You should see a swathe of ventral temporal activity, particularly on right,
1255 % with two clusters that survive correction for multiple comparisons.
1256 % For GRDS, you get smaller, more anterior fusiform clusters; for EEG
1257 % you get much more posterior occipital clusters; for fused results, you
1258 % get a more focal lateral midfusiform. Note that the differences across
1259 % these SPMs are likely more apparent than real, given thresholding and
1260 % the relatively small number of subjects (see Henson et al, 2007,
1261 % Neuroimage and Taylor & Henson, in prep, for more information). MSP
1262 % should also benefit more from group inversion (below).
1263 %
1264 % For the IID inversions (in the "Inv2" directory), the results are more
1265 % consistent across modalities, but also more superficial (eg lateral);
1266 % a well known bias of standard minimum norm.
1267
1268 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1269 %%%%%%%%%%%%%%%%% 3D PPMs
1270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271
1272 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1273 % which do not have any mutliple comparison issues (eg no need for RFT)
1274 % However, it will take a long time...
1275
1276 %for val=1:Nval
1277 % valdir = fullfile(stdir,sprintf('Inv%d',val));
1278 % cd(valdir);
1279 % for sen = find(invertedflag)
1280 % outdir = fullfile(valdir,sennam{sen})
1281 % cd(outdir)
1282 % load('SPM.mat');
1283 % spm_spm_Bayes(SPM);
1284 % end
1285 %end
1286 %disp('!!!Review the group 3D PPMs if you wish!!!')
1287
1288 % Again you need to know a bit about the PPM Results interface, but press
1289 % Results and select an SPM.mat file for one of the modalities and enter
1290 % a new T-contrast [1 -1], but this time, select "Bayesian" rather than
1291 % "classical" when prompted. Then change the default effect size threshold
1292 % to 0, and probability threshold to 0.99. This will show voxels with a
1293 % >99% probability of the increase for faces relative to scrambled being
1294 % greater than zero.
1295
1296
1297 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1298 % Group inversion (Litvak & Friston, 2008)
1299 % This pools MSP priors across subjects, in theory to make localisations
1300 % more consistent (it has no effect on IID inversions).
1301
1302 grp_mod_evd=[];
1303 if grpinvertflag
1304 for sen = find(invertedflag)
1305 clear DD
1306 for sub=1:Nsub
1307 DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
1308 end
1309 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1310 for sub=1:Nsub
1311 Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
1312 DD{sub}.val = val;
1313 DD{sub}.inv{val}.comment = [DD{sub}.inv{val}.comment 'Group'];
1314 DD{sub}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1315 DD{sub}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1316 % DD{sub}.inv{val}.inverse.NrTOL = -1;
1317 DD{sub}.inv{val}.inverse.pQe = Ce{sen};
1318 DD{sub}.inv{val}.inverse.pQp = [];
1319 end
1320 DD = spm_eeg_invert(DD);
1321
1322 % Save outputs
1323 for sub=1:Nsub
1324 D = DD{sub};
1325 D.fname = [D.fname(1:(end-4)) '_grp.mat'];
1326 grppreavgnam{sen,sub} = fullfile(D.path, D.fname);
1327 save(fullfile(D.path, D.fname), 'D');
1328 grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F(end);
1329
1330 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1331
1332 D.inv{val}.contrast.woi = contwin{val}{1};
1333 D.inv{val}.contrast.fboi = contwin{val}{2};
1334 D.inv{val}.contrast.type = conengy{val};
1335 D.val = val;
1336 D = spm_eeg_inv_results(D);
1337 spm_eeg_inv_results_display(D);
1338
1339 if write3Dflag
1340 D.inv{val}.contrast.smooth = source_smooth_spm;
1341 D.inv{val}.contrast.rms = 1;
1342 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1343 D = spm_eeg_inv_Mesh2Voxels(D);
1344 GrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1345 end
1346 save(fullfile(D.path, D.fname), 'D');
1347 end % end sub
1348 end % end val
1349 end % end sen
1350
1351 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1352 % Fuse any Group inversion
1353 if fuseinvertflag
1354 for sub=1:Nsub
1355 Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
1356 clear DD;
1357 for sen = find(invertflag)
1358 DD{sen} = spm_eeg_ldata(grppreavgnam{sen,sub});
1359 end
1360
1361 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1362 for sen = find(invertflag)
1363 DD{sen}.val = val;
1364 DD{sen}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1365 DD{sen}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1366 % DD{sen}.inv{val}.inverse.NrTOL = -1;
1367 DD{sen}.inv{val}.inverse.pQe = Ce{sen};
1368 DD{sen}.inv{val}.inverse.OutExt = '_fused';
1369 DD{sen}.inv{val}.inverse.pQp = {DD{sen}.inv{val}.inverse.QP}; % Use group source priors...
1370 end
1371 DD{1}.inv{val}.inverse.grpopt = 0; % ...so no need to (re)optimise source priors
1372
1373 D = spm_eeg_invert_fuse(DD);
1374 grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F(end);
1375
1376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1377
1378 D.inv{val}.contrast.woi = contwin{val}{1};
1379 D.inv{val}.contrast.fboi = contwin{val}{2};
1380 D.inv{val}.contrast.type = conengy{val};
1381 D.val = val;
1382 D = spm_eeg_inv_results(D);
1383 spm_eeg_inv_results_display(D);
1384
1385 if write3Dflag
1386 D.inv{val}.contrast.smooth = source_smooth_spm;
1387 D.inv{val}.contrast.rms = 1;
1388 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1389 D = spm_eeg_inv_Mesh2Voxels(D);
1390 GrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1391 end
1392 save(fullfile(D.path, D.fname), 'D');
1393 end % end val
1394
1395 grppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1396 end % end sub
1397 sennam{Nsen+1} = 'fused';
1398 end
1399 end
1400
1401 % Review model-evidences (note cannot compare fused to individual modality
1402 % inversions, since data differs)
1403 for sub = 1:Nsub
1404 subnum = dosubs(sub);
1405 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1406 for val=1:Nval
1407 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(grp_mod_evd(sub,val,:))))))
1408 end
1409 end
1410
1411 eval(sprintf('save preproc%d.mat subsphere nrejects Nsubchan nbadchan nevents allsssbad dosubs preavgnam sennam invertflag SourceImgs Nval mod_evd fid_err Ce grp_mod_evd GrpSourceImgs grppreavgnam',length(dosubs)))
1412
1413 %%%%%%%%%%%%%%%%% Quick look at mean over subjects...
1414
1415 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1416 S.fig = lastfig;
1417 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1418 S.val = val;
1419 for sen = find(invertedflag)
1420 S.fig = S.fig+1;
1421 S.P = strvcat(grppreavgnam{sen,:});
1422 meg_inv_display_con(S);
1423 drawnow
1424 end
1425 end
1426 lastfig = S.fig;
1427
1428 % The group MSP results look similar to before, but with higher maximal values.
1429 % The IID results are (as expected) unchanged by group inversion, apart
1430 % from when the group inversions are fused (when there is a small change,
1431 % as expected owing to the third inversion step).
1432
1433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1434 %%%%%%%%%%%%%%%%% 3D SPMs
1435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1436
1437 stdir = fullfile(wd,'GrpSourceSPMs');
1438 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1439
1440 clear grpsubs imgfiles grpcons
1441 grpsubs{1} = [1:Nsub];
1442 % If want to split subjects, eg into two groups, eg:
1443 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1444 grpcons = [1:2];
1445 Ngrp = length(grpsubs);
1446 for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1447 valdir = fullfile(stdir,sprintf('Inv%d',val));
1448 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1449 for sen = find(invertedflag)
1450 outdir = fullfile(valdir,sennam{sen})
1451 for grp = 1:Ngrp
1452 for sub = 1:length(grpsubs{grp});
1453 subnum = grpsubs{grp}(sub);
1454 imgfiles{grp}{sub} = GrpSourceImgs{sen}{val}{subnum}(grpcons,:);
1455 end
1456 end
1457 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1458 end
1459 end
1460 cd(stdir)
1461 disp('!!!Review the group 3D stats if you wish!!!')
1462
1463 %return
1464
1465 % Again you need to know a bit about the SPM Results interface, but press
1466 % Results and select an SPM.mat file for one of the modalities under MSP
1467 % (GS) and enter a new T-contrast [1 -1] to find voxels in which greater source
1468 % amplitude for faces than scrambled faces at 0.001 uncorrected. Results
1469 % are similar, though a bit more reliable, after group inversion than
1470 % previously (though additional orbitofrontal in MAGS and GRDS unclear!).
1471 % The IID results will not have changed (except for the fused results).
1472
1473 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1474 %%%%%%%%%%%%%%%%% 3D PPMs
1475 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1476
1477 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1478 % which do not have any mutliple comparison issues (eg no need for RFT)
1479 % However, it will take a long time...
1480
1481 %for val=1:Nval % NB: Group inversion makes no difference for IID (val=2 here)
1482 % valdir = fullfile(stdir,sprintf('Inv%d',val));
1483 % cd(valdir);
1484 % for sen = find(invertedflag)
1485 % outdir = fullfile(valdir,sennam{sen})
1486 % cd(outdir)
1487 % load('SPM.mat');
1488 % spm_spm_Bayes(SPM);
1489 % end
1490 %end
1491 %disp('!!!Review the group 3D PPMs if you wish!!!')
1492
1493 % Again you need to know a bit about the PPM Results interface, but press
1494 % Results and select an SPM.mat file for one of the modalities and enter
1495 % a new T-contrast [1 -1], but this time, select "Bayesian" rather than
1496 % "classical" when prompted. Then change the default effect size threshold
1497 % to 0, and probability threshold to 0.99. This will show voxels with a
1498 % >99% probability of the increase for faces relative to scrambled being
1499 % greater than zero.
1500
1501 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1502 % adding fMRI priors (Henson et al, submitted)
1503
1504 fmri_grp_mod_evd=[];
1505 if fmriinvertflag
1506
1507 % Create fMRI priors from an SPM image....
1508 % Here we assume a group prior in MNI space, but priors could be
1509 % defined separately for each subject based on their own fMRI data
1510
1511 clear S
1512 S.D = grppreavgnam{1,1}; % Take first subject/modality; doesnt really matter
1513 S.fmri = '/imaging/rh01/Methods/fMRIPriors/fMRI_ANOVA_flip0/UFvsS_05_cor_15.img';
1514 S.gm = '/imaging/local/spm/spm5/canonical/c1single_subj_T1.nii';
1515 S.space = 1; % SPM and GM images are in MNI space
1516 S.hthr = 0.5; S.ethr = 0;
1517 S.ncomp = Inf; S.smooth = 0.2;
1518 S.pflag = 1;
1519
1520 D0 = spm_eeg_inv_fmripriors(S);
1521 fMRI = D0.inv{D0.val}.fmri.priors; % Priors (cell of vectors)
1522
1523 % Note that here, we simply add the fMRI priors to the previous
1524 % group-optimised source priors. One could also invert each subject
1525 % afresh, re-optimising the MSP priors at the same time as the fMRI
1526 % priors (either individually, or as a group), by using the preavgnam
1527 % files (not grppreavgnam ones), and not setting grpopt=0.
1528
1529 for sen = find(invertflag)
1530 clear DD
1531 for sub=1:Nsub
1532 % DD{sub} = spm_eeg_ldata(preavgnam{sen,sub});
1533 DD{sub} = spm_eeg_ldata(grppreavgnam{sen,sub});
1534 end
1535 for val=1:Nval
1536 for sub=1:Nsub
1537 Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
1538 DD{sub}.inv{val}.comment = [DD{sub}.inv{val}.comment 'Group+fMRI priors'];
1539 DD{sub}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1540 DD{sub}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1541 % DD{sub}.inv{val}.inverse.NrTOL = -1;
1542 DD{sub}.val = val;
1543 DD{sub}.inv{val}.inverse.pQe = Ce{sen};
1544 DD{sub}.inv{val}.inverse.pQp = [];
1545 end
1546 DD{1}.inv{val}.inverse.pQp = fMRI;
1547 DD{1}.inv{val}.inverse.pQp{end+1} = DD{1}.inv{val}.inverse.QP;
1548 % DD{1}.inv{val}.inverse.grpopt = 1; % Do (re)optimise source priors
1549 DD{1}.inv{val}.inverse.grpopt = 0; % Do not (re)optimise source priors
1550
1551 DD = spm_eeg_invert(DD);
1552
1553 % Save outputs
1554 for sub=1:Nsub
1555 D = DD{sub};
1556 % D.fname = [D.fname(1:(end-4)) '_grp_fmri.mat'];
1557 D.fname = [D.fname(1:(end-4)) '_fmri.mat'];
1558 fmrigrppreavgnam{sen,sub} = fullfile(D.path, D.fname);
1559 save(fullfile(D.path, D.fname), 'D');
1560 fmri_grp_mod_evd(sub,val,sen) = D.inv{val}.inverse.F1(end);
1561
1562 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1563
1564 D.inv{val}.contrast.woi = contwin{val}{1};
1565 D.inv{val}.contrast.fboi = contwin{val}{2};
1566 D.inv{val}.contrast.type = conengy{val};
1567 D.val = val;
1568 D = spm_eeg_inv_results(D);
1569 spm_eeg_inv_results_display(D);
1570
1571 if write3Dflag
1572 D.inv{val}.contrast.smooth = source_smooth_spm;
1573 D.inv{val}.contrast.rms = 1;
1574 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1575 D = spm_eeg_inv_Mesh2Voxels(D);
1576 FmriGrpSourceImgs{sen}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1577 end
1578 save(fullfile(D.path, D.fname), 'D');
1579 end % end val
1580 end % end sub
1581 end % end sen
1582
1583 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1584 % Fuse using Group and fMRI priors
1585 if fuseinvertflag
1586 for sub=1:Nsub
1587 Ce{3} = {speye(Nsubchan(sub,3))}; % EEG (no noise estimates yet)
1588 clear DD;
1589 for sen = find(invertflag)
1590 DD{sen} = spm_eeg_ldata(fmrigrppreavgnam{sen,sub});
1591 end
1592
1593 for val=1:Nval
1594 for sen = find(invertflag)
1595 DD{sen}.val = val;
1596 DD{sen}.inv{val}.inverse.Nm = []; % Reset any parameters from last inversion
1597 DD{sen}.inv{val}.inverse.Nr = []; % Reset any parameters from last inversion
1598 % DD{sen}.inv{val}.inverse.NrTOL = -1;
1599 DD{sen}.inv{val}.inverse.OutExt = '_fused';
1600 DD{sen}.inv{val}.inverse.pQe = Ce{sen};
1601 DD{sen}.inv{val}.inverse.pQp = {DD{sen}.inv{val}.inverse.QP}; % Use group source priors...
1602 end
1603 DD{1}.inv{val}.inverse.grpopt = 0; % ...so no need to (re)optimise source priors
1604
1605 D = spm_eeg_invert_fuse(DD);
1606 fmri_grp_mod_evd(sub,val,Nsen+1) = D.inv{val}.inverse.F1(end);
1607
1608 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contrast
1609
1610 D.inv{val}.contrast.woi = contwin{val}{1};
1611 D.inv{val}.contrast.fboi = contwin{val}{2};
1612 D.inv{val}.contrast.type = conengy{val};
1613 D.val = val;
1614 D = spm_eeg_inv_results(D);
1615 spm_eeg_inv_results_display(D);
1616
1617 if write3Dflag
1618 D.inv{val}.contrast.smooth = source_smooth_spm;
1619 D.inv{val}.contrast.rms = 1;
1620 % D.inv{val}.contrast.scalefactor = 1; % if want to compare across dif size windows
1621 D = spm_eeg_inv_Mesh2Voxels(D);
1622 FmriGrpSourceImgs{Nsen+1}{val}{sub} = strvcat(D.inv{val}.contrast.sfname{:});
1623 end
1624
1625 save(fullfile(D.path, D.fname), 'D');
1626 end % end val
1627
1628 sennam{Nsen+1} = 'fused';
1629 fmrigrppreavgnam{Nsen+1,sub} = fullfile(D.path,D.fname);
1630 end
1631 end
1632 end
1633
1634
1635 % Review model-evidences (note cannot compare fused to individual modality
1636 % inversions, since data differs)
1637 for sub = 1:Nsub
1638 subnum = dosubs(sub);
1639 disp(sprintf('\tSub %d (%d)...',subnum,sub))
1640 for val=1:Nval
1641 disp(sprintf('\t\tModel-evidences (a.u.): %s',mat2str(round(squeeze(fmri_grp_mod_evd(sub,val,:))))))
1642 end
1643 end
1644
1645 %%%%%%%%%%%%%%%%% Quick look at mean over subjects...
1646
1647 clear S; S.con = [1 -1]; S.boi{1} = [140 190]; S.boi{2} = []; % boi{2}=[] is all freqs
1648 S.fig = lastfig;
1649 for val=1:Nval
1650 S.val = val;
1651 for sen = find(invertedflag)
1652 S.fig = S.fig+1;
1653 S.P = strvcat(fmrigrppreavgnam{sen,:});
1654 meg_inv_display_con(S);
1655 drawnow
1656 end
1657 end
1658 lastfig = S.fig;
1659
1660 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1661 %%%%%%%%%%%%%%%%% 3D SPMs
1662 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1663
1664 stdir = fullfile(wd,'fMRIGrpSourceSPMs');
1665 try cd(stdir); catch eval(sprintf('!mkdir %s',stdir)); end
1666
1667 clear grpsubs imgfiles grpcons
1668 grpsubs{1} = [1:Nsub];
1669 % If want to split subjects, eg into two groups, eg:
1670 % grpsubs{1} = [1:4]; grpsubs{2} = [5:Nsub];
1671 grpcons = [1:2];
1672 Ngrp = length(grpsubs);
1673 for val=1:Nval
1674 valdir = fullfile(stdir,sprintf('Inv%d',val));
1675 try cd(valdir); catch eval(sprintf('!mkdir %s',valdir)); end
1676 for sen = find(invertedflag)
1677 outdir = fullfile(valdir,sennam{sen})
1678 for grp = 1:Ngrp
1679 for sub = 1:length(grpsubs{grp});
1680 subnum = grpsubs{grp}(sub);
1681 imgfiles{grp}{sub} = FmriGrpSourceImgs{sen}{val}{subnum}(grpcons,:);
1682 end
1683 end
1684 meg_batch_anova(imgfiles,outdir,'/imaging/local/spm/spm5/EEGtemplates/mni_mesh_cortex_8196.nii');
1685 end
1686 end
1687 cd(stdir)
1688 disp('!!!Review the group+fMRI 3D stats if you wish!!!')
1689
1690 % Again you need to know a bit about the SPM Results interface, but press
1691 % Results and select an SPM.mat file for one of the modalities under MSP
1692 % (GS) and enter a new T-contrast [1 -1] to find voxels in which greater source
1693 % amplitude for faces than scrambled faces at 0.001 uncorrected. MSP results
1694 % are not changed much by the addition of fMRI priors (one possibility
1695 % being that the MSP are already optimal; Henson et al, submitted). The
1696 % fMRI priors have a more noticeable effect on the IID results, moving
1697 % them closer to the fMRI priors.
1698
1699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1700 %%%%%%%%%%%%%%%%% 3D PPMs
1701 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1702
1703 % This uses SPMs above to calculate Posterior Probability Maps (PPMs),
1704 % which do not have any mutliple comparison issues (eg no need for RFT)
1705 % However, it will take a long time...
1706
1707 %for val=1:Nval
1708 % valdir = fullfile(stdir,sprintf('Inv%d',val));
1709 % cd(valdir);
1710 % for sen = find(invertedflag)
1711 % outdir = fullfile(valdir,sennam{sen})
1712 % cd(outdir)
1713 % load('SPM.mat');
1714 % spm_spm_Bayes(SPM);
1715 % end
1716 %end
1717 %disp('!!!Review the group+fMRI 3D PPMs if you wish!!!')
1718
1719 % Again you need to know a bit about the PPM Results interface, but press
1720 % Results and select an SPM.mat file for one of the modalities and enter
1721 % a new T-contrast [1 -1], but this time, select "Bayesian" rather than
1722 % "classical" when prompted. Then change the default effect size threshold
1723 % to 0, and probability threshold to 0.99. This will show voxels with a
1724 % >99% probability of the increase for faces relative to scrambled being
1725 % greater than zero.
1726
1727 % Phew! You got to here! Well done!
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.