[spm操作] 什么是ROI，如何做ROI以及批量提取ROI的%signal change的示例程序

本帖作为《用Matlab和SPM批量处理被试的经验总结》的一部分
目录贴请见http://home.52brain.com/forum.ph ... =1&extra=#pid158525

ROI (region of interest, 感兴趣区)，顾名思义，就是研究者感兴趣的地方。通常在脑成像数据处理中，会把某些ROI中的信号提取出来，来和行为或是其它脑区数据做相关。做静息态的朋友也应该了解，所谓的功能连接(functional connectivity )，就是把某个脑区做成ROI，也叫种子点(seed region)，提取出其信号随时间点变化的序列（timecourse），然后用这个序列去和全脑体素做相关（据我所知，至少有时候是简单的皮尔逊积差相关），相关较强的脑区就是与这个种子点功能连接显著的脑区。
在某些有任务的研究中，研究者会提取某个脑区的信号（beta weights， contrast value 或者% signal change），然后像对待行为数据那样使用它们。


如何选择或定义ROI
选择ROI一般基于（1）特定的脑结构，
（2）前人文献中显著的脑区，
（3）元分析结果，
（4）定位run （localizer run，同一批被试在另一个任务相同或相似的单独的run）中激活的脑区

可以参考Poldrack, R. A. (2007). Region of interest analysis for fMRI. Social cognitive and affective neuroscience, 2(1), 67-70.

基于脑结构的（比如杏仁核）叫做 structural ROI。其它的叫做functional ROI。具体例子可以看MarsBar的手册。

在实际的定义时，对于结构ROI，可以根据现成的模板（template）来定义，比如AAL的各个脑区。
对于功能ROI，可以提取激活峰值点的坐标，做成球形ROI（spherical ROI），也可以根据激活团块提取激活区ROI。实际上，你可以定义任何形状的ROI（当然还要考虑体素本身的颗粒感）。具体的定义方式可以参考MarsBar的手册。关于提取ROI内的信号，其实ROI就是一个小一点的mask，mask就是大一点的ROI，或者说它们本质上一样。ROI就是一个大的三维矩阵，除了对应于你的感兴趣的脑区之外的地方都是0，ROI之内都是1（在SPM中看的时候ROI边缘介于0和1，但实际上图像本身只有1和0两个值。）。
如下图就是一个球形ROI的示例：
 

定义ROI或进行分析常用的工具包
我个人最常用的是MarsBar
此外还有WFU PickAtlas
这两个工具包都有详细的说明和示例。


做ROI以及提取% signal change和contrast value的批处理示例程序
（基于Marsbar，所以要把MarsBar的工具包纳入到matlab的搜索路径中）
#########################################################
% MarsBaR的批处理程序，将结果保存为txt和mat两种格式。文件题目标明roi是取对比值还是信号百分比，半径，中心和保存文件的时间（以便于用统计软件进行批处理）
clear;

marsbar('on');   % 打开marsbar
subjroot = 'I:\kongliliushuang\SPM_MRIresults_Analysis\StruNorm_tTest\';   % 数据所在的路径
roi_dir = ('I:\kongliliushuang \SPM_MRIresults_Analysis\ROI_StruNorm_tTest\');% 存放ROI的路径
roiimg_dir=fullfile(roi_dir,'img');%存放ROI图像的路径
spm('defaults', 'fmri');     %设置SPM默认参数
sphere_matr=[3,6,10];%三种球形ROI的半径，分别为3mm，6mm以及10mm，你可以只用一种，也可多试几种
              %以下输入要取的ROI的中心（可以是其他研究中激活点的峰值，或是定位run的激活脑区峰值）
fociMatrix=[63 -54 39;
             51 -48 27;
             63 -54 24];
    cd (roi_dir);%保存结果文件的地方
    commonstring='myROI'; %保存的几个文件共有的开头的名字，可以随意更改
    save_name=[commonstring,datestr(now,'yyyy_mm_dd_HH_MM_SS'),'.mat'];  %保存综合文件的名字
    save_name1=[commonstring,'_perc_',datestr(now,'yyyy_mm_dd_HH_MM_SS'),'.txt'];  %保存信号百分比文件的名字
    save_name2=[commonstring,'_contr_',datestr(now,'yyyy_mm_dd_HH_MM_SS'),'.txt'];  %保存对比值文件的名字
    save_name3=[commonstring,'_centermatrix_',datestr(now,'yyyy_mm_dd_HH_MM_SS'),'.txt'];  %保存ROI球心矩阵的文件的名字


  summary_sphere_result_perc=[]; %设置一个空矩阵来装所有的结果
  summary_sphere_result_con=[];  

for loop_sphere=1:3   %使用三种半径

summary_result_perc=[]; %设置一个空矩阵来装所有的结果
  summary_result_con=[];


for loop=1:size(fociMatrix,1)

%%%%%%%%%%%%%%%%%%%%%%%%%%设置球形ROI%%%%%%%%%%%%%%%%%%%%%%%

sphere_centre = fociMatrix(loop,:); %球形ROI的中心  此处必须是MNI坐标，如果不是，需要转化为MNI坐标
sphere_radius= sphere_matr(loop_sphere);         %球形ROI的半径
sphere_roi = maroi_sphere(struct('centre', sphere_centre,'radius', sphere_radius));
% 给roi命名
sphere_roi = label(sphere_roi,sprintf( 'sphere%2dmm_%2d_%2d_%2d_roi',sphere_radius,sphere_centre(1),sphere_centre(2),sphere_centre(3)) );
roi_name=sprintf('sphere%2dmm_%2d_%2d_%2d_roi.mat',sphere_radius,sphere_centre(1),sphere_centre(2),sphere_centre(3));
saveroi(sphere_roi, fullfile(roi_dir,roi_name));% 将ROI保存为MarsBaR的ROI文件
% 将ROI保存为图像
save_as_image(sphere_roi, fullfile(roiimg_dir,sprintf('sphere%2dmm_%2d_%2d_%2d_roi.img',sphere_radius,sphere_centre(1),sphere_centre(2),sphere_centre(3))));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

roi=sphere_roi;

%% 
nsub=20;      %被试的个数
result_perc=[]; %设置一个空矩阵来装所有的结果
result_con=[];
  for sub=1:nsub
            if sub==3|| sub==7 || sub==17;   %跳过有问题的被试
            continue;
           end
    dir=[ subjroot sprintf('subject%02d',sub)];

    D = mardo(fullfile(dir, 'MyResults', 'SPM.mat'));   %此处最重要，将SPM的设计矩阵等信息载入进来

    Y = get_marsy(roi, D, 'mean'); % 从模型中取出ROI中的数据
    E = estimate(D, Y); % MarsBaR 进行估计
           %以下取%signal
        %change+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            [e_specs,e_names]=event_specs(E);  %定义模型中所有的事件，即条件
             n_events=size(e_specs,2);
             dur=9;  %所提取事件的持续时间

              for e_s=1:n_events
                 pct_ev(e_s)=event_signal(E,e_specs(:,e_s),dur) %返回% signal change
              end

             perc(1)=(pct_ev(2)+pct_ev(7)+pct_ev(12)+pct_ev(17))/4;  %由于每个run中都有我要研究的条件，将每个run中的此条件提取出%信号变化，然后求均值。2,7,12,17表示在整个设计矩阵中这个条件的排序，可参考http://home.52brain.com/thread-25155-1-1.html 中的设计矩阵图片
             perc(2)=(pct_ev(3)+pct_ev(8)+pct_ev(13)+pct_ev(18))/4;
             perc(3)=(pct_ev(4)+pct_ev(9)+pct_ev(14)+pct_ev(19))/4;
             perc(4)=(pct_ev(5)+pct_ev(10)+pct_ev(15)+pct_ev(20))/4;
             perc(5)=(pct_ev(6)+pct_ev(11)+pct_ev(16)+pct_ev(21))/4;
            result=perc;
            result_perc=[result_perc; result]
            %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            %以下取constrast值++++++++++++++++++++++++++++++++++++++++++++++++++++

            [E Ic]=add_contrasts(E,D); %取出模型中的所有对比
           stat_struct=compute_contrasts(E,Ic)
             conva(1)=stat_struct.con(1);   % 取出所有前四个对比的contrast值
             conva(2)=stat_struct.con(2);
             conva(3)=stat_struct.con(3);
             conva(4)=stat_struct.con(4);
             resultc=conva;
             result_con=[result_con;resultc]
             %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  end

  summary_result_perc=[summary_result_perc,[fociMatrix(loop,:),[sphere_radius*ones(1,2)];result_perc]];
  summary_result_con=[summary_result_con,[fociMatrix(loop,:),[sphere_radius*ones(1,1)];result_con]];

end

summary_sphere_result_perc=[summary_sphere_result_perc;summary_result_perc];
summary_sphere_result_con=[summary_sphere_result_con;summary_result_con];
end
cd (roi_dir);%保存结果文件的地方
save(save_name,'fociMatrix','summary_sphere_result_perc','summary_sphere_result_con');
save(save_name1,'summary_sphere_result_perc','-ASCII')
save(save_name2,'summary_sphere_result_con','-ASCII') 

save(save_name3,'fociMatrix','-ASCII')   % 保存结果。可以直接拷贝到excel或用其他软件处理。输出的txt文档是一张大表，你要搞清楚不同ROI和半径大小是怎么排列的。

你把生成的txt文件中的东西贴到spreadsheet中看起来就是下图这个样子的（当然这个截图只是一个小部分）。为了便于查看，我把有ROI坐标的那一行的格式换成了整数并且高亮成了黄色。可以看出，我从坐标为[21 12 31]的ROI中提取出了6个条件的信号（横着看有6列），并且表格中21 12 21 后面的3 3 3表示ROI半径是3mm。而且能看出我提取了17个被试的数据。当然整个表格可以非常大，你也可以根据自己的目的来修改代码已生成更符合你需要的排列方式。