MRI学习笔记-合并多个session的fMRI数据

原因：

实验共有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');