引言
先前简单介绍了PPI分析的基本原理,提到了ROI选取的问题,本文就ROI选取以及mask制作的流程进行详细阐述。为便于以后查阅,也将之前的PPI分析批处理code放在一起。
一、ROI选取
在SPM12的示例数据中,针对单个被试的PPI处理,选取ROI“V2”时,是采用的脑激活图与V2的peak点坐标为核心的半径为6mm的球的交集。
%示例数据VOI提取code
matlabbatch{5}.spm.util.voi.name = 'V2';
matlabbatch{5}.spm.util.voi.roi{1}.spm.spmmat = {''};
matlabbatch{5}.spm.util.voi.roi{1}.spm.contrast = 3;
matlabbatch{5}.spm.util.voi.roi{1}.spm.threshdesc = 'FWE';
matlabbatch{5}.spm.util.voi.roi{1}.spm.thresh = 0.05;
matlabbatch{5}.spm.util.voi.roi{1}.spm.extent = 0;
matlabbatch{5}.spm.util.voi.roi{2}.sphere.centre = [15 -78 -9];
matlabbatch{5}.spm.util.voi.roi{2}.sphere.radius = 6;
matlabbatch{5}.spm.util.voi.roi{2}.sphere.move.local.spm = 1;
matlabbatch{5}.spm.util.voi.expression = 'i1 & i2';
不过由于示例数据只是做了单个被试,而我们分析的是两组,具有对照组。所以参考吴士豪的博士论文,我们决定将两组组分析结果的激活图进行重叠,重叠后的区域再和选定的brodmann功能区或者说AAl解剖区域进行重叠,最后得到的mask作为感兴趣区域。
二、ROImask的制作
1.将RW组与SD组正激活、负激活的cluster在DPABI中分别保存为mask,得到文件如下:
2.这里演示两组正激活图像重叠code过程:
%-------------------------------
%时间:20220922
%程序员:wishing
%目的:求两张激活图cluster.nii的交集,
%--------------------------------------------
clear,clc;
cd 'D:\desk\speed\data_1\Dp_spm\PPI\ROImask\mask_make';
%重叠RW与SD两组的激活图
%读取两个nii文件信息
W1=spm_vol('rwred_mask.nii');
W2=spm_vol('sdred_mask.nii');
%读取两个nii文件的体积
Y1=spm_read_vols(W1);
Y2=spm_read_vols(W2);
%展示nii文件信息
%imshow(Y1(:,:,30));
%取交集
Y2=fillmissing(Y2,'constant',0);%将Y2中的nan值替换为0,若不存在则不替换
YY=Y1&Y2;
%命名,保存为mask.nii
V=W1;
V.fname='RSred_mask.nii';
spm_write_vol(V,YY);
3.根据经验选取BA9区域,提取MRIcron中Brodmann模板的BA9区
1).打开DPABI的view界面
2).overlay选择MRIcron中的BA模板
3).设置阈值,设置9就是BA9区,然后保存cluster就能得到BA9的mask
4.改变BA9mask文件的维度,使其和SPM处理结果的文件维度相同、
1)这里需要准备一个SPM处理结果的文件,这里我选取了GLM分析后的beta_0009.nii.
2)使用函数处理。
P1='BA9.nii';P2='beta_0009.nii';
spm_mask(P1,P2);
%这样P1的维度就和P2一样了,而且生成了一个mP2文件,
%注意此函数格式为: FORMAT spm_mask(P1, P2, thresh)
% P1 - matrix of input image filenames from which
% to compute the mask.
%P2 - matrix of input image filenames on which
% to apply the mask.
%thresh - optional threshold(s) for defining the mask.
%The masked images are prepended with the prefix `m'.
%有时需要设置thresh为1
5.再将重叠后的图像与BA9(根据经验选取的模板)重叠,这里与第2步的code类似,注意BA9图像中可能有空值,这时就需要用到Y2=fillmissing(Y2,'constant',0);这个函数。
6.最后得到生成的mask文件。
三、批量生成批处理需要的文件夹
%时间:20220919
%程序员:wishing
%优化后的批量制作文件夹
%--------------------------------------------------------------------------
% 设置循环
clear,clc;
%————————————————————————————————
for i=1:37
p(i) = string(num2str(i, "%02d"));%生成“01”“02”一直到“37”的字符串
end
%--------------------------------------------------------------------------
% 设置文件夹
%————————————————————————————————
for n=1:length(p)
filename=string(append('D:\desk\speed\data_1\Dp_spm\PPI\RSBA18\RW\sub_1',p(n)));
mkdir(filename);
end
四、PPI批处理流程code
%-----------------------------------------------------------------------
% Job saved on 19-Sep-2022 11:16:50 by cfg_util (rev $Rev: 7345 $)
% spm SPM - SPM12 (7771)
% cfg_basicio BasicIO - Unknown
%时间:20220919
%程序员:wishing
%目的:实现任务态脑区激活分析后的PPI全流程循环处理,包含VOI提取,PPI变量制作,PPIGLM模型设计与contrast设计
%-----------------------------------------------------------------------
% Initialise SPM
%--------------------------------------------------------------------------
clear,clc;
spm('Defaults','fMRI');
spm_jobman('initcfg');
% spm_get_defaults('cmdline',true);
%--------------------------------------------------------------------------
% 设置循环
%—————————————————————————————————————
for i=1:37
p(i) = string(num2str(i, "%02d"));%生成“01”“02”一直到“37”的字符串
end
for n=1:length(p)
% 设置文件路径
%--------------------------------------------------------------------------
filename_out=char(append('D:\desk\speed\data_1\Dp_spm\PPI\RSBA18\RW\sub_1',p(n)));
filename_spm=char(append('D:\desk\speed\data_1\Dp_spm\RW\FunImg\sub_1',p(n),'\SPM.mat'));
filename_in=char(append('D:\desk\speed\data_1\Dp_spm\RW\FunImg\sub_1',p(n),'\'));
filename_rp=char(append('D:\desk\speed\data_1\Dp_spm\RW\FunImg\sub_1',p(n),'\rp_1',p(n),'_speed_00001.txt'));
file_ppi=char(append('D:\desk\speed\data_1\Dp_spm\RW\FunImg\sub_1',p(n),'\PPI_s-b_RSBA18.mat'));
file_mask='D:\desk\speed\data_1\Dp_spm\PPI\ROImask\mask_make\mask\RSred_ba18_mask.nii,1';
%--------------------------------------------------------------------------
% VOI提取
%--------------------------------------------------------------------------
matlabbatch{1}.spm.util.voi.spmmat = {filename_spm};
matlabbatch{1}.spm.util.voi.adjust = 0;
matlabbatch{1}.spm.util.voi.session = 1;
matlabbatch{1}.spm.util.voi.name = 'RSBA18';
% matlabbatch{1}.spm.util.voi.roi{1}.spm.spmmat = {''};
% matlabbatch{1}.spm.util.voi.roi{1}.spm.contrast = 1;
% matlabbatch{1}.spm.util.voi.roi{1}.spm.threshdesc = 'FWE';
% matlabbatch{1}.spm.util.voi.roi{1}.spm.thresh = 0.05;
% matlabbatch{1}.spm.util.voi.roi{1}.spm.extent = 0;
matlabbatch{1}.spm.util.voi.roi{1}.mask.image = {file_mask};
matlabbatch{1}.spm.util.voi.roi{1}.mask.threshold = 0.5;
matlabbatch{1}.spm.util.voi.expression = 'i1';
%--------------------------------------------------------------------------
% PPI交互作用变量生成
%--------------------------------------------------------------------------
matlabbatch{2}.spm.stats.ppi.spmmat = {filename_spm};
matlabbatch{2}.spm.stats.ppi.type.ppi.voi(1) = cfg_dep('Volume of Interest: VOI mat File', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','voimat'));
matlabbatch{2}.spm.stats.ppi.type.ppi.u = [1 1 1
2 1 -1];%这里还不确定,第一列是condition数,第二列是con_name,第三列是权重
matlabbatch{2}.spm.stats.ppi.name = 's-b_RSBA18';
matlabbatch{2}.spm.stats.ppi.disp = 0;
%--------------------------------------------------------------------------
% GLM PPImodel设计
%--------------------------------------------------------------------------
matlabbatch{3}.spm.stats.fmri_spec.dir = {filename_out};
matlabbatch{3}.spm.stats.fmri_spec.timing.units = 'secs';
matlabbatch{3}.spm.stats.fmri_spec.timing.RT = 2;
matlabbatch{3}.spm.stats.fmri_spec.timing.fmri_t = 35;
matlabbatch{3}.spm.stats.fmri_spec.timing.fmri_t0 = 18;
%session
f = spm_select('FPList',filename_in, '^swr.*\.nii$');
matlabbatch{3}.spm.stats.fmri_spec.sess.scans = cellstr(f);
%regressors PPI交互项和头动参数
matlabbatch{3}.spm.stats.fmri_spec.sess.cond = struct('name', {}, 'onset', {}, 'duration', {}, 'tmod', {}, 'pmod', {}, 'orth', {});
matlabbatch{3}.spm.stats.fmri_spec.sess.multi = {''};
matlabbatch{3}.spm.stats.fmri_spec.sess.regress = struct('name', {}, 'val', {});
matlabbatch{3}.spm.stats.fmri_spec.sess.multi_reg = {...
file_ppi;...
filename_rp};
matlabbatch{3}.spm.stats.fmri_spec.sess.hpf = 192;
matlabbatch{3}.spm.stats.fmri_spec.fact = struct('name', {}, 'levels', {});
matlabbatch{3}.spm.stats.fmri_spec.bases.hrf.derivs = [0 0];
matlabbatch{3}.spm.stats.fmri_spec.volt = 1;
matlabbatch{3}.spm.stats.fmri_spec.global = 'None';
matlabbatch{3}.spm.stats.fmri_spec.mthresh = 0.8;
matlabbatch{3}.spm.stats.fmri_spec.mask = {''};
matlabbatch{3}.spm.stats.fmri_spec.cvi = 'AR(1)';
%--------------------------------------------------------------------------
% GLM PPImodel估计
%--------------------------------------------------------------------------
matlabbatch{4}.spm.stats.fmri_est.spmmat(1) = cfg_dep('fMRI model specification: SPM.mat File', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
matlabbatch{4}.spm.stats.fmri_est.write_residuals = 0;
matlabbatch{4}.spm.stats.fmri_est.method.Classical = 1;
%--------------------------------------------------------------------------
% PPI contrast
%--------------------------------------------------------------------------
matlabbatch{5}.spm.stats.con.spmmat(1) = cfg_dep('Model estimation: SPM.mat File', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','spmmat'));
matlabbatch{5}.spm.stats.con.consess{1}.tcon.name = 'PPI_interaction';
matlabbatch{5}.spm.stats.con.consess{1}.tcon.weights = [1 0 0 0];
matlabbatch{5}.spm.stats.con.consess{1}.tcon.sessrep = 'none';
matlabbatch{5}.spm.stats.con.delete = 0;
%--------------------------------------------------------------------------
% SPM run
%--------------------------------------------------------------------------
spm_jobman('run',matlabbatch);
end
五、批量移动结果,准备统计分析
%时间:20220919
%程序员:wishing
%批量复制结果文件
%--------------------------------------------------------------------------
% 设置循环
clear,clc;
%————————————————————————————————
for i=[2:4,6:7,11:13,15:22,24:31,34:37]
p(i) = string(num2str(i, "%02d"));%生成“01”“02”一直到“37”的字符串,剔除1, 5, 8, 9, 10, 14, 23, 32, 33)
end
p = rmmissing(p);%删除跳过产生的缺失值
%----------------------------------------------------------------
% 设置文件夹
%————————————————————————————————
for n=1:length(p)
sourcepath=string(append('D:\desk\speed\data_1\Dp_spm\PPI\RSBA18\RW\sub_1',p(n),'\con_0001.nii'));
targetpath=string(append('D:\desk\speed\data_1\Dp_spm\PPI\RSBA18\RW\RW28CON\sub_1',p(n),'.nii'));
copyfile(sourcepath,targetpath);
end