571 lines
22 KiB
Matlab
571 lines
22 KiB
Matlab
function[job] = spm_bsm(job)
|
|
%
|
|
% Copyright J.B.C. Marsman
|
|
% NeuroImaging Center
|
|
% University Medical Center Groningnen
|
|
%
|
|
% September 2015
|
|
%
|
|
% Cisler 2014
|
|
% Rissman 2004
|
|
% Mumford 2012
|
|
|
|
% Camara 2008 (reward processing)
|
|
|
|
display('== == == Beta Series Method: Running analysis == == ==')
|
|
|
|
N = numel(job.spms);
|
|
|
|
warning('OFF', 'MATLAB:MKDIR:DirectoryExists');
|
|
|
|
for spm_i = 1:N
|
|
spmfile = job.spms{spm_i};
|
|
%% load original SPM file
|
|
load(spmfile);
|
|
|
|
%% check for canonical HRF :
|
|
if (SPM.xBF.order ~=1 )
|
|
warning([sprintf('Problem found in File : %s\n', spmfile) ...
|
|
'No canonical HRF was used to model the data\n' ...
|
|
'This model will be skipped']);
|
|
continue;
|
|
end
|
|
|
|
if (length(SPM.Sess) > 1 )
|
|
warning([sprintf('Problem found in File : %s\n', spmfile) ...
|
|
sprintf('Multiple sessions are not supported.\n') ...
|
|
sprintf('This model will be skipped\n')]);
|
|
continue;
|
|
end
|
|
|
|
SPM_orig = SPM;
|
|
U_orig = SPM.Sess.U;
|
|
|
|
target_folder = [fileparts(spmfile) filesep 'BSM'];
|
|
|
|
%% create output folder
|
|
mkdir(target_folder);
|
|
|
|
%% specify active columns:
|
|
|
|
all_columns = 1:length(SPM.Sess(1).U);
|
|
if ~isempty(job.acfl)
|
|
passive_columns = eval(['all_columns( job.acfl );']);
|
|
else
|
|
passive_columns = [];
|
|
end
|
|
active_columns = setdiff(all_columns, passive_columns);
|
|
|
|
|
|
job.active_columns = active_columns;
|
|
job.SPM_orig = SPM_orig;
|
|
job.U_orig = U_orig;
|
|
job.target_folder = target_folder;
|
|
|
|
%% loop for active columns:
|
|
|
|
methods = {'Traditional BSM Rissman et al', 'Adapted BSM Rissman et al', 'Fast approach (Mumford et al)', 'Adapted fast approach','none'};
|
|
display(sprintf('== == == Beta Series Method: Analysis Type : %s == == ==', methods{job.type}));
|
|
|
|
switch(job.type)
|
|
case 1, %% Traditional BSM Rissman et al
|
|
BSM_type1(job);
|
|
|
|
case 2, %% Adapted Rissman et al
|
|
BSM_type2(job);
|
|
|
|
case 3, %% Fast approach (Mumford et al)
|
|
BSM_type3(job);
|
|
|
|
case 4, %% Adapted Fast approach
|
|
BSM_type4(job);
|
|
|
|
end
|
|
pwd
|
|
keyboard
|
|
save('IntermediateState','job');
|
|
display('== == == Beta Series Method : Concatenating betas == == ==')
|
|
BSM_concat_betas(job)
|
|
return; %for now don't do the last few steps
|
|
display('== == == Beta Series Method : Regressing ROIs== == ==')
|
|
BSM_Regress_Regions(job);
|
|
display('== == == Beta Series Method : Contrasting second level betas == == ==')
|
|
BSM_Contrast_manager(job);
|
|
display('== == == Beta Series Method : Done == == ==')
|
|
end
|
|
|
|
end
|
|
|
|
|
|
function[] = BSM_type1(job)
|
|
%% Traditional BSM by Rissman et al 2004
|
|
|
|
active_columns = job.active_columns;
|
|
|
|
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
column = active_columns(c);
|
|
other_columns = setdiff(1:length(active_columns), column);
|
|
|
|
display(sprintf('*) Exploding column %d...', column));
|
|
|
|
target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
|
|
mkdir(target_folder_c);
|
|
|
|
pos_index =1;
|
|
|
|
U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
|
|
|
|
for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
|
|
|
|
SPM = job.SPM_orig; % copy settings from original SPM model
|
|
target_folder_i = [target_folder_c filesep sprintf('trial_%04d', i)];
|
|
mkdir(target_folder_i);
|
|
SPM.swd = target_folder_i;
|
|
|
|
%% add event i from 'column'
|
|
U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(i) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(i);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
%% add all other events but i from 'column' to next column
|
|
pos_index = pos_index +1;
|
|
other_trials = setdiff(1:length(job.U_orig(column).ons), i);
|
|
|
|
U_new(pos_index).name{1} = sprintf('%s_oc', SPM.Sess.U(column).name{1});
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(other_trials) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(other_trials);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
pos_index = pos_index + 1;
|
|
|
|
for oc = 1:length(other_columns)
|
|
other_column = other_columns(oc);
|
|
|
|
U_new(pos_index).name{1} = 'other conds';
|
|
U_new(pos_index).ons = [U_new(pos_index).ons job.SPM_orig.Sess.U(other_column).ons];
|
|
|
|
if length(job.SPM_orig.Sess.U(other_column).dur) == 1
|
|
U_new(pos_index).dur = [U_new(pos_index).dur repmat(job.SPM_orig.Sess.U(other_column).dur, 1, length(job.SPM_orig.Sess.U(other_column).ons))];
|
|
else
|
|
U_new(pos_index).dur = [U_new(pos_index).dur job.SPM_orig.Sess.U(other_column).dur];
|
|
end
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
|
|
end
|
|
|
|
%% store the new U structure
|
|
SPM.Sess.U = U_new;
|
|
|
|
BSM_finalize_model(job, SPM);
|
|
|
|
end %% end loop for trial i
|
|
end %% end loop for active columns
|
|
end
|
|
|
|
function BSM_type2(job)
|
|
%% Adapted Rissman et al:
|
|
active_columns = job.active_columns;
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
column = active_columns(c);
|
|
other_columns = setdiff(1:length(active_columns), column);
|
|
|
|
display(sprintf('*) Exploding column %d...', column));
|
|
|
|
target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
|
|
mkdir(target_folder_c);
|
|
|
|
pos_index =1;
|
|
|
|
U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
|
|
|
|
for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
|
|
|
|
SPM = job.SPM_orig; % copy settings from original SPM model
|
|
target_folder_i = [target_folder_c filesep sprintf('trial_%04d', i)];
|
|
mkdir(target_folder_i);
|
|
SPM.swd = target_folder_i;
|
|
|
|
%% add event i from 'column'
|
|
U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(i) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(i);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
%% add all other events but i from 'column' to next column
|
|
pos_index = pos_index +1;
|
|
other_trials = setdiff(1:length(job.U_orig(column).ons), i);
|
|
|
|
U_new(pos_index).name{1} = sprintf('%s_oc', SPM.Sess.U(column).name{1});
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(other_trials) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(other_trials);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
pos_index = pos_index + 1;
|
|
|
|
for oc = 1:length(other_columns)
|
|
other_column = other_columns(oc);
|
|
|
|
U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(other_column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(other_column).ons;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(other_column).dur;
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
|
|
|
|
pos_index = pos_index+ 1;
|
|
end
|
|
|
|
%% store the new U structure
|
|
SPM.Sess.U = U_new;
|
|
|
|
BSM_finalize_model(job, SPM);
|
|
end %% end loop for trial i
|
|
|
|
end %% end loop for active columns
|
|
end
|
|
|
|
|
|
|
|
function BSM_type3(job)
|
|
%% Fast approach
|
|
|
|
active_columns = job.active_columns;
|
|
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
column = active_columns(c);
|
|
other_columns = setdiff(1:length(active_columns), column);
|
|
|
|
display(sprintf('*) Exploding column %d...', column));
|
|
|
|
target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
|
|
mkdir(target_folder_c);
|
|
|
|
pos_index =1;
|
|
|
|
SPM = job.SPM_orig;
|
|
SPM.swd = target_folder_c;
|
|
U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
|
|
|
|
for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
|
|
|
|
%% add event i from 'column'
|
|
U_new(pos_index).name{i} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(i) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(i);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
%% add all other events but i from 'column' to next column
|
|
pos_index = pos_index +1;
|
|
end
|
|
|
|
for oc = 1:length(other_columns)
|
|
other_column = other_columns(oc);
|
|
|
|
U_new(pos_index).name{1} = sprintf('%s_%d', SPM.Sess.U(other_column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(other_column).ons;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(other_column).dur;
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
|
|
|
|
pos_index = pos_index+ 1;
|
|
end
|
|
|
|
%% store the new U structure
|
|
SPM.Sess.U = U_new;
|
|
|
|
BSM_finalize_model(job, SPM);
|
|
|
|
%% end loop for active columns
|
|
end
|
|
end
|
|
|
|
function BSM_type4(job)
|
|
active_columns = job.active_columns;
|
|
|
|
%% Adapted Fast approach
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
column = active_columns(c);
|
|
other_columns = setdiff(1:length(active_columns), column);
|
|
|
|
display(sprintf('*) Exploding column %d...', column));
|
|
|
|
target_folder_c = [job.target_folder filesep sprintf('col_%04d', column)];
|
|
mkdir(target_folder_c);
|
|
|
|
pos_index =1;
|
|
|
|
SPM = job.SPM_orig;
|
|
SPM.swd = target_folder_c;
|
|
U_new = struct('name',[],'ons', [], 'dur',[], 'P',[], 'dt', [], 'u', [],'pst',[]);
|
|
|
|
for i = 1:length(job.U_orig(column).ons) %% walk over all onsets within 'column'
|
|
|
|
%% add event i from 'column'
|
|
U_new(pos_index).name{i} = sprintf('%s_%d', SPM.Sess.U(column).name{1} , i);
|
|
U_new(pos_index).ons = job.SPM_orig.Sess.U(column).ons(i) ;
|
|
U_new(pos_index).dur = job.SPM_orig.Sess.U(column).dur(i);
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(column).dt;
|
|
|
|
%% add all other events but i from 'column' to next column
|
|
pos_index = pos_index +1;
|
|
end
|
|
|
|
for oc = 1:length(other_columns)
|
|
other_column = other_columns(oc);
|
|
|
|
U_new(pos_index).name{1} = 'other conds';
|
|
U_new(pos_index).ons = [U_new(pos_index).ons job.SPM_orig.Sess.U(other_column).ons];
|
|
|
|
if length(job.SPM_orig.Sess.U(other_column).dur) == 1
|
|
U_new(pos_index).dur = [U_new(pos_index).dur repmat(job.SPM_orig.Sess.U(other_column).dur, 1, length(job.SPM_orig.Sess.U(other_column).ons))];
|
|
else
|
|
U_new(pos_index).dur = [U_new(pos_index).dur M_orig.Sess.U(other_column).dur];
|
|
end
|
|
U_new(pos_index).P = job.SPM_orig.Sess.U(other_column).P;
|
|
U_new(pos_index).dt = job.SPM_orig.Sess.U(other_column).dt;
|
|
end
|
|
|
|
|
|
%% store the new U structure
|
|
SPM.Sess.U = U_new;
|
|
|
|
BSM_finalize_model(job, SPM);
|
|
|
|
%% end loop for active columns
|
|
end
|
|
end
|
|
|
|
function BSM_finalize_model(job, SPM)
|
|
|
|
U = spm_get_ons(SPM,1);
|
|
SPM.Sess(1).U = U;
|
|
[X,Xn,Fc] = spm_Volterra(SPM.Sess.U,SPM.xBF.bf,SPM.xBF.Volterra);
|
|
|
|
k = SPM.nscan(1);
|
|
fMRI_T = SPM.xBF.T;
|
|
fMRI_T0 = SPM.xBF.T0;
|
|
|
|
X = X((0:(k - 1))*fMRI_T + fMRI_T0 + 32,:);
|
|
% and orthogonalise (within trial type)
|
|
for i = 1:length(Fc)
|
|
X(:,Fc(i).i) = spm_orth(X(:,Fc(i).i));
|
|
end
|
|
|
|
Xname={};
|
|
for i = 1:size(X,2)
|
|
Xname{end + 1} = [sprintf('Sn(%i) ',1) SPM.Sess(1).U(i).name{1}];
|
|
end
|
|
|
|
%% reinstate motion parameters
|
|
if (job.keepr)
|
|
C = job.SPM_orig.Sess(1).C.C;
|
|
Cname = job.SPM_orig.Sess(1).C.name;
|
|
|
|
X = [X spm_detrend(C)];
|
|
Xname = {Xname{:} Cname{:}};
|
|
end
|
|
|
|
|
|
%% reinstate baseline column
|
|
X(:,end+1) = 1;
|
|
SPM.xX.X = X;
|
|
|
|
SPM.xX.iB = size(X,2);
|
|
SPM.xX.iC = 1:size(X,2)-1;
|
|
Xname{end+1} = [sprintf('Sn(%i) ',1) 'constant'];
|
|
SPM.xX.name = Xname;
|
|
|
|
try
|
|
SPM.xX = rmfield(SPM.xX, 'W');
|
|
catch
|
|
end
|
|
save([SPM.swd filesep 'SPM.mat'], 'SPM');
|
|
|
|
%% estimate model:
|
|
matlabbatch{1}.spm.stats.fmri_est.spmmat = {[SPM.swd filesep 'SPM.mat']};
|
|
matlabbatch{1}.spm.stats.fmri_est.method.Classical = 1;
|
|
spm_jobman('run', matlabbatch);
|
|
end
|
|
|
|
%% Concatenate betas per condition in a Beta series 4D file
|
|
function BSM_concat_betas(job)
|
|
active_columns = job.active_columns;
|
|
target_folder = job.target_folder;
|
|
if (job.type < 3)
|
|
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
|
|
column = active_columns(c);
|
|
target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
|
|
beta_files = spm_select('FPListRec', target_folder_c, 'beta_0001.nii');
|
|
|
|
matlabbatch{1}.spm.util.cat.vols = cellstr(beta_files);
|
|
matlabbatch{1}.spm.util.cat.name = [target_folder_c filesep 'betas.nii'];
|
|
matlabbatch{1}.spm.util.cat.dtype = 4;
|
|
spm_jobman('run', matlabbatch);
|
|
|
|
if (job.keepfl == 0)
|
|
|
|
to_delete = spm_select('FPListRec', target_folder_c, '^*.*');
|
|
items = ~regexp(to_delete, 'betas.nii');
|
|
to_delete = to_delete(items);
|
|
spm_unlink(to_delete); %% throw away all files
|
|
end
|
|
end
|
|
|
|
else
|
|
|
|
for c = 1:length(active_columns) %% split out all columns but baseline
|
|
column = active_columns(c);
|
|
beta_for_cond_c = 1:length(job.U_orig(column).ons);
|
|
target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
|
|
beta_files = {}; %% empty beta_files
|
|
for i =1:length(job.U_orig(column).ons)
|
|
beta_files{i} = spm_select('FPList', target_folder_c, sprintf('^beta_%04d.nii', i));
|
|
end
|
|
|
|
matlabbatch{1}.spm.util.cat.vols = beta_files';
|
|
matlabbatch{1}.spm.util.cat.name = [target_folder_c filesep 'betas.nii'];
|
|
matlabbatch{1}.spm.util.cat.dtype = 4;
|
|
spm_jobman('run', matlabbatch);
|
|
|
|
if (job.keepfl == 0)
|
|
display('*) Cleaning up ...');
|
|
to_delete = cellstr(spm_select('FPListRec', target_folder_c, '^*.*'));
|
|
items_to_retain = ~cellfun(@isempty,regexp(to_delete, 'betas'));
|
|
to_delete(find(items_to_retain)) = [];
|
|
|
|
for f = 1:length(to_delete)
|
|
spm_unlink(to_delete{f}); %% throw away all files
|
|
end
|
|
end
|
|
|
|
end
|
|
end
|
|
end
|
|
|
|
%% Step 3: Use Betas series files for GLM on a ROI timecourse (average Beta series extracted for a ROI)
|
|
%% Organized now as follows : / BSM / <ROI_filename> / col_000n
|
|
function[] = BSM_Regress_Regions(job)
|
|
target_folder = job.target_folder;
|
|
active_columns = job.active_columns;
|
|
|
|
for r = 1:length(job.rois)
|
|
ROI(r).volume = spm_read_vols(spm_vol(job.rois{r}));
|
|
ROI(r).ind = find(ROI(r).volume);
|
|
[ROI(r).x ROI(r).y ROI(r).z] = ind2sub(size(ROI(r).volume), ROI(r).ind);
|
|
end
|
|
for r = 1:length(job.rois)
|
|
|
|
if (isempty(ROI(r).ind))
|
|
ROI(r) = [];
|
|
warning([sprintf(' %d ) Problem found for ROI : ''%s'', no active voxels.\n',r, job.rois{r}) ...
|
|
sprintf('This ROI will be skipped\n')]);
|
|
else
|
|
display(sprintf(' %d ) ROI ''%s'' contains %d voxels', r, job.rois{r}, length(ROI(r).ind)));
|
|
end
|
|
|
|
end
|
|
|
|
for r = 1:length(job.rois)
|
|
for c = 1:length(job.active_columns)
|
|
column = active_columns(c);
|
|
target_folder_c = [target_folder filesep sprintf('col_%04d', column)];
|
|
|
|
betas = cellstr(spm_select('ExtFPList', target_folder_c, 'betas.nii', 1:999));
|
|
B = zeros(1, length(betas));
|
|
|
|
for b = 1:length(betas)
|
|
H = spm_vol(betas{b});
|
|
B(b) = mean(spm_sample_vol(H, ROI(r).x, ROI(r).y, ROI(r).z,0));
|
|
end
|
|
|
|
[base roiname ext]= fileparts(job.rois{r});
|
|
target_folder_r = [target_folder filesep roiname filesep sprintf('col_%04d', column)];
|
|
mkdir(target_folder_r);
|
|
|
|
matlabbatch{1}.spm.stats.factorial_design.dir = { target_folder_r};
|
|
matlabbatch{1}.spm.stats.factorial_design.des.t1.scans = betas;
|
|
matlabbatch{1}.spm.stats.factorial_design.cov.c = B;
|
|
matlabbatch{1}.spm.stats.factorial_design.cov.cname = roiname;
|
|
matlabbatch{1}.spm.stats.factorial_design.cov.iCFI = 1;
|
|
matlabbatch{1}.spm.stats.factorial_design.cov.iCC = 1;
|
|
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(1) = cfg_dep('Factorial design specification: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
|
|
matlabbatch{2}.spm.stats.fmri_est.write_residuals = 0;
|
|
matlabbatch{2}.spm.stats.fmri_est.method.Classical = 1;
|
|
|
|
matlabbatch{3}.spm.stats.con.spmmat(1) = cfg_dep('Model estimation: SPM.mat File', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
|
|
matlabbatch{3}.spm.stats.con.consess{1}.tcon.name = 'r';
|
|
matlabbatch{3}.spm.stats.con.consess{1}.tcon.weights = [0 1];
|
|
matlabbatch{3}.spm.stats.con.consess{1}.tcon.sessrep = 'none';
|
|
matlabbatch{3}.spm.stats.con.consess{2}.tcon.name = '-r';
|
|
matlabbatch{3}.spm.stats.con.consess{2}.tcon.weights = [0 -1];
|
|
matlabbatch{3}.spm.stats.con.consess{2}.tcon.sessrep = 'none';
|
|
matlabbatch{3}.spm.stats.con.delete = 0;
|
|
|
|
spm_jobman('run', matlabbatch);
|
|
end
|
|
end
|
|
|
|
end
|
|
|
|
%% Step 4: Contrast Second level Betas for each ROI
|
|
function[] = BSM_Contrast_manager(job)
|
|
target_folder = job.target_folder;
|
|
for r = 1:length(job.rois)
|
|
[base roiname ext]= fileparts(job.rois{r});
|
|
target_folder_final = [target_folder filesep roiname filesep 'results'];
|
|
mkdir(target_folder_final);
|
|
|
|
for c1 = 1:length(job.active_columns)
|
|
for c2 = setdiff(1:length(job.active_columns),c1)
|
|
|
|
betafile_c1 = spm_select('FPList', [target_folder filesep roiname filesep sprintf('col_%04d', c1)], 'beta_0002.nii');
|
|
betafile_c2 = spm_select('FPList', [target_folder filesep roiname filesep sprintf('col_%04d', c2)], 'beta_0002.nii');
|
|
|
|
col1 = job.active_columns(c1);
|
|
col2 = job.active_columns(c2);
|
|
|
|
matlabbatch{1}.spm.util.imcalc.input = { betafile_c1; betafile_c2 };
|
|
matlabbatch{1}.spm.util.imcalc.output = sprintf('con_%04d_minus_%04d.nii', col1, col2);
|
|
matlabbatch{1}.spm.util.imcalc.outdir = { target_folder_final};
|
|
matlabbatch{1}.spm.util.imcalc.expression = 'i1 - i2';
|
|
matlabbatch{1}.spm.util.imcalc.var = struct('name', {}, 'value', {});
|
|
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;
|
|
|
|
matlabbatch{2}.spm.util.imcalc.input = { betafile_c1; betafile_c2 };
|
|
matlabbatch{2}.spm.util.imcalc.output = sprintf('con_%04d_minus_%04d.nii', col2, col1);
|
|
matlabbatch{2}.spm.util.imcalc.outdir = { target_folder_final};
|
|
matlabbatch{2}.spm.util.imcalc.expression = 'i2 - i1';
|
|
matlabbatch{2}.spm.util.imcalc.var = struct('name', {}, 'value', {});
|
|
matlabbatch{2}.spm.util.imcalc.options.dmtx = 0;
|
|
matlabbatch{2}.spm.util.imcalc.options.mask = 0;
|
|
matlabbatch{2}.spm.util.imcalc.options.interp = 1;
|
|
matlabbatch{2}.spm.util.imcalc.options.dtype = 4;
|
|
|
|
spm_jobman('run', matlabbatch);
|
|
end
|
|
end
|
|
display(sprintf(' *) Results for roi %s are stored in : %s', roiname, job.target_folder_final));
|
|
end
|
|
display(' *) These are difference in beta images and can be viewed via mricron or xjview');
|
|
|
|
end |