Contents
% First level (individual subject) analysis % Christian Kaul, July 2007 % presented @ Matlab for Cogntive Neuroscience 12/07/07 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % this script assumes a folder structure like this: % % exp_name/ % exp_name/s1 % exp_name/s1/structural % the structural folder contains % structural, fieldmaps, segmentation,... % exp_name/s1/structural/ROI %the ROI folder contains all the area defs, ergo put your amygdala def here % % exp_name/s1/localizer_data % exp_name/s1/localizer_data/run1 % contains functional files % exp_name/s1/localizer_data/run2 % exp_name/s1/localizer_stats % contains the SPM, betas, con-images % % exp_name/s1/experiment % exp_name/s1/experiment/run1 % contains functional files % exp_name/s1/experiment/run2 % exp_name/s1/experiment_stats % contains the SPM, betas, con-images % this batch script with perform the following: % 1) setting up a SPM model with a number of conditions, onsets and % durations % 2) estimate this model % 3) set up contrasts within the estimated SPM model % clear all % Specify root directory % The intended output (SPM results) directory for each subject is statsDirName = 'exp_stats';
Specify root directory
origDir = 'D:\V1load_motion2'; dataFolder = 'experiment\run'; structuralFolder = 'structural'; % subjects and the number of volumes collected for each run and the name of % the structural file subjects = { ... '1', ones(1,48)*48-3, 'ED-MA04903.img'; ... '2', ones(1,48)*48-3, 'GV-MA04904.img'; ... '3', ones(1,48)*48-3, 'HR-sMA04830.img'; ... '4', ones(1,48)*48-3, 'DB-meansMA04771.img'; ... '5', ones(1,48)*48-3, 'AD.img'; ... }; % ========================================================================= % ========================================================================= % ========================================================================= % =========================== Set up model ================================= % ========================================================================= TR = 1.3; hpcutoff = 128; N_con = {'task_and_motion', 'fixation'}; % conditions of MT localizer % ========================================================================= cd(origDir); spm fmri for s0 = 1 : size(subjects,1) runs = length(subjects{s0,2}); % go to subjects directory cd([origDir '\s' subjects{s0,1}]); % if the intended stats directory does not exits, make it and go into it if exist(statsDirName) ~= 7 mkdir(statsDirName); cd(statsDirName); else % if it does exist, go into it and delete its contents. cd(statsDirName); delete('*.*'); end % setup job structure for model specification % --------------------------------------------------------------------- % timing variables jobs{1}.stats{1}.fmri_spec.timing.units = 'scans'; jobs{1}.stats{1}.fmri_spec.timing.RT = TR; jobs{1}.stats{1}.fmri_spec.timing.fmri_t = 16 jobs{1}.stats{1}.fmri_spec.timing.fmri_t0 = 1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for rnum = 1:runs % spm Directory <-- intended output directory spmDir = [origDir '\s' subjects{s0,1} '\' statsDirName]; jobs{1}.stats{1}.fmri_spec.dir = {spmDir}; % specify number of volumes in the directory nscans = subjects{s0,2}(rnum); % define epi files in the run epiDir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(rnum)]; % select scans and assign to job f = spm_select('List', epiDir, '^rf.*\.img$'); fs = cellstr([repmat([epiDir '\'],size(f,1),1) f]); jobs{1}.stats{1}.fmri_spec.sess(rnum).scans = fs; % clear temporary variables for next run f = []; fs = []; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % You should define your conditions beforehand --- % % See note from SPM5 help attached to the END OF THIS DOCUMENT for % % an explanation of what the condition file is and how you should % % construct it. % exemplar structure for 2 conditions in each run if subjects{s0,2}(rnum) == 45 % this if ensures that the right runs are selected onsets(1,:) = [4]; % actual start condition 1 onsets(2,:) = [28]; % actual start condition 2 (rest period) runDuration = [24 15]; % in scans (can be seconds if spm is set up like that) % elseif subjects{s0,2}(rnum) == 200 % if different run types this % is an option to distinguish between them end for c = 1 : length(N_con) jobs{1}.stats{1}.fmri_spec.sess(rnum).cond(c).name = N_con{c}; jobs{1}.stats{1}.fmri_spec.sess(rnum).cond(c).onset = onsets(c,:); jobs{1}.stats{1}.fmri_spec.sess(rnum).cond(c).duration = ones(size(onsets(c,:))) * runDuration(c); jobs{1}.stats{1}.fmri_spec.sess(rnum).cond(c).tmod = 0; % no time modulation jobs{1}.stats{1}.fmri_spec.sess(rnum).cond(c).mod = struct('name', {}, 'param', {}, 'poly', {}); end onsets = []; runDuration = []; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % motion correction parameters % first find the motion correction text file's full name exp_run_folder_dir = [origDir '\s' subjects{s0,1} '\' dataFolder int2str(rnum)]; mCorrFile = spm_select('List',exp_run_folder_dir,'^rp.*\.txt$'); % then add it to its path and assign it to SPM jobs mCorrFilePath = [exp_run_folder_dir '\' mCorrFile]; jobs{1}.stats{1}.fmri_spec.sess(rnum).multi_reg = {mCorrFilePath}; % factorial design jobs{1}.stats{1}.fmri_spec.fact = struct('name', {}, 'levels', {}); % high pass filter jobs{1}.stats{1}.fmri_spec.sess(rnum).hpf = hpcutoff; end % end of run loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % basis functions: neither time nor dispersion derivative jobs{1}.stats{1}.fmri_spec.bases.hrf.derivs = [0 0]; % model interactions (Volterra) OPTIONS: 1|2 = order of convolution jobs{1}.stats{1}.fmri_spec.volt = 1; % global normalisation jobs{1}.stats{1}.fmri_spec.global = 'None'; % explicit masking jobs{1}.stats{1}.fmri_spec.mask = {''}; % serial correlations jobs{1}.stats{1}.fmri_spec.cvi = 'AR(1)'; % run model specification % --------------------------------------------------------------------- cd(spmDir); % save and run job save specify.mat jobs disp(['RUNNING model specification for subject number ' subjects{s0,1}]); spm_jobman('run','specify.mat'); clear jobs % ========================================================================= % =========================== Estimate ================================= % ========================================================================= % setup job structure for model estimation and estimate % --------------------------------------------------------------------- jobs{1}.stats{1}.fmri_est.spmmat = {[spmDir,'\SPM.mat']}; save estimate.mat jobs disp(['RUNNING model estimation for subject ' subjects{s0,1}]) spm_jobman('run','estimate.mat'); clear jobs end % end of subject loop cd(origDir); %========================================================================== % Attached note from SPM5 help about multiple conditions and how to % construct them: % % If you have multiple conditions then entering the details a % condition at a time is very inefficient. This option can be % used to load all the required information in one go. You will % first need to create a *.mat file containing the relevant % information. % % This *.mat file must include the following cell arrays (each 1 % x n): names, onsets and durations. eg. names=cell(1,5), % onsets=cell(1,5), durations=cell(1,5), then % names{2}='SSent-DSpeak', onsets{2}=[3 5 19 222], % durations{2}=[0 0 0 0], contain the required details of the % second condition. These cell arrays may be made available by % your stimulus delivery program, eg. COGENT. The duration % vectors can contain a single entry if the durations are % identical for all events. % % Time and Parametric effects can also be included. For time % modulation include a cell array (1 x n) called tmod. It should % have a have a single number in each cell. Unused cells may % contain either a 0 or be left empty. The number specifies the % order of time modulation from 0 = No Time Modulation to 6 = % 6th Order Time Modulation. eg. tmod{3} = 1, modulates the 3rd % condition by a linear time effect. % % For parametric modulation include a structure array, which is % up to 1 x n in size, called pmod. n must be less than or equal % to the number of cells in the names/onsets/durations cell % arrays. The structure array pmod must have the fields: name, % param and poly. Each of these fields is in turn a cell array % to allow the inclusion of one or more parametric effects per % column of the design. The field name must be a cell array % containing strings. The field param is a cell array containing % a vector of parameters. Remember each parameter must be the % same length as its corresponding onsets vector. The field poly % is a cell array (for consistency) with each cell containing a % single number specifying the order of the polynomial expansion % from 1 to 6. % % Note that each condition is assigned its corresponding entry % in the structure array (condition 1 parametric modulators are % in pmod(1), condition 2 parametric modulators are in pmod(2), % etc. Within a condition multiple parametric modulators are % accessed via each fields cell arrays. So for condition 1, % parametric modulator 1 would be defined in pmod(1).name{1}, % pmod(1).param{1}, and pmod(1).poly{1}. A second parametric % modulator for condition 1 would be defined as pmod(1).name{2}, % pmod(1).param{2} and pmod(1).poly{2}. If there was also a % parametric modulator for condition 2, then remember the first % modulator for that condition is in cell array 1: % pmod(2).name{1}, pmod(2).param{1}, and pmod(2).poly{1}. If % some, but not all conditions are parametrically modulated, % then the non-modulated indices in the pmod structure can be % left blank. For example, if conditions 1 and 3 but not % condition 2 are modulated, then specify pmod(1) and pmod(3). % Similarly, if conditions 1 and 2 are modulated but there are 3 % conditions overall, it is only necessary for pmod to be a 1 x % 2 structure array. % % EXAMPLE: % Make an empty pmod structure: % pmod = struct('name',{''},'param',{},'poly',{}); % Specify one parametric regressor for the first condition: % pmod(1).name{1} = 'regressor1'; % pmod(1).param{1} = [1 2 4 5 6]; % pmod(1).poly{1} = 1; % Specify 2 parametric regressors for the second condition: % pmod(2).name{1} = 'regressor2-1'; % pmod(2).param{1} = [1 3 5 7]; % pmod(2).poly{1} = 1; % pmod(2).name{2} = 'regressor2-2'; % pmod(2).param{2} = [2 4 6 8 10]; % pmod(2).poly{2} = 1; % % The parametric modulator should be mean corrected if % appropriate. Unused structure entries should have all fields % left empty. % % ========================================================================= % ========================================================================= % ========================================================================= % =========================== Set up contrasts ============================ % ========================================================================= % set up a various number of contrasts clc cd(origDir); spm fmri % contrasts = {'motion_vs_static' ,'motionVSstatic_interaction'}; % contrasts = {'motion_vs_static_2baselines' ,'motionVSstatic_interaction_2baselines'}; contrasts = {'motion_vs_static_2baselines_no_motioncorr', 'motionVSstatic_interaction_2baselines_no_motioncorr'}; for cont_nr = 1:length(contrasts) contrasttype = contrasts{cont_nr}; for s0 = 1 : size(subjects,1) % enter your contast here i.e. like that % contr_input = [-1 2 -1 0 0 0 0 0 0]; contr_input = get_sensible_design_matrix_input(['s' subjects{s0,1}],contrasttype,0); % setup job structure for contrasts to set up % tcontrasts here % ----------------------------------------------------------------- jobs{1}.stats{1}.con.spmmat = {[origDir '\s' subjects{s0,1} '\' statsDirName '\SPM.mat']}; jobs{1}.stats{1}.con.consess{1}.tcon.name = contrasttype; jobs{1}.stats{1}.con.consess{1}.tcon.convec = contr_input; jobs{1}.stats{1}.con.consess{1}.tcon.sessrep = 'none'; % fcontrast here % jobs{1}.stats{1}.con.consess{2}.fcon.name = 'NAME OF CONTRAST'; % jobs{1}.stats{1}.con.consess{2}.fcon.convec = [cell with contrast]; % jobs{1}.stats{1}.con.consess{2}.fcon.sessrep = 'none'; % very important little setting. if 1 then all existing contrasts % are deleted jobs{1}.stats{1}.con.delete = 0; % save and run job save contrasts.mat jobs disp(['RUNNING contrast specification for subject number ' subjects{s0,1}]); spm_jobman('run','contrasts.mat'); disp(['contasts created']); clear jobs end % end of subject loop end % end of contrasts loop cd(origDir);