Commit 0275c38c authored by Julian Kosciessa's avatar Julian Kosciessa
Browse files

update analysis scripts

parent 42b63b97
function SpecifyIndividualRegressorMats_4tardis_woutmot_PM
% Create SPM experimental condition info for each subject for later
% batching with tardis
% v3: create Dim-Att regressors, include cue onset, probe onset regressors
% use FAST modelinstead of AR(1): Corbin, N., Todd, N., Friston, K. J., & Callaghan, M. F. (2018). Accurate modeling of temporal correlations in rapidly sampled fMRI time series. Human Brain Mapping, 26(3), 839?14. http://doi.org/10.1002/hbm.24218
% - AP to change back
% v4: AP remove -1 from design timing, but add durations to cue (1 volume) and
% probe (2).
% N = 44 YA + 53 OA;
% AP remove 2131, 2237, 1215 (see oNe Note)
% JQK: set up 2 regressor matrices for F-tests: load + stimulus target feature
IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
'1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';...
'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281';...
'2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';'2130';...
'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';'2149';'2157';...
'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';'2211';'2213';'2214';...
'2215';'2216';'2217';'2219';'2222';'2224';'2226';'2227';'2236';'2238';...
'2241';'2244';'2246';'2248';'2250';'2251';'2252';'2258';'2261'};
%Path to Base Dir of local output
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
mkdir([pn.root, 'B_data/D_batchFiles1stLevelGLM_v2/']);
%addpath([pn.root, 'T_tools/spm12/']); % add spm functions
for indID = 1:numel(IDs)
disp(['Processing ', IDs{indID}]);
matlabbatch = cell(1);
% specify general parameters
matlabbatch{1}.spm.stats.fmri_spec.dir = {[pn.root, 'B_data/F_GLM1stLevelout_v1/' IDs{indID} '/']};
matlabbatch{1}.spm.stats.fmri_spec.timing.units = 'scans';
matlabbatch{1}.spm.stats.fmri_spec.timing.RT = 0.645;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t = 16;
matlabbatch{1}.spm.stats.fmri_spec.timing.fmri_t0 = 8;
matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 0];
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 = {''};
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)';
% subjects with different run number
if strcmp(IDs{indID}, '2131') || strcmp(IDs{indID}, '2237')
numOfRuns = 2;
else
numOfRuns = 4;
end
for indSession = 1:numOfRuns
% get regressor data
pn.regressors = [pn.root, 'B_data/A_regressors/'];
load([pn.regressors, IDs{indID}, '_Run',num2str(indSession),'_regressors.mat'])
% change regressor matrix, such that column 4 represents trial-by-trial
% indices of dimensionality
blockOnsets=find(Regressors(:,1) == 1);
for indBlock = 1:numel(blockOnsets)
if indBlock == numel(blockOnsets)
Regressors(blockOnsets(indBlock):end,4) = ...
Regressors(blockOnsets(indBlock),4);
else
Regressors(blockOnsets(indBlock):blockOnsets(indBlock+1)-1,4) = ...
Regressors(blockOnsets(indBlock),4);
end
end
% For the scans, we need to specify separate 3D images via comma seperator
% Under the hood: design-specified time points are extracted from 2D voxel*time matrix
basefile = ['/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/B3_PLS_mean_v3/B_data/BOLDin/',IDs{indID},'_run-',num2str(indSession),'_feat_detrended_bandpassed_FIX_MNI3mm_std.nii'];
allFiles={};
if strcmp(IDs{indID}, '2132') && indSession == 2
disp(['Change #volumes for run2']);
for indVol = 1:1040
allFiles{indVol,1} = [basefile, ',', num2str(indVol)];
end
else
for indVol = 1:1054
allFiles{indVol,1} = [basefile, ',', num2str(indVol)];
end
end
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).scans = allFiles;
PMDimDur=[];
PMDimOnsets=[];
PM=[];
% stimulus viewing condition
indCond = 1;
for indDim = 1:4
DimOnsets = find(Regressors(:,3) == 1 & Regressors(:,4) == indDim);
PMDimDur=cat(1,PMDimDur,repmat(4, numel(DimOnsets), 1));
PMDimOnsets=cat(1,PMDimOnsets,DimOnsets);
PM=cat(1,PM,repmat(indDim, numel(DimOnsets), 1));
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).name = 'StimOnset';
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).onset = PMDimOnsets;
%onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(indCond).onset;
%matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(indCond).duration = repmat(1, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).duration = PMDimDur;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).tmod = 0;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).pmod.name = 'Stim Load';
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).pmod.param = PM;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).pmod.poly = 1;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(1).orth = 1;
%indCond = indCond + 1;
end
% add cue regressor
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).name = 'CueOnset';
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).onset = find(Regressors(:,2) == 1); % IMPORTANT: SPM starts counting at 0
onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).onset;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).duration = repmat(3, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).tmod = 0;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).pmod = struct('name', {}, 'param', {}, 'poly', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(2).orth = 0;
indCond = indCond + 1;
% add probe regressor
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).name = 'ProbeOnset';
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).onset = find(Regressors(:,11) == 1); % IMPORTANT: SPM starts counting at 0
onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).onset;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).duration = repmat(2, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).tmod = 0;
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).pmod = struct('name', {}, 'param', {}, 'poly', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).cond(3).orth = 0;
%% add regressors
% MotConfoundFile=[pn.root, 'B_data/MotionParameters/' IDs{indID} '/' 'sess-' int2str(indSession) '/' IDs{indID} '_motionout_scol.txt'];
% matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).multi_reg={MotConfoundFile
% [pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indSession) '_motion_6dof.txt']};
% matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).multi_reg={[pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indSession) '_motion_6dof.txt']};
%% add general session information
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).multi = {''};
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).regress = struct('name', {}, 'val', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(indSession).hpf = 128;
matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
matlabbatch{1}.spm.stats.fmri_spec.bases.hrf.derivs = [1 0];
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 = {[pn.root, 'B_data/E_standards/GM_MNI3mm_mask.nii']};
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)';
end % session loop (i.e. runs)
%% Estimate Model
matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model 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;
% %% Now contrasts: Simple Contrast Effects (for now)
% matlabbatch{3}.spm.stats.con.spmmat(1) = cfg_dep('fMRI Contrast Manager: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
%
% % With only one HRF derivative
%
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.name = 'Stimulus Condition';
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.weights = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.sessrep = 'replsc';
%
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.name = 'Load PM 1';
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.weights = [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.sessrep = 'replsc';
%
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.name = 'Load PM 2';
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.weights = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.sessrep = 'replsc';
% %
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.name = 'Dim 3';
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.weights = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.sessrep = 'replsc';
% %
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.name = 'Dim 4';
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.weights = [0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.sessrep = 'replsc';
%
%
% % matlabbatch{4}.spm.stats.results.spmmat = cfg_dep('fMRI Results Report: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
% % matlabbatch{4}.spm.stats.results.conspec.titlestr = 'Dim 1 < Dim2-4';
% % matlabbatch{4}.spm.stats.results.conspec.contrasts = 5;
% % matlabbatch{4}.spm.stats.results.conspec.threshdesc = 'FWE';
% % matlabbatch{4}.spm.stats.results.conspec.thresh = 0.05;
% % matlabbatch{4}.spm.stats.results.conspec.extent = 20;
% % matlabbatch{4}.spm.stats.results.conspec.conjunction = 1;
% % matlabbatch{4}.spm.stats.results.conspec.mask.none = 1;
% % matlabbatch{4}.spm.stats.results.units = 1;
% % matlabbatch{4}.spm.stats.results.export{1}.pdf = true;
%
save([pn.root, 'B_data/D_batchFiles1stLevelGLM_v2/',IDs{indID},'_SPM1stBatchGLM.mat'], 'matlabbatch');
end % subject loop
end % finish completely
\ No newline at end of file
% Run SPM GLM analysis (SPM 6906)
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
addpath([pn.root, 'T_tools/spm12/']);
spm fmri
w = warning ('on','all');
\ No newline at end of file
% Copy DVARS into GLM folder for inclusion as SPM regressor
IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
'1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';...
'1215';'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281';...
'2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';'2130';...
'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';'2149';'2157';...
'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';'2211';'2213';'2214';...
'2215';'2216';'2217';'2219';'2222';'2224';'2226';'2227';'2236';'2238';...
'2241';'2244';'2246';'2248';'2250';'2251';'2252';'2258';'2261'};
pn.DVARSdata = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/preproc/B_data/D_preproc/';
pn.GLMpath = '/Users/kosciessa/Desktop/beegfs/StateSwitch/WIP/G_GLM/B_data/M_MotionParameters/';
for indID = 1:numel(IDs)
for indRun = 1:4
try
curpath = [pn.DVARSdata, IDs{indID}, '/preproc2/run-', num2str(indRun),'/motionout_post/',IDs{indID},'_run-', num2str(indRun), '_feat_detrended_highpassed_denoised_dvars.txt'];
target = [pn.GLMpath,IDs{indID}, '/',IDs{indID},'_run-', num2str(indRun), '_feat_detrended_highpassed_denoised_dvars.txt'];
copyfile(curpath, target);
catch
display('Error');
end
end
end
\ No newline at end of file
function SpecifyIndividualRegressorMats_noMotion_withFeature
function A_SpecifyIndividualRegressorMats_withFeature
% Create SPM experimental condition info for each subject for later
% batching with tardis
......@@ -9,31 +9,20 @@ function SpecifyIndividualRegressorMats_noMotion_withFeature
% v4: AP remove -1 from design timing, but add durations to cue (1 volume) and
% probe (2).
% N = 44 YA + 53 OA;
% AP remove 2131, 2237, 1215 (see oNe Note)
% JQK: set up 2 regressor matrices for F-tests: load + stimulus target feature
IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
'1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';...
'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281';...
'2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';'2130';...
'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';'2149';'2157';...
'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';'2211';'2213';'2214';...
'2215';'2216';'2217';'2219';'2222';'2224';'2226';'2227';'2236';'2238';...
'2241';'2244';'2246';'2248';'2250';'2251';'2252';'2258';'2261'};
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
addpath([pn.root, 'T_tools/spm12/']);
spm fmri
'1215';'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281'};
%pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
pn.root = '/home/beegfs/kosciessa/StateSwitch/WIP/G_GLM/';
%Path to Base Dir of local output
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
pn.batchPath = [pn.root, 'B_data/D_batchFiles1stLevelGLM_v3/'];
mkdir(pn.batchPath);
%addpath([pn.root, 'T_tools/spm12/']); % add spm functions
pn.filePathRoot = '/Users/kosciessa/Desktop/beegfs/StateSwitch/WIP/G_GLM/';
pn.batchPath = [pn.filePathRoot, 'B_data/D_batch1stLevel_v2/'];
addpath([pn.filePathRoot, 'T_tools/spm12/']); % add spm functions
spm fmri
for indID = 1:numel(IDs)
......@@ -54,8 +43,7 @@ for indID = 1:numel(IDs)
matlabbatch{1}.spm.stats.fmri_spec.global = 'None';
matlabbatch{1}.spm.stats.fmri_spec.mthresh = -Inf;
matlabbatch{1}.spm.stats.fmri_spec.mask = {''};
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)';
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'FAST';
% subjects with different run number
......@@ -68,7 +56,7 @@ for indID = 1:numel(IDs)
for indRun = 1:numOfRuns
% get regressor data
pn.regressors = [pn.root, 'B_data/A_regressors/'];
pn.regressors = [pn.filePathRoot, 'B_data/A_regressors/'];
load([pn.regressors, IDs{indID}, '_Run',num2str(indRun),'_regressors.mat'])
% change regressor matrix, such that column 4 represents trial-by-trial
......@@ -86,7 +74,7 @@ for indID = 1:numel(IDs)
% For the scans, we need to specify separate 3D images via comma seperator
% Under the hood: design-specified time points are extracted from 2D voxel*time matrix
basefile = ['/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/B3_PLS_mean_v3/B_data/BOLDin/',IDs{indID},'_run-',num2str(indRun),'_feat_detrended_bandpassed_FIX_MNI3mm_std.nii'];
basefile = [pn.root, 'B_data/BOLDin/',IDs{indID},'_run-',num2str(indRun),'_regressed.nii'];
allFiles={};
if strcmp(IDs{indID}, '2132') && indRun == 2
......@@ -105,16 +93,16 @@ for indID = 1:numel(IDs)
% stimulus viewing condition
% To get relevant-feature-specific conditions, we have to create 4
% stim conditions: color, direction, size & luminance
AttributeNames = {'color'; 'direction'; 'size'; 'luminance'};
AttributeNames = {'1Target'; '2Targets'; '3Targets'; '4Targets'};
for indAtt = 1:4 % here: contained in target set
AttOnsets = find(Regressors(:,3) == 1 & Regressors(:,6+indAtt) == 1);
PM = Regressors(AttOnsets,4);
AttOnsets = find(sum(Regressors(:,7:10),2) == indAtt);
PM = Regressors(AttOnsets,6); % add sequence position as parametric regressor
AttDur = repmat(4, numel(AttOnsets),1);
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).name = AttributeNames{indAtt};
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).onset = AttOnsets;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).duration = AttDur;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).tmod = 0;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.name = 'Stim Load';
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.name = 'SequencePosition';
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.param = PM;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).pmod.poly = 1;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(indAtt).orth = 1;
......@@ -129,25 +117,24 @@ for indID = 1:numel(IDs)
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).pmod = struct('name', {}, 'param', {}, 'poly', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(5).orth = 0;
% add probe regressor
% add probe regressor with parametric RT modulation
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).name = 'ProbeOnset';
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).onset = find(Regressors(:,11) == 1); % IMPORTANT: SPM starts counting at 0
onsets = matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).onset;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).duration = repmat(2, numel(onsets), 1); clear onsets; % duration of 1 (VarToolbox) vs 0 (SPM convention)
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).tmod = 0;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.name = 'RT';
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.param = Regressors((Regressors(:,11) == 1),13);
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod.poly = 1;
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).pmod = struct('name', {}, 'param', {}, 'poly', {});
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).cond(6).orth = 0;
%% add 24 Friston motion regressors, and DVARS
%% add regressors
% MotConfoundFile=[pn.root, 'B_data/MotionParameters/' IDs{indID} '/' 'sess-' int2str(indRun) '/' IDs{indID} '_motionout_scol.txt'];
% matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi_reg={MotConfoundFile
% [pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indRun) '_motion_6dof.txt']};
% matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi_reg={[pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indRun) '_motion_6dof.txt']};
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi_reg={...
[pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_sess-' int2str(indRun) '_motion_24dof.txt'];...
[pn.root, 'B_data/M_MotionParameters/' IDs{indID} '/' IDs{indID} '_run-' int2str(indRun) '_feat_detrended_highpassed_denoised_dvars.txt']};
%% add general session information
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).multi = {''};
matlabbatch{1}.spm.stats.fmri_spec.sess(indRun).regress = struct('name', {}, 'val', {});
......@@ -157,8 +144,8 @@ for indID = 1:numel(IDs)
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 = {[pn.root, 'B_data/E_standards/GM_MNI3mm_mask.nii']};
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'AR(1)';
matlabbatch{1}.spm.stats.fmri_spec.mask = {[pn.root, 'B_data/E_standards/mni_icbm152_gm_tal_nlin_sym_09c_MNI_3mm.nii']};
matlabbatch{1}.spm.stats.fmri_spec.cvi = 'FAST';
end % session loop (i.e. runs)
%% Estimate Model
......@@ -166,44 +153,7 @@ for indID = 1:numel(IDs)
matlabbatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model 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;
% %% Now contrasts: Simple Contrast Effects (for now)
% matlabbatch{3}.spm.stats.con.spmmat(1) = cfg_dep('fMRI Contrast Manager: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
%
% % With only one HRF derivative
%
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.name = 'Stimulus Condition';
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.weights = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% matlabbatch{3}.spm.stats.con.consess{1}.tcon.sessrep = 'replsc';
%
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.name = 'Load PM 1';
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.weights = [0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% matlabbatch{3}.spm.stats.con.consess{2}.tcon.sessrep = 'replsc';
%
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.name = 'Load PM 2';
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.weights = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.sessrep = 'replsc';
% %
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.name = 'Dim 3';
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.weights = [0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{3}.tcon.sessrep = 'replsc';
% %
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.name = 'Dim 4';
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.weights = [0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0];
% % matlabbatch{3}.spm.stats.con.consess{4}.tcon.sessrep = 'replsc';
%
%
% % matlabbatch{4}.spm.stats.results.spmmat = cfg_dep('fMRI Results Report: SPM.mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
% % matlabbatch{4}.spm.stats.results.conspec.titlestr = 'Dim 1 < Dim2-4';
% % matlabbatch{4}.spm.stats.results.conspec.contrasts = 5;
% % matlabbatch{4}.spm.stats.results.conspec.threshdesc = 'FWE';
% % matlabbatch{4}.spm.stats.results.conspec.thresh = 0.05;
% % matlabbatch{4}.spm.stats.results.conspec.extent = 20;
% % matlabbatch{4}.spm.stats.results.conspec.conjunction = 1;
% % matlabbatch{4}.spm.stats.results.conspec.mask.none = 1;
% % matlabbatch{4}.spm.stats.results.units = 1;
% % matlabbatch{4}.spm.stats.results.export{1}.pdf = true;
%
save([pn.batchPath, IDs{indID},'_SPM1stBatchGLM.mat'], 'matlabbatch');
end % subject loop
......
#!/bin/bash
RootPath=/home/beegfs/kosciessa/StateSwitch/WIP/G_GLM
DATADIR=${RootPath}/B_data/D_batch1stLevel_v2
#subjList="1117 1118 1120 1124 1125 1126 1131 1132 1135 1136 1151 1160 1164 1167 1169 1172 1173 1178 1182 1214 1215 1216 1219 1223 1227 1228 1233 1234 1237 1239 1240 1243 1245 1247 1250 1252 1257 1261 1265 1266 1268 1270 1276 1281"
#subjList="1117"
subjList="1215"
for subj in $subjList; do
echo "#PBS -N STSWD_SPM_${subj}" > job
echo "#PBS -l nodes=1:ppn=1" >> job
echo "#PBS -l mem=16gb" >> job
echo "#PBS -l walltime=30:00:00" >> job
echo "#PBS -j oe" >> job
echo "#PBS -o ${RootPath}/Y_logs/" >> job
echo "#PBS -m n" >> job
echo "#PBS -d ." >> job
echo "module load spm12 ; run_spm12.sh /mcr run ${DATADIR}/${subj}_SPM1stBatchGLM.mat" >> job
qsub job
rm job
done
% Specify design model at first level
% use concatenation across runs and correct for concatenation in design matrix
% N = 44 YA + 53 OA;
% N = 44 YA
IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
'1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';'1215';...
'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281';...
'2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';'2130';...
'2131';'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';'2149';'2157';...
'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';'2211';'2213';'2214';...
'2215';'2216';'2217';'2219';'2222';'2224';'2226';'2227';'2236';'2237';'2238';...
'2241';'2244';'2246';'2248';'2250';'2251';'2252';'2258';'2261'};
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281'};
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
pn.root = '/Users/kosciessa/Desktop/beegfs/StateSwitch/WIP/G_GLM/';
addpath([pn.root, 'T_tools/spm12/']); % add spm functions
for indID = 1:numel(IDs)
......@@ -20,11 +15,11 @@ for indID = 1:numel(IDs)
spm('defaults','fmri')
spm_jobman('initcfg');
% specify SPM job
load([pn.root, 'B_data/D_batchFiles1stLevelGLM/',IDs{indID},'_SPM1stBatchGLM_cat.mat'], 'matlabbatch', 'N');
% run job
spm_jobman('run', matlabbatch);
load([pn.root, 'B_data/D_batch1stLevel_v2/',IDs{indID},'_SPM1stBatchGLM.mat'], 'matlabbatch');
% % run job
% spm_jobman('run', matlabbatch);
% correct matrix for independent runs
% https://en.wikibooks.org/wiki/SPM/Concatenation
spm_fmri_concatenate([matlabbatch{1}.spm.stats.fmri_spec.dir{1},'SPM.mat'], N);
spm_fmri_concatenate([matlabbatch{1}.spm.stats.fmri_spec.dir{1},IDs{indID},'_SPM1stBatchGLM.mat'], [1054, 1054, 1054, 1054]);
clear matlabbatch N;
end
% List of open inputs
nrun = X; % enter the number of runs here
jobfile = {'/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/A_scripts/B_runVarTBX_job.m'};
jobs = repmat(jobfile, 1, nrun);
inputs = cell(0, nrun);
for crun = 1:nrun
end
spm('defaults', 'FMRI');
spm_jobman('run', jobs, inputs{:});
%-----------------------------------------------------------------------
% Job saved on 07-Aug-2018 12:36:45 by cfg_util (rev $Rev: 6460 $)
% spm SPM - SPM12 (6906)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
% N = 44 YA + 53 OA;
IDs = {'1117';'1118';'1120';'1124';'1125';'1126';'1131';'1132';'1135';'1136';...
'1151';'1160';'1164';'1167';'1169';'1172';'1173';'1178';'1182';'1214';'1215';...
'1216';'1219';'1223';'1227';'1228';'1233';'1234';'1237';'1239';'1240';'1243';...
'1245';'1247';'1250';'1252';'1257';'1261';'1265';'1266';'1268';'1270';'1276';'1281';...
'2104';'2107';'2108';'2112';'2118';'2120';'2121';'2123';'2125';'2129';'2130';...
'2131';'2132';'2133';'2134';'2135';'2139';'2140';'2145';'2147';'2149';'2157';...
'2160';'2201';'2202';'2203';'2205';'2206';'2209';'2210';'2211';'2213';'2214';...
'2215';'2216';'2217';'2219';'2222';'2224';'2226';'2227';'2236';'2237';'2238';...
'2241';'2244';'2246';'2248';'2250';'2251';'2252';'2258';'2261'};
pn.root = '/Volumes/LNDG/Projects/StateSwitch/dynamic/data/mri/task/analyses/G_GLM/';
addpath([pn.root, 'T_tools/spm12/']); % add spm functions
for indID = 1%:numel(IDs)
% initialize SPM batch processing
spm('defaults','fmri')
spm_jobman('initcfg');
% specify SPM job
matlabbatch{1}.spm.tools.variability.modeltype = 'lss';
matlabbatch{1}.spm.tools.variability.modelmat = {[pn.root, 'B_data/D_batchFiles1stLevelGLM_v3/',IDs{indID},'_SPM1stBatchGLM.mat']};
matlabbatch{1}.spm.tools.variability.metric = 'sd';
matlabbatch{1}.spm.tools.variability.resultprefix = ['sd_', IDs{indID}];
matlabbatch{1}.spm.tools.variability.resultdir = {[pn.root, 'B_data/F_GLM1stLevelout_v1/A_SPMfiles']};
% run job
spm_jobman('run', matlabbatch);