%
% This is a tutorial on scripting an analysis pipeline with SPM in
% Matlab. It follows the SPM course 2018 held in Institute of Systems
% Neuroscience at UKE.
%
% Summary of the pipeline:
% In short, we convert both anatomical and functional dicom images to NifTi
% and merge them into 4D files. Surface representations from anatomical
% images are computed with cat12 toolbox, which also generates deformation
% fields and segmented tissues. Using the segmented tissues we
% skull-strip the anatomical image. Functional images are motion corrected,
% coregistered to skull-striped anatomical image and resliced. The mean EPI
% image that is written during motion correction is segmented and the
% resulting tissue images are used to create a binary mask in native space.
% This mask is used during first-level analysis. The resulting beta images
% are normalized with the deformation fields coming from the
% Surface/Segment routing of the anatomical images and smoothed for second
% level analysis.
%
% Audience;
% People with a basic understanding of programming in Matlab. You need to
% know what is a variable, structure, cell array etc.
%
% SPM version:
% I have used for this spm12 (version: 6685). I suppose any spm12 version
% would be OK to use it. Don't forget to add spm12 to Matlab's path.
%
% Operating System:
% If you are on Windows, for path definitions you will have to either change
% by hand (or with Find and Replace) the / slashes to back slashes \. Or
% alternatively use the filesep function Matlab provides, try
% sprintf('%s',filesep) to see.
%
% Parallelization:
% In many situations, you can take advantage of parallelizing batches
% across multiple CPUs. Currently, this is not optimal: Replicing the for
% loops with parfor loops might or might not always work. If you are
% insisting on this, you may want to compute your matlabbatch{n} variables
% inside a normal for loop, and re-run a second loop using parfor, example:
% parfor k = 1:N;spm_jobman('run', matlabbatch{k});end
%
% Perspective:
% Most of the code that you can read below is simply snippets that are
% taken from our group's analysis and preprocessing pipeline, which I wrote
% with object-oriented programming. If you master batching SPM pipelines
% with the code below, you may find it useful to have a look at my toolbox
% FancyCarp at https://github.com/selimonat/fancycarp
%
% This analysis script assumes that your data is organized as
% project/data/subjectX/runY. For future, it is best to align with the BIDS
% standards.
%
%
Declare Project-level variables that we will use later.
total_subject = 2; %number of participants.
total_run = 3; %number of runs
project_folder = '/home/onat/Desktop/spm-kurs/my_project/data/'; %location of the project data.
dummies_to_delete = 4; %number of dummy scans to delete
TR = 0.99; %in unit of seconds
tpm_dir = '/home/onat/Documents/Code/Matlab/spm12-6685/tpm/'; % path to TPM templates
path_tpm = @(n) sprintf('%s/TPM.nii,%i',tpm_dir,n) % this is an inline function that returns the Nth TPM template
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% PREPROCESSING OF ANATOMICAL IMAGES%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(1.1) DICOM CONVERSION FOR ANATOMICAL IMAGES
%this will convert all dicom files and generate one NifTi file called
%data.nii.
matlabbatch = [];
for subject_number = 1:total_subject
fprintf('Running Subject: %03d\n',subject_number)
%folder where dicom files are located.
destination = sprintf('%ssub%03d/run000/mrt/',project_folder,subject_number)
%find all files in destimation that starts with MR prefix (you may want
%to change MR in line with your needs.
anatomical_files = spm_select('FPListRec',destination,'^MR');
%create batch cell array
matlabbatch{1}.spm.util.import.dicom.data = cellstr(anatomical_files);
matlabbatch{1}.spm.util.import.dicom.root = 'flat';
matlabbatch{1}.spm.util.import.dicom.outdir = {destination};
matlabbatch{1}.spm.util.import.dicom.protfilter = '.*';
matlabbatch{1}.spm.util.import.dicom.convopts.format = 'nii';
matlabbatch{1}.spm.util.import.dicom.convopts.icedims = 0;
%run the batch: will create sTRIO files.
spm_jobman('run', matlabbatch);
%delete dicom files
files = cellstr(anatomical_files);
delete(files{:});
% rename sTRIO... file to data.nii
files = spm_select('FPListRec',destination,'^sTRIO');
movefile(files,regexprep(files,sprintf('%ssTRIO.*$',filesep),sprintf('%sdata.nii',filesep)));%rename it to data.nii
end
(1.2) Runs CAT12 Segment Surface routine.
If you would like to do the conversion via the GUI, you need to first run the CAT12 toolbox from the SPM window, SPM GUI > toolboxes > CAT12 > Preprocessing > Segment Data.
The results of this step are needed for skull stripping.
This is a time consuming part, you may want to parallelize the computation if you have the parallel computing toolbox. I tested this in my computer using parfor, and it worked.
Think thoroughly if you need this step, as it takes considerable time.
Will write to the disk: mri and report folders mri/iy_data.nii mri/p1data.nii mri/p2data.nii mri/wmdata.nii mri_y_data.nii
for subject_number = 1:total_subject
%OR: parfor subject_number = 1:total_subject
fprintf('Running Subject: %03d\n',subject_number)
path_hr = sprintf('%ssub%03d/run000/mrt/data.nii',project_folder,subject_number);
matlabbatch = [];
matlabbatch{1}.spm.tools.cat.estwrite.data = {spm_select('expand',path_hr)};
matlabbatch{1}.spm.tools.cat.estwrite.nproc = 0;
matlabbatch{1}.spm.tools.cat.estwrite.opts.tpm = {sprintf('%sTPM.nii',tpm_dir)};
matlabbatch{1}.spm.tools.cat.estwrite.opts.affreg = 'mni';
matlabbatch{1}.spm.tools.cat.estwrite.extopts.APP = 1;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.LASstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.gcutstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.cleanupstr = 0.5;
matlabbatch{1}.spm.tools.cat.estwrite.extopts.darteltpm = {'/home/onat/Documents/Code/Matlab/spm12-6685/toolbox/cat12/templates_1.50mm/Template_1_IXI555_MNI152.nii'};
matlabbatch{1}.spm.tools.cat.estwrite.extopts.vox = 1.5;
matlabbatch{1}.spm.tools.cat.estwrite.output.surface = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.native = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.GM.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.native = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.mod = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.WM.dartel = 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.bias.warped = 1;
matlabbatch{1}.spm.tools.cat.estwrite.output.jacobian.warped= 0;
matlabbatch{1}.spm.tools.cat.estwrite.output.warps = [1 1];
%
spm_jobman('run', matlabbatch);
end
(1.3) Skull Strip anatomical image
Needs results of SegmentSurface, will produce a skullstripped version of hr (ss_data.nii). It will also automatically create a normalized version (w_ss_data.nii).
%path to p1 and p2 images created by SegmentSurface
for subject_number = 1:total_subject
fprintf('Running Subject: %03d\n',subject_number)
path_hr = sprintf('%ssub%03d/run000/mrt/data.nii', project_folder,subject_number);
dir_hr = sprintf('%ssub%03d/run000/mrt/', project_folder,subject_number);
c1 = sprintf('%ssub%03d/run000/mrt/mri/p1data.nii',project_folder,subject_number);
c2 = sprintf('%ssub%03d/run000/mrt/mri/p2data.nii',project_folder,subject_number);
path_skullstrip = sprintf('%ssub%03d/run000/mrt/ss_data.nii', project_folder,subject_number);
deformation_field = sprintf('%ssub%03d/run000/mrt/mri/y_data.nii',project_folder,subject_number);
if exist(c1) && exist(c2) %only strip if C1 and C2 images are present.
matlabbatch = [];
matlabbatch{1}.spm.util.imcalc.input = cellstr(strvcat(path_hr,c1,c2)); %3 volumes that are given as arguments to the calculator
matlabbatch{1}.spm.util.imcalc.output = path_skullstrip;
matlabbatch{1}.spm.util.imcalc.outdir = {dir_hr};
matlabbatch{1}.spm.util.imcalc.expression = 'i1.*((i2+i3)>0.2)'; % operations to be conducted with the 3 images
matlabbatch{1}.spm.util.imcalc.options.dmtx = 0;
matlabbatch{1}.spm.util.imcalc.options.mask = 0;
matlabbatch{1}.spm.util.imcalc.options.interp = 1;
matlabbatch{1}.spm.util.imcalc.options.dtype = 4;
%spm_jobman('run', matlabbatch);
%normalize the skull stripped image as soon as it is
%created.
VolumeNormalize(deformation_field,path_skullstrip);
else
fprintf('Need to run segment first...\n')
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% PREPROCESSING OF FUNCTIONAL IMAGES %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(2.1) DICOM CONVERSION FOR FUNCTIONAL IMAGES and REMOVAL OF DUMMY SCANS
%this will convert all dicom files, delete dummies and generate a 4D NifTi
%file called data.nii.
matlabbatch = [];
for subject_number = 1:total_subject
for run_number = 1:total_run
fprintf('Running Subject: %03d, Run: %03d\n',subject_number,run_number)
% find dicom files in the destination folder
destination = sprintf('%ssub%03d/run%03d/mrt/',project_folder,subject_number,run_number);
% you may need to change ^MR according to your needs.
dicom_files = spm_select('FPListRec',destination,'^MR');
% create the batch.
matlabbatch = [];
matlabbatch{1}.spm.util.import.dicom.data = cellstr(dicom_files);
matlabbatch{1}.spm.util.import.dicom.root = 'flat';
matlabbatch{1}.spm.util.import.dicom.outdir = {destination};
matlabbatch{1}.spm.util.import.dicom.protfilter = '.*';
matlabbatch{1}.spm.util.import.dicom.convopts.format = 'nii';
matlabbatch{1}.spm.util.import.dicom.convopts.icedims = 0;
% run batch: will create fTRIO..nii files.
spm_jobman('run', matlabbatch);
% delete first N dummies.
fprintf('Deleting dummy files\n')
files = spm_select('FPListRec',destination,'^fTRIO');%change fTRIO if your files has a different prefix
files_to_delete = cellstr(files(1:dummies_to_delete,:));
delete(files_to_delete{:});
% merge to 4D
fprintf('Merging fTRIO files into 4D\n')
files = spm_select('FPListRec',destination,'^fTRIO');%change fTRIO if your files has a different prefix
matlabbatch = [];
matlabbatch{1}.spm.util.cat.vols = cellstr(files);
matlabbatch{1}.spm.util.cat.name = 'data.nii';
matlabbatch{1}.spm.util.cat.dtype = 0;
% run batch: will create a file called data.nii
spm_jobman('run', matlabbatch);
% delete all dicom files
fprintf('Deleting dicom images files\n')
dicom_files = cellstr(dicom_files);
delete(dicom_files{:});
% delete all 3D fTRIO files
fprintf('Deleting fTRIOimages files\n')
files = cellstr(files);
delete(files{:});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(2.2) Realignment and coregistration
We will make a double pass realignment and save only the mean image. This will write rp_data.txt files on the disk for every run. We will than find the best transformation of this mean image to the anatomical one. And finally apply this transformation to all other EPIs and reslice them for further analysis. This will create r_data.nii image.
for subject_index = 1:total_subject
%
fprintf('Running Subject: %03d, Run: %03d\n',subject_index)
% get the EPIs from all runs for this subject
epi_run = [];
c = 0;
for run_index = 1:total_run
c = c +1;
path_epi = sprintf('%ssub%03d/run%03d/mrt/data.nii',project_folder,subject_index,run_index)
files = spm_select('expand',path_epi);
epi_run{c} = cellstr(files);
end
%
mean_epi = sprintf('%ssub%03d/run001/mrt/meandata.nii',project_folder,subject_index);
anatomical = sprintf('%ssub%03d/run000/mrt/data.nii',project_folder,subject_index);
%or this: if you had already an skull striped anatomical image:"
%anatomical = sprintf('%ssub%03d/run000/mrt/ss_data.nii',project_folder,subject_index);
matlabbatch = [];
% double-pass realign EPIs and reslice the mean image only.
matlabbatch{1}.spm.spatial.realign.estwrite.data = epi_run;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.quality = 0.9;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.sep = 4;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.fwhm = 5;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.rtm = 1;%double pass
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.interp = 2;
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.weight = '';
% we would like to have the mean_image saved for coreg
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.which = [0 1];%write only the mean image for later
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.interp = 4;
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.mask = 1;
matlabbatch{1}.spm.spatial.realign.estwrite.roptions.prefix = 'r';
% coregister EPIs to anatomical (only the affine matrix is modified)
% this procedure works best when the skull stripped image is used as
% anatomical reference. You may want to run that batch first below, and
% then adjust above the path to anatomical by adding the ss_ prefix.
matlabbatch{2}.spm.spatial.coreg.estimate.ref = cellstr(anatomical); %the one that is kept immobile
matlabbatch{2}.spm.spatial.coreg.estimate.source = cellstr(mean_epi); %the one that is moved around to best fit
matlabbatch{2}.spm.spatial.coreg.estimate.other = vertcat(epi_run{:}); %all others that are moved around using the best fitting transformation.
matlabbatch{2}.spm.spatial.coreg.estimate.eoptions.cost_fun = 'nmi';
matlabbatch{2}.spm.spatial.coreg.estimate.eoptions.sep = [4 2];
matlabbatch{2}.spm.spatial.coreg.estimate.eoptions.tol = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
matlabbatch{2}.spm.spatial.coreg.estimate.eoptions.fwhm = [7 7];
% write EPIs
matlabbatch{3}.spm.spatial.realign.write.data = vertcat(epi_run{:});
matlabbatch{3}.spm.spatial.realign.write.roptions.which = [2 1];%all images as well as the mean image.
matlabbatch{3}.spm.spatial.realign.write.roptions.interp = 4;
matlabbatch{3}.spm.spatial.realign.write.roptions.wrap = [0 0 0];
matlabbatch{3}.spm.spatial.realign.write.roptions.mask = 1;
matlabbatch{3}.spm.spatial.realign.write.roptions.prefix = 'r';
spm_jobman('run', matlabbatch);
end
to have a look at the motion parameters change to the first run of participant N
a = load('rp_data.txt')
subplot(1,2,1)
plot(a(:,1:3))
subplot(1,2,2)
plot(a(:,4:6))
(2.3) Skull Strip mean-EPI image but first segment mean EPI
Same as what we did with the anatomical image. The skull stripped mean image will be handy to use it as a mask for the first-level analysis. Will make it faster. However to compute a skull strip mean EPI we need the results of segmentation. We will not run this time CAT12's segmentation, but instead spm's default segmentation.
(2.3A) Runs the new segment of SPM12 on the mean EPI image.
Will write to the disk: meandata_seg8.mat iy_meandata.nii c{1-5}meandata.nii y_meandata.nii
for subject_index = 1:total_subject
path_meanepi = sprintf('%ssub%03d/run001/mrt/meandata.nii',project_folder,subject_index);
matlabbatch = [];
matlabbatch{1}.spm.spatial.preproc.channel.vols = cellstr(path_meanepi);
matlabbatch{1}.spm.spatial.preproc.channel.biasreg = 0.001;
matlabbatch{1}.spm.spatial.preproc.channel.biasfwhm = 60;
matlabbatch{1}.spm.spatial.preproc.channel.write = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(1).tpm = {path_tpm(1)};
matlabbatch{1}.spm.spatial.preproc.tissue(1).ngaus = 1;
matlabbatch{1}.spm.spatial.preproc.tissue(1).native = [1 0];
matlabbatch{1}.spm.spatial.preproc.tissue(1).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(2).tpm = {path_tpm(2)};
matlabbatch{1}.spm.spatial.preproc.tissue(2).ngaus = 1;
matlabbatch{1}.spm.spatial.preproc.tissue(2).native = [1 0];
matlabbatch{1}.spm.spatial.preproc.tissue(2).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(3).tpm = {path_tpm(3)};
matlabbatch{1}.spm.spatial.preproc.tissue(3).ngaus = 2;
matlabbatch{1}.spm.spatial.preproc.tissue(3).native = [1 0];
matlabbatch{1}.spm.spatial.preproc.tissue(3).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(4).tpm = {path_tpm(4)};
matlabbatch{1}.spm.spatial.preproc.tissue(4).ngaus = 3;
matlabbatch{1}.spm.spatial.preproc.tissue(4).native = [1 0];
matlabbatch{1}.spm.spatial.preproc.tissue(4).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(5).tpm = {path_tpm(5)};
matlabbatch{1}.spm.spatial.preproc.tissue(5).ngaus = 4;
matlabbatch{1}.spm.spatial.preproc.tissue(5).native = [1 0];
matlabbatch{1}.spm.spatial.preproc.tissue(5).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(6).tpm = {path_tpm(6)};
matlabbatch{1}.spm.spatial.preproc.tissue(6).ngaus = 2;
matlabbatch{1}.spm.spatial.preproc.tissue(6).native = [0 0];
matlabbatch{1}.spm.spatial.preproc.tissue(6).warped = [0 0];
matlabbatch{1}.spm.spatial.preproc.warp.mrf = 1;
matlabbatch{1}.spm.spatial.preproc.warp.cleanup = 1;
matlabbatch{1}.spm.spatial.preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
matlabbatch{1}.spm.spatial.preproc.warp.affreg = 'mni';
matlabbatch{1}.spm.spatial.preproc.warp.fwhm = 0;
matlabbatch{1}.spm.spatial.preproc.warp.samp = 3;
matlabbatch{1}.spm.spatial.preproc.warp.write = [1 1];
spm_jobman('run', matlabbatch);
end
(2.3B) skullstrip the mean EPIs, creating a mask in native space
for subject_index = 1:total_subject
path_meanepi_segmented = ''
for n = 1:3
path_meanepi_segmented = strvcat(path_meanepi_segmented, sprintf('%ssub%03d/run001/mrt/c%dmeandata.nii',project_folder,subject_index,n));
end
Vout = sprintf('%ssub%03d/run001/mrt/ss_meandata.nii',project_folder,subject_index); %output volume that will be written to disk
spm_imcalc(path_meanepi_segmented,Vout,'(i1+i2+i3)>0')
end
%check out the results using Mango, you must see a nice binary volume in native space
%that we can use for masking first level analysis.
(2.3) First-level analysis
This section relies on the data that is stored in runYYY/design/modelXX
folder. Each modelXX folder contains one specific model. Organizing
models in separate folders will help you encapsulate results
in one folder specific for one analysis. This is because spm will write
the results of an analysis using modelXX in a folder names spm/modelXX/
What is contained in the modelXX is folder is data.mat file that
contains information about stimulus onsets. When you load data.mat to
the workspace, this spawns a structure array called cond. The fields
names of the cond are organized exactly the same way SPM expects onset
information. Therefore, this Matlab structure can be directly fed into
the matlabbatch variable without further massaging.
What is important also is that the motion parameters that we computed
above will be inserted as nuisance covariates.
This batch will write the SPM.mat file together with beta images.
model_num = 1;
path_skullstrip_meanepi = sprintf('%ssub%03d/run001/mrt/ss_meandata.nii',project_folder,subject_index);
for subject_number = 1:total_subject
spm_dir = sprintf('%ssub%03d/run001/spm/model%02d/',project_folder,subject_number,model_num);
path_spm = sprintf('%ssub%03d/run001/spm/model%02d/SPM.mat',project_folder,subject_number,model_num);
deformation_field = sprintf('%ssub%03d/run000/mrt/mri/y_data.nii',project_folder,subject_number);
matlabbatch = [];
matlabbatch{1}.spm.stats.fmri_spec.dir = {spm_dir};
matlabbatch{1}.spm.stats.fmri_spec.timing.units = 'scans';
matlabbatch{1}.spm.stats.fmri_spec.timing.RT = TR;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t = 16;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = 1;
for session = 1%:total_run
%location of EPIs
path_epi = sprintf('%ssub%03d/run%03d/mrt/rdata.nii',project_folder,subject_number,session);
matlabbatch{1}.spm.stats.fmri_spec.sess(session).scans = cellstr(spm_select('expand',path_epi));%use always the realigned data.
%location of onsets
path_model = sprintf('%ssub%03d/run%03d/design/model%02d/data.mat',project_folder,subject_number,session,model_num);
dummy = load(path_model);
matlabbatch{1}.spm.stats.fmri_spec.sess(session).cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(session).cond = dummy.cond;
%load nuissance parameters
path_param_motion = sprintf('%ssub%03d/run%03d/mrt/rp_data.txt',project_folder,subject_number,session);
nuis = load(path_param_motion);
nuis = zscore([nuis [zeros(1,size(nuis,2));diff(nuis)] nuis.^2 [zeros(1,size(nuis,2));diff(nuis)].^2 ]);
for nNuis = 1:size(nuis,2)
matlabbatch{1}.spm.stats.fmri_spec.sess(session).regress(nNuis).val = nuis(:,nNuis);
matlabbatch{1}.spm.stats.fmri_spec.sess(session).regress(nNuis).name = mat2str(nNuis);
end
%
matlabbatch{1}.spm.stats.fmri_spec.sess(session).multi = {''};
matlabbatch{1}.spm.stats.fmri_spec.sess(session).multi_reg = {''};
matlabbatch{1}.spm.stats.fmri_spec.sess(session).hpf = 128;
end
matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [0 0];%we have [0 0], [ 1 0] or [ 1 1] for 1, 2, or 3 regressors.
matlabbatch{1}.spm.stats.fmri_spec.volt = 1;
matlabbatch{1}.spm.stats.fmri_spec.global = 'None';
matlabbatch{1}.spm.stats.fmri_spec.mthresh = -Inf;%
matlabbatch{1}.spm.stats.fmri_spec.mask = {path_skullstrip_meanepi}; %use the native mask we computed above.
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'none';
spm_jobman('run', matlabbatch);%create SPM file first
%estimation
matlabbatch = [];
matlabbatch{1}.spm.stats.fmri_est.spmmat = {path_spm};
matlabbatch{1}.spm.stats.fmri_est.method.Classical = 1;
spm_jobman('run', matlabbatch);
% As you must have remarked, we did not yet smoothened the EPI images,
% which is slightly odd with typical processing pipelines.
% The idea is that smoothing is a linear operator. And linear operators are
% commutative, whether you do first A and then B, is the same as B and then
% A if A and B operations are linear. Given that modelling and smoothing
% are linear operations, we can first fit a model and then to the smoothing
% and normalization for the second-level. This saves tremendous amounts of
% disk space.
% Have a look at the following thread in the SPM maillist:
%https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind00&L=spm&P=R30456&1=spm&9=A&I=-3&J=on&X=3C2B8009325E6AC405&Y=sonat%40uke.de&d=No+Match%3BMatch%3BMatches&z=4
% normalize, apply deformation field to beta images.
total_beta = 10;% Bulk of beta images are about nuissance images, which we are not interesting.
% we will smooth and normalize only the
% first TOTAL_BETA images. This will
% different depending on your model.
beta_weights = '';
for n = 1:total_beta
beta_weights = strvcat(beta_weights , sprintf('%sbeta_%04d.nii',spm_dir,n));
end
VolumeNormalize(deformation_field,beta_weights)
% and smooth.
beta_weights = '';
for n = 1:total_beta
beta_weights = strvcat(beta_weights , sprintf('%sw_beta_%04d.nii',spm_dir,n));
end
VolumeSmooth(beta_weights,8);
end
(2.4) Second-level analysis using ANOVA
2nd level analysis only makes sense with smoothened and normalized images, thus we seek for beta images with a prefix s_w_.
The result of this analysis are stored in project_path/spm/modelXX in the form of beta images. You can then select the SPM.mat file from SPM GUI and make statistical analyses using the contrast GUI.
total_beta = 10;
model_num = 1;
spm_dir = sprintf('%sseconlevel/model%02d/',project_folder,model_num);
matlabbatch = [];
for beta_index = 1:total_beta
%collect the beta image for a given condition for all participants
beta_images = '';
for subject_index = 1:total_subject
beta_images = strvcat(beta_images,sprintf('%ssub%03d/run001/spm/model%02d/s_w_beta_%04d.nii',project_folder,subject_number,model_num,beta_index));
end
%and store it in an ANOVA cell.
matlabbatch{1}.spm.stats.factorial_design.des.anova.icell(beta_index).scans = cellstr(beta_images );
end
%
matlabbatch{1}.spm.stats.factorial_design.dir = cellstr(spm_dir);
matlabbatch{1}.spm.stats.factorial_design.des.anova.dept = 0;
matlabbatch{1}.spm.stats.factorial_design.des.anova.variance = 1;
matlabbatch{1}.spm.stats.factorial_design.des.anova.gmsca = 0;
matlabbatch{1}.spm.stats.factorial_design.des.anova.ancova = 0;
matlabbatch{1}.spm.stats.factorial_design.cov = struct('c', {}, 'cname', {}, 'iCFI', {}, 'iCC', {});
matlabbatch{1}.spm.stats.factorial_design.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {});
matlabbatch{1}.spm.stats.factorial_design.masking.tm.tm_none = 1;
matlabbatch{1}.spm.stats.factorial_design.masking.im = 1;
matlabbatch{1}.spm.stats.factorial_design.masking.em = {''};
matlabbatch{1}.spm.stats.factorial_design.globalc.g_omit = 1;
matlabbatch{1}.spm.stats.factorial_design.globalm.gmsca.gmsca_no = 1;
matlabbatch{1}.spm.stats.factorial_design.globalm.glonorm = 1;
%
matlabbatch{2}.spm.stats.fmri_est.spmmat = {[spm_dir '/SPM.mat']};
matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1;
%
spm_jobman('run', matlabbatch);
