原因:
实验共有12个session,在某些session中有些condition不存在,在一阶分析和二阶分析时会报错,因此在预处理后将12个session的数据合并为1个进行后续一阶分析和二阶分析。
文件名设定:
1.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri路径下包括所有处理步骤后生成的数据
2.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri\data_04_analysis2路径下为预处理后每个被试的数据
3.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri\data_04_analysis2\01路径下分为各个session的数据
4.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri\data_05_1st_level路径下建立模型1(Model1,M1)的文件夹,用于放置分成12sessions进行一阶分析的结果
5.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri\data_05_1st_level\M1路径下为每个被试(01,02,03...)的一阶分析结果
6.D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri\data_05_1st_level2为合并数据后放置一阶分析结果的路径,层级设置同 data_05_1st_level
第一步:
先按照分成12sessions进行一阶分析,结果保存在 data_05_1st_level 中,用于后续读取 SPM.mat 中的 SPM.xX.K(1,nsess).X0。
第二步:
这一步是为了将12个sessions中的SPM.xX.K(1,nsess).X0,头动参数,runs(不知道这个是啥,可能是代表哪个session的意思,因为在对应session中填充数字1)拼接在一起,在后续一阶分析时使用。
修改以下脚本中注明需修改的部分,修改好后运行。
clc;
clear;
projectdir = 'D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri';%最上层文件路径,需修改
odatadir = 'data_04_analysis2';%original data folder,需修改
fdatadir = 'data_05_1st_level';%1st data folder,需修改
numsess = 12;%session数量,需修改
%session number,每个被试data_04_analysis2中session文件命名,这里为session01/session02.../session12,可以根据自己的文件命名修改
for nsess = 1:numsess
if nsess < 10
sessions{nsess} = strcat('session0',num2str(nsess));
else
sessions{nsess} = strcat('session',num2str(nsess));
end
end
%subject number,被试编号为1:27(需修改),被试的编号为01/02/03.../27,可以根据自己的被试编号命名修改
for nsub = 1:27
if nsub < 10
subjects{nsub} = strcat('0',num2str(nsub));
else
subjects{nsub} = num2str(nsub);
end
%读取data_05_1st_level一阶分析的SPM.mat中的SPM.xX.K(1,nsess).X0参数
trendfile = fullfile(projectdir,fdatadir, 'M1', subjects{nsub}, 'SPM.mat');
SPM = importdata(trendfile);% loads SPM object
for nsess = 1:numsess
trendsa{nsess} = SPM.xX.K(1,nsess).X0;
end
trends = blkdiag(trendsa{1},trendsa{2},trendsa{3},trendsa{4},trendsa{5},trendsa{6},trendsa{7},trendsa{8},trendsa{9},...
trendsa{10},trendsa{11},trendsa{12});%有多少session就按照顺序写多少个trendsa{n},需修改
runs = blkdiag(ones(size(trendsa{1},1),1),ones(size(trendsa{2},1),1),ones(size(trendsa{3},1),1),ones(size(trendsa{4},1),1)...
,ones(size(trendsa{5},1),1),ones(size(trendsa{6},1),1),ones(size(trendsa{7},1),1),ones(size(trendsa{8},1),1)...
,ones(size(trendsa{9},1),1),ones(size(trendsa{10},1),1),ones(size(trendsa{11},1),1),ones(size(trendsa{12},1),1));%有多少session就按照顺序写多少个ones(size(trendsa{n},1),1),需修改
%读取head motion
for nsess = 1:numsess
rpfname = fullfile(projectdir,odatadir, subjects{nsub}, sessions{nsess}, '*.txt');
rpfname = dir(rpfname);
rpfname = rpfname.name;
rpdata = importdata(fullfile(projectdir,odatadir, subjects{nsub}, sessions{nsess}, rpfname));
rpdatac{nsess} = rpdata - repmat(mean(rpdata), size(rpdata,1),1);
end
movements = blkdiag(rpdatac{1},rpdatac{2},rpdatac{3},rpdatac{4},rpdatac{5},rpdatac{6},rpdatac{7},rpdatac{8},rpdatac{9},...
rpdatac{10},rpdatac{11},rpdatac{12});%有多少session就按照顺序写多少个rpdatac{n},需修改
regressors{nsub} = [trends movements runs];%将刚才生成的trends movements runs拼接起来
spy(regressors{nsub})%最后会生成一个图,也没啥用,看看就行
% save the trend, motion and run regressors
dlmwrite(fullfile(projectdir,odatadir, subjects{nsub}, ['trend_motion_' subjects{nsub} '.txt']), regressors{nsub}, 'delimiter', '\t');%会在data_04_analysis2每个被试的文件夹下生成一个trend_motion_01.txt的文件,01代表被试编号
end
第三步:
在进行新的一阶分析前需要修改onset,因为将多个sessions拼接后,图像在原来时间轴上的位置变化了。
第一个session的onset无需修改,从第二个session开始修改。在原来onset的基础上增加 前面所包含图像的个数*TR。比如session01的图像数量为200,session02的图像数量为220,TR = 2 s,则session02的onset应该在原有基础上增加 200*2 = 200 s,则session03的onset应该在原有基础上增加 (200+220)*2 = 840 s。
第四步:
重新生成conditions文件,见之前一阶分析教程
第五步:
重新进行一阶分析,与原有一阶分析脚本有稍许不同
修改好后运行脚本,将在jobfiles文件夹中生成的single_sess_01(01为被试编号)全选进SPM的batch中运行完成一阶分析。
% Create jobfiles for fixed effects analysis
clc;
clear;
spm_defaults
%% FFX variables
conditions = {'CertainGoAll','CertainStopLeft','CertainStopRight','UncertainGo',...
'UncertainStopLeftSuccess','UncertainStopLeftFailed','UncertainStopRightSuccess','UncertainStopRightFailed',...
'UncertainStopAllSuccess','UncertainStopAllFailed','errors'};%可更改,conditions名称,注意顺序,后续顺序都要一致
% Script options
est = 1; % Include model estimation job? 1=yes,可修改
% conman = 0; % Include contrast manager job?,可修改
conman = 1;
% Timing parameters
units = 'secs'; % 'scans' or 'secs'
TR = 2; % TR in seconds,可修改
tbins = 62; % Microtime bins,可修改
tbinref = 32; % Reference time bin,可修改
nmodel = 1;%可修改
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check correct version of SPM defaults loaded before running jobs %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Directory names
rootdir = 'D:\LLYdata\motor_inhibition_data2\motor_inhibition_fmri';%可修改
imagedir = 'data_04_analysis2';%可修改
ffxdir = strcat('data_05_1st_level2\M',num2str(nmodel));% Directory to save analysis,可修改 (created under subject directory)
rfxdir = strcat('data_06_2st_level\M',num2str(nmodel));%可更改,rfx所在文件夹
rpexp = 'trend_*.txt';%head motion
ffxname = 'single_sess';
%
for sub = 1:27%被试编号,可修改
if sub < 10
subdir{sub} = strcat('0',num2str(sub));
else
subdir{sub} = num2str(sub);
end
end
nsub = length(subdir);
sessdir = {'session01'}; % Directories for functional images (ones per session)
nsess = length(sessdir);
% Functional image file prefixes (row=subject, col=session)
functprefix = cell(nsub,nsess);
for sub = 1:nsub
for sess = 1:nsess
functprefix{sub,sess} = 'swa';%一阶分析所用.nii文件的前缀,可修改
end
end
%%% Only functional images to be included in the model should match functprefix in sessdir
%%% Move any others to 'dummy_scans'
% Multiple conditions files (row=subject, col=session)
multicond = cell(nsub,nsess);
for sub = 1:nsub
% for sess = 1:nsess
multicond{sub,sess} = fullfile(rootdir,ffxdir,'conditions_files',['conditions_',subdir{sub},'.mat']);
% multicond{sub,sess} = fullfile(rootdir,ffxdir,'conditions_files',[subdir{sub},'_session',num2str(sess),'_conditions.mat']);
% end
end
% Explicit mask files (row = subject)
exmask = cell(nsub,1);
for sub = 1:nsub
exmask{sub,1} = ''; % Explicit mask filename
end
% Temporal / dispersion derivatives
derivs = [0 0]; % None: [0 0], time derivative: [1 0], time & dispersion: [1 1]
% High pass filter
hpf = Inf; % Value in seconds, or disable using 'Inf' without quotes
%% Contrasts
delcon=1; % Delete existing contrasts
% Work out contrast vector for each condition
sesscols = length(conditions)+6;
cols = sesscols*nsess+nsess;
sessind = 1:nsess;
colind = (sessind-1)*sesscols+1;
for cond = 1:length(conditions)
thiscolind = colind+cond-1;
eval(sprintf('%s=zeros(1,%d);',conditions{cond},cols))
eval(sprintf('%s(thiscolind)=1;',conditions{cond}))
end
null=zeros(1,cols);
% T-contrasts
consess = [];
conpos = [];
conneg = [];
%各条件与null互相比较
% for cond = 1:length(conditions)
for cond = 1:length(conditions)
consess{1,end+1}.tcon.name = eval(['strcat(conditions{1,',num2str(cond),'} ,''-null'')']);%1:2:2*length(conditions)
consess{1,end+1}.tcon.name = eval(['strcat(''null-'',conditions{1,',num2str(cond),'})']);%2:2:2*length(conditions)
conpos{end+1} = conditions{1,cond};% 生成rfx需要的Positive condition
conneg{end+1} = 'null';% 生成rfx需要的Negative condition
end
% 条件之间比较
contrast = [1,2;1,3;1,4;2,3;4,5;4,6;4,7;4,8;4,9;4,10;5,6;5,7;5,9;6,8;6,10;7,8;7,9;8,10;9,10;2,5;3,7];%可更改,所需要比较的条件在conditions中的列数
for ncontrast = 1:length(contrast)
a = contrast(ncontrast,1);
b = contrast(ncontrast,2);
consess{1,end+1}.tcon.name = strcat(conditions{1,a},'-',conditions{1,b});
consess{1,end+1}.tcon.name = strcat(conditions{1,b},'-',conditions{1,a});
conpos{end+1} = conditions{1,a};% 生成rfx需要的Positive condition
conneg{end+1} = conditions{1,b};% 生成rfx需要的Negative condition
end
%生成的Positive condition和Negative condition保存到rfx的文件夹中
save(fullfile(rootdir,rfxdir,'conpos.mat'),'conpos');
save(fullfile(rootdir,rfxdir,'conneg.mat'),'conneg');
% Evaluate contrasts from names
confile=fopen(fullfile(rootdir,ffxdir,'jobfiles','contrast_numbers.txt'),'w');
for c = 1:size(consess,2)
consess{1,c}.tcon.convec = eval(consess{1,c}.tcon.name);
fprintf(confile,'Contrast %02d: %s\n',c,consess{1,c}.tcon.name);
end
fclose(confile);
%% Create directories, find files
ffxpath = cell(1,nsub);
scanfiles = cell(nsub,nsess);
multireg = cell(nsub,nsess);
if ~exist(fullfile(rootdir,ffxdir,'jobfiles'),'dir')
mkdir(fullfile(rootdir,ffxdir,'jobfiles'));
end
for sub = 1:nsub
ffxpath{1,sub} = fullfile(rootdir,ffxdir,subdir{sub});
% Create ffx directory if needed
if ~exist(ffxpath{1,sub},'dir')
mkdir(ffxpath{1,sub});
end
% Scan filenames
blockdir = {'session01','session02','session03','session04','session05','session06','session07','session08',...
'session09','session10','session11','session12'}; % Directories for functional images (ones per session),有多少session就写多少,可修改
nblk = length(blockdir);
for blk = 1:nblk
raffilesa{blk} = dir(fullfile(rootdir,imagedir,subdir{sub},blockdir{blk},'swa*.nii'));%swa可修改
raffilesacell{blk} = cell(size(raffilesa{blk}));
[raffilesacell{blk}{:}] = deal(raffilesa{blk}.name);
for h = 1:size(raffilesa{blk},1)
raffilesacell{blk}{h} = fullfile(rootdir,imagedir,subdir{sub},blockdir{blk},raffilesa{blk}(h).name);
end
end
functfiles = [raffilesacell{1,1,:};raffilesacell{1,2,:};raffilesacell{1,3,:};raffilesacell{1,4,:};raffilesacell{1,5,:};raffilesacell{1,6,:};raffilesacell{1,7,:};raffilesacell{1,8,:};...
raffilesacell{1,9,:};raffilesacell{1,10,:};raffilesacell{1,11,:};raffilesacell{1,12,:}];%有多少session就写多少raffilesacell{1,n,:}
scanfiles{sub,sess} = functfiles;
% Multiple regressor filenames
rpfile = dir(fullfile(rootdir,imagedir,subdir{sub},rpexp));
if size(rpfile,1) == 0
error('Multiple regressor file not found using %s', fullfile(rootdir,imagedir,subdir{sub},rpexp))
elseif size(rpfile,1) > 1
error('More than one multiple regressor file found using %s',fullfile(rootdir,imagedir,subdir{sub},rpexp))
else
multireg{sub,sess} = fullfile(rootdir,imagedir,subdir{sub},rpfile.name);
end
end
%% Create jobfiles
jobs = cell(1,1);
jobfile = cell(nsub,1);
for sub = 1:nsub
% Specification
jobs{1,1}.stats{1}.fmri_spec.dir = {ffxpath{1,sub}};
jobs{1,1}.stats{1}.fmri_spec.timing.units = units;
jobs{1,1}.stats{1}.fmri_spec.timing.RT = TR;
jobs{1,1}.stats{1}.fmri_spec.timing.fmri_t = tbins;
jobs{1,1}.stats{1}.fmri_spec.timing.fmri_t0 = tbinref;
for sess=1:nsess
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).scans = scanfiles{sub,sess};
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).cond = struct([]);
% matlabbatch{1}.spm.stats.fmri_spec.sess(sess).cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {});
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).multi = {multicond{sub,sess}};
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).regress = struct([]);
%matlabbatch{1}.spm.stats.fmri_spec.sess(sess).regress = struct('name', {}, 'val', {});
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).multi_reg = {multireg{sub,sess}};
jobs{1,1}.stats{1}.fmri_spec.sess(1,sess).hpf = hpf;
end
jobs{1,1}.stats{1}.fmri_spec.fact = struct([]);
%matlabbatch{1}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
jobs{1,1}.stats{1}.fmri_spec.bases.hrf.derivs = derivs;
jobs{1,1}.stats{1}.fmri_spec.volt = 1;
jobs{1,1}.stats{1}.fmri_spec.global = 'None';
jobs{1,1}.stats{1}.fmri_spec.mask = {exmask{sub,1}};
jobs{1,1}.stats{1}.fmri_spec.cvi = 'AR(1)';
% Estimation
if est==1
jobs{1,1}.stats{1,2}.fmri_est.spmmat = {fullfile(ffxpath{1,sub},'SPM.mat')};
jobs{1,1}.stats{1,2}.fmri_est.method.Classical = 1;
end
% Contrast manager
if conman==1
jobs{1,1}.stats{1,3}.con.spmmat = {fullfile(ffxpath{1,sub},'SPM.mat')};
jobs{1,1}.stats{1,3}.con.consess = consess;
jobs{1,1}.stats{1,3}.con.delete = delcon;
end
% Save jobfile
jobfilename = fullfile(rootdir,ffxdir,'jobfiles',[ffxname,'_',subdir{sub},'.mat']);
jobfile{sub,1} = jobfilename;
save(jobfile{sub,1}, 'jobs')
end
spm_defaults;
spm('chmod','fmri');