PPI（Psycho-Physiologic Interaction）分析是检验某个
种子点脑区
和
其它脑区
的功能连接的强度是否受到不同实验条件影响（在不同的实验条件下不同）的一种较复杂的分析方法。
结合我之前帖子里的情绪面孔的实验，可以这样具体说明：我想看杏仁核和其它哪些脑区的功能连接强度会在
看正性面孔时
以及
看中性面孔时
有不同。
或者说我想看正性情绪是否调节了杏仁核和其它哪些脑区的连接。这个分析得出的结果是全脑的激活，如果内侧前额叶有激活，就表明，
正性情绪时
杏仁核
和
内侧前额叶
的连接强度
不同于
中性情绪时
杏仁核
和
内侧前额叶
的连接强度
。至于如何不同，这个要画图才能准确看出来。具体的步骤和画图操作可以看《SPM8 Manual》的第33章。这一章中的例子是多个run的情况，手册建议将多个run连接在一起并加上代表每个run的regressor，便于PPI的分析。
之前我的帖子介绍了批量进行第一层统计分析的代码，这个PPI分析要使用第一层分析生成的文件，所以要先跑完第一层的统计分析才能运行下面的PPI分析的代码。
%%%%%%%%%%%%%%%%% 处理程序开始 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%本程序用来进行SPM中第一层的PPI的分析
%by 空里流霜， Aug, 2019
% 先说说我的文件夹的基本安排。我有个文件夹叫kongliliushuang,里面有子文件夹scripts，BehavioralData和ImagingData；%scripts文件夹里有imaging文件夹，用来放处理脑成像数据分析的代码，里面也有做好的一些mask；
%BehavioralData里面有个文件夹preprocessed，存放整理好的整洁的行为数据。
%imagingData文件夹中包含preprocessing文件夹，放置
预处理的各步骤的文件夹，其中就有
文件夹FunImgSWRAF，里面存放每个被试编号命名的文件夹，每个被试命名的文件夹里是文件夹RUN1
(假设每个被试只有一个run)
，再里面就是swraf打头的预处理过的脑成像文件。
clear;
subjIDs=load('/kongliliushuang/subjectID_list.txt');
%我有一个txt文件，里面是实验中所有被试编号的列表（只有数字，不能有文字或字母）。同一个文件夹下还有txt文件，里面装的是我想删除的被试列表
subjIDs_exc1=load('/kongliliushuang/subjectID_list_fMRIExcludeBadImaging.txt');
% load all the subjects you want to exclude
subjIDs_exc2=load('/kongliliushuang/subjectID_list_fMRIExcludeHeadMotion.txt');
% load all the subjects you want to exclude
subjIDs_exc3=load('/kongliliushuang/subjectID_list_fMRIExcludeBadBehavioral.txt');
% load all the subjects you want to exclude
subjIDs_exc4=load('/kongliliushuang/subjectID_list_excluded.txt');
% load all the subjects you want to exclude
%下面用几个for语句逐一删除不需要的被试
for m=1:size(subjIDs_exc1,1)
indexRow=find (subjIDs(:,1)==subjIDs_exc1(m));
subjIDs(indexRow,:)=[];
end
for m=1:size(subjIDs_exc2,1)
indexRow=find (subjIDs(:,1)==subjIDs_exc2(m));
subjIDs(indexRow,:)=[];
end
for m=1:size(subjIDs_exc3,1)
indexRow=find (subjIDs(:,1)==subjIDs_exc3(m));
subjIDs(indexRow,:)=[];
end
for m=1:size(subjIDs_exc4,1)
indexRow=find (subjIDs(:,1)==subjIDs_exc4(m));
subjIDs(indexRow,:)=[];
end
%%%####################################
filename1='myEmotionExperiment';
%实验的名字
dirName='positiveEmotion';
%具体想做的分析的名字
addpath('/kongliliushuang/scripts/imaging/');
%预处理用的代码放置的位置，里面也有事先制作的mask，在下面会用到。
mkdir(['/kongliliushuang/ImagingData/analysis/',filename1,'-',dirName]);
cwd=['/kongliliushuang/ImagingData/analysis/',filename1,'-',dirName,'/'];
datadir='/kongliliushuang/ImagingData/preprocessing/FunImgSWRAF/';
% 预处理的文件SWRAF的位置
maskfile='/kongliliushuang/ImagingData/brainmask0.2.img';
%explicit mask的路径，这个是事先自己制作好的，制作方法可以参考我这个帖子http://home.52brain.com/thread-26070-1-1.html。当然也可以选择默认的implicit的mask，但那样会遮挡部分大脑边缘脑区。所以这里用了explicit mask。
mulcondpath=['/kongliliushuang/ImagingData/fMRI_mat_files/',filename1];
% 制作设计矩阵需要的每个事件的names, onsets和 durations，以及contrast所在的位置，事先要做好，代码可以在这里找到http://home.52brain.com/forum.ph ... =1&extra=#pid163180
subExpID=subjIDs(:,1);
% 已经筛选好的要进行第一层分析的被试
%-----------分析各步骤的开关------
creatFolder=1;
showIndiResult=0;
%%%%%%%%%%%%%%%%%%上面是第一层PPI分析的预备步骤%%%%%%%%%
%%%%%%%%%%%%%%% 生成放置第一层分析结果的文件夹 %%%%%%%%%%%%%%%%%
if creatFolder
cd(cwd)
mkdir('PPIsubResults')
%
生成放置PPI分析结果的文件夹，名字要和传统的第一层分析不同
cd ('PPIsubResults')
for wnumfolder=1:size(subExpID,1)
mkdir (num2str(subExpID(wnumfolder)));
end
end
%end of switch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spm_get_defaults
% set defaults
global defaults
spm('defaults','FMRI')
nses=1; %
run的个数，这里只有一个。如果有多个run，需要连接（concatenate）在一起，并且加入哑变量来作为regressor。这种情况下进行第一层分析的编程时要调整每个事件的onset，并且头动也要连接在一起。具体的操作并不复杂。
spm_jobman('initcfg');
for sub=1:size(subExpID,1)
resultpath1=[cwd,'subResults/',num2str(subExpID(sub)),'/'];
%传统第一层分析的路径，里面是每个被试的文件夹
resultpath2=[cwd,'PPIsubResults/',num2str(subExpID(sub)),'/'];
%PPI第一层分析的路径，里面是每个被试的文件夹
cd (resultpath2);
%%%%%%%%%%%%%为PPI分析设置contrasts %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jobs{1}.stats{1}.con.spmmat=cellstr(fullfile(resultpath1,'SPM.mat'));
%打开传统第一层分析生成的SPM.mat
jobs{1}.stats{1}.con.consess{1}.tcon.name ='PositiveEmotion';
%感兴趣的效应
jobs{1}.stats{1}.con.consess{1}.tcon.convec ={[1 0 -1]};
% 第一层的设计矩阵中第一列是正性表情，第三列是中性表情
jobs{1}.stats{1}.con.consess{1}.tcon.sessrep = 'none';
jobs{1}.stats{1}.con.consess{2}.fcon.name ='Perception';
%这里设置一个f contrast，方便下面PPI调整数据
jobs{1}.stats{1}.con.consess{2}.fcon.convec ={[1 0 0; 0 1 0; 0 0 1]};
% 第一层的设计矩阵中第一列是正性表情，第二列是负性表情，第三列是中性表情
jobs{1}.stats{1}.con.consess{2}.fcon.sessrep = 'none';
jobs{1}.stats{1}.con.delete = 1;
if 0
jobs{1}.stats{2}.results.spmmat=cellstr(fullfile(resultpath1,'SPM.mat'));
jobs{1}.stats{2}.results.conspec.contrasts=Inf;
jobs{1}.stats{2}.results.conspec.threshdesc='none';
% 多重比较校正'none','FWE'or'FDR'
jobs{1}.stats{2}.results.conspec.thresh=0.01;
% p值
jobs{1}.stats{2}.results.conspec.extent=0;
%show cluster with how many number if voxels
end
spm_jobman('run',jobs);
clear jobs;
%%%%%%%%%%%%% 设置VOI并提取信号 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jobs{1}.util{1}.voi.spmmat =cellstr(fullfile(resultpath1,'SPM.mat'));
% 打开传统第一层分析生成的SPM.mat
jobs{1}.util{1}.voi.adjust = 1;
% 如果不调整效应（就是控制设计矩阵中的各个事件的影响），此处写0;
% 如果根据所有其它条件调整，此处写NaN；
% 如果有选择地调整，就需要用到上面设置的F-contrast。这里写的数字是这个F contrast在SPM.mat中的顺序，可以通过SPM的results按钮打开查看。上面生成的F-contrast应该是SPM.mat中保存的唯一的一个contrast (因为jobs{1}.stats{1}.con.delete = 1;)，所以此处可以写1.
% 如果你的F contrast 是[1 0 0; 0 1 0; 0 0 1], 你只控制设计矩阵中前三个事件
% 如果你的F contrast 是 [1 0 0;0 0 1], 你只控制第一个和第三个事件
jobs{1}.util{1}.voi.session = 1;
jobs{1}.util{1}.voi.name = 'Amg';
jobs{1}.util{1}.voi.roi{1}.mask.image = {'/kongliliushuang/ImagingData/analysis/ROI/img/AmygdalaROI.img,1'};
%事先做好的ROI；在SPM里，此处可以设置做球形ROI还是使用图像。因为我事先做好了ROI的图像，所以直接使用了
jobs{1}.util{1}.voi.roi{1}.mask.threshold = 0.5;
%ROI图像里ROI的部分是1，所以0.5的阈值也是可以的。
jobs{1}.util{1}.voi.expression = 'i1';
%大于0.5的部分都会成为VOI，并将要从中提取信号
spm_jobman('run',jobs);
clear jobs;
%%%%%%%%%%%%% 制作PPI交互作用文件 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jobs{1}.stats{1}.ppi.spmmat = cellstr(fullfile(resultpath1,'SPM.mat'));
%打开传统第一层分析生成的SPM.mat
jobs{1}.stats{1}.ppi.type.ppi.voi = cellstr(fullfile(resultpath1,'VOI_Amg_1.mat'));
% 上一步生成的VOI文件
jobs{1}.stats{1}.ppi.type.ppi.u = [1 1 1; 3 1 -1];
%这里的contrast形式有点复杂。
%第一列的1和3表示设计矩阵的第一列（正性表情）和第三列（中性表情）是我们要比较的两个条件；
%第二列表示某一个事件的名字在此事件中的顺序，如果没有parametric modulator，一个事件中只有一个名字，所以此处是1 1；
%第三列是contrast weight，因为之前设置的t contrast是[1 0 -1], 此处便写成 1 -1。
jobs{1}.stats{1}.ppi.name = dirName;
jobs{1}.stats{1}.ppi.disp = 1;
spm_jobman('run',jobs);
clear jobs;
%%%%%%%%%%%%% PPI的GLM分析 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jobs{1}.stats{1}.fmri_spec.dir=cellstr(resultpath2);
%放置PPI结果的路径
jobs{1}.stats{1}.fmri_spec.timing.units='secs';
%seconds='secs'; number of TR=scans
jobs{1}.stats{1}.fmri_spec.timing.RT=2; %repetition time
jobs{1}.stats{1}.fmri_spec.timing.fmri_t = 34;
jobs{1}.stats{1}.fmri_spec.timing.fmri_t0 = 9;
%可以参考下面我的预处理代码的注释
%For data collected from SIEMENS,
%the temporal middle point is 33 and 34. The spatial middle point is 17 and 18.
%It should correspond to the "microtime resolution"(also called t) of the 1st level statistical analysis
%and "microtime onset"(also called t0).If there are 33 number of slices, t0=17.
% When reference for slice timing is 33, the 33rd slice is the 17th to be collected . If the reference slice is 17, then t0=9.
%%%% 设置PPI分析的regressors %%%%%
muldirFiles1=[datadir,num2str(subExpID(sub)),'/RUN1/'];
%预处理过的图像所在的路径
mulfiles=spm_get('Files',muldirFiles1,'swr*.img') ;
% 预处理过的图像
mulmotionfiles=spm_get('Files',muldirFiles1,'rp*.txt') ;
%头动
myPPI=load([resultpath1,'PPI_',dirName,'.mat']);
% PPI交互项
jobs{1}.stats{1}.fmri_spec.sess(1).scans=cellstr(mulfiles);
jobs{1}.stats{1}.fmri_spec.sess(1).regress(1).name = 'PPI-interaction';
jobs{1}.stats{1}.fmri_spec.sess(1).regress(1).val = myPPI.PPI.ppi;
jobs{1}.stats{1}.fmri_spec.sess(1).regress(2).name = 'Amg-BOLD';
jobs{1}.stats{1}.fmri_spec.sess(1).regress(2).val = myPPI.PPI.Y;
jobs{1}.stats{1}.fmri_spec.sess(1).regress(3).name = ['Psych_',dirName];
jobs{1}.stats{1}.fmri_spec.sess(1).regress(3).val = myPPI.PPI.P;
jobs{1}.stats{1}.fmri_spec.sess(1).multi_reg = {''};
jobs{1}.stats{1}.fmri_spec.sess(1).hpf = 128;
jobs{1}.stats{1}.fmri_spec.bases.hrf.derivs = [0 0];
jobs{1}.stats{1}.fmri_spec.volt = 1;
