fMRI数据处理课件

功能磁共振成像（fMRI）数据分析方法

大连海事大学信息科学技术学院冯士刚前言在处理fMRI数据处理之前，我们必须修过以下课程或者有相应基础：（按重要性排列）1.概率论和统计学如果你不懂概率论和统计学，那么fMRI数据处理对你来说简直就是“天书”。（关于这一点，只在本科学过一点数理统计，远远不够用）2.计算机编程如果没有强大的计算机编程技能是不可能成为一个fMRI数据的有效使用者的。存在很多种有助于fMRI分析的计算机语言：matlab、unix、vb。（曾经雄心勃勃地打算自学matlab，可是那本厚重的<matlab宝典>还静静地躺在书架上......，自学很难）3.线性代数从数理统计到图像处理，线性代数对于fMRI分析的各个方面都和重要。想要更深刻地理解fMRI分析，必须要有线性代数的基本知识。（只知道一般线性模型不够！）4磁共振成像要完全理解fMRI数据，就要知道数据从何而来，实际上测量的是什么。这在理解MRI伪迹如何影响fMRI数据时尤为重要。（原子核、电磁场等物理学知识要涉猎！）5神经生理学和生物物理学

fMRI信号来自于对单个神经元活动的间接测量，神经元如何编码信息、从而引起血流信号的改变等等.....（fMRI是基于BOLDsignal）6信号和图像处理（高等数学）要有对信号和图像处理方法的基本了解，特别是对傅里叶分析的理解。（傅里叶分析啊....总是在papers中看到...可惜不懂......）

精心钻研十年！才能成为一个领域的专家！在脑功能成像（fMRI）分析领域尤是如此！fMRI研究框架实验设计 被试招募与 扫描科学问题结果解释

实验假设 数据统计 分析SPM,AFNI,FSL,VoxBo功能磁共振成像技术物理基础

——核磁共振（nuclearmagneticresonance）生物物理机制

——血氧水平依赖（bloodoxygenleveldependent,BOLD）核磁共振原理

半数以上的原子核具有自旋，旋转时产生一小磁场。当加一外磁场，这些原子核的能级将分裂，既塞曼效应。在外磁场B0中塞曼分裂图：氢原子+-原子核电子质子人体组织内的质子存在状态质子的运动：进动频率0=090射频脉冲纵向弛预自旋-晶格弛预T1弛预横向弛预自旋-自旋弛预T2弛预人体正常脑组织的T1、T2驰预时间驰预时间(ms)脑白质脑灰质脑脊液颅板板障T1

7809203000-260T290100300-84T1WIPDWIT2WI脑功能成像技术Deoxy-HbDecreaseBOLD原理合氧血红蛋白和脱氧血红蛋白的比率变化SPM软件包简介主要的脑功能分析软件：SPM：http://www.fil.ion.bpmf.ac.uk/spm

运行平台：MATLAB(WINDOWS&Linux)

数据格式：Analyze，MINC，ecatT

注：在最新的SPM5\spm8中将只采用NIFTI格式AFNI：/afni

运行平台：Linux，Unix，MacOSX

数据格式：AFNI，BRIK，HEADSPM以其标准化的操作，较为完备的理论成为国际上通行的脑功能分析软件。SPM软件包简介

SPM是由它的理论和思想最初是由英国的KarlFriston在英国医学研究理事会加速器中心对PET功能神经显像数据进行常规统计分析时提出并形成的。SPMClassicSPM94SPM95SPM5SPM96SPM99SPM2……SPM8SPM软件包安装AddFolder

(andnotwithSubfolders)RemoveolderSPMversionsfromthepathSPM软件包简介空间预处理面板模型建立和统计推断面板一般功能面板SPM软件包简介SPM的处理流程：校准标准化平滑统计参数图一般线形模型数据处理前准备数据格式转换。脑图的左右显示问题。spm_flip_analyze_images.m中flip=1or0Spm以前的版本，spm_defouts.m文件中

defaults.analyze.flip=1or0SliceTiming运行matlab,执行命令：spm

fmri，这样将打开spm8的操作界面。点击“SliceTiming”.在spm8，每一步处理都采用了直观的“树形结构”面板,

如果一个分支项右面有<-X号，你必须为之指定选项(否则不能运行该tree)，分支项的选项在其右侧面板指定，而帮助信息则在下面的面板中显示。Data，点击data并在下面的面板中点击“newsession”,这样在data下会出现“session”的分支项，选中该项并点击面板下方的“selectfiles”，然后用spm文件选择器选择你要处理的数据，最后点击“down”。NumberofSlices，我们输入“36”TR，我们输入“2”TA，我们输入“2-2/36”Sliceorder，我们输入“1:2:35,2:2:36”ReferenceSlice,我们输入“35”FilenamePrefix，采用默认设置。最后点击面板上方的向右的绿色三角即开始运行。运行完后将会生成一系列a*.img文件，这就是时间校准后的数据。

Realignment我们在预处理面板校准选项中选择“Realign(Est&Res)”，出现如左侧对话框，我们按下面设置进行.选中“data”,选择“NewSession”,然后选中data下出现的“Session”选项。点击“SpecifyFiles”，用spm文件选择器选择刚做完时间校准的图像（a*.img）。点击上方绿色的三角开始运行。

Normalise

在预处理面板标准化选项中选择“Normalise:Estimate&Write”。选中“data”——“newsbject”，在data下新出现的“subject”选项中作如下设置，“sourceimage”选择校准步骤中生成的mean文件，“imagetowrite”选择所有刚进行完校准的文件“ra*.img”，“templateimage”我们选择“EPI.nii”，其余采用默认设置，点绿三角运行。

注意，头动信息和空间标准化的图形文件将以spm_“data”.ps的形式保存于matlab的工作目录下，如我们是2009年4月30日处理的数据，则将以spm_2009Apr30.ps文件存于matlab的work目录下。Spm8进行统计分析时会改变工作目录，而spm8产生的其它结果文件也会以相同的名称存放在工作目录中。如果我们同一天逐个处理多个被试，请处理完一个被试后记住将该文件拷贝至被试数据目录的result文件夹中进行备份。Smoothing我们在“imagetosmooth”选项中选择所有刚进行完标准化的文件“wra*.img”然后点绿三角运行即可。FWHM我们采用默认设置“888”。

TimeBOLDsignalTimesinglevoxeltimeseriesspm建模分析ModelspecificationParameterestimationHypothesisStatisticSPMModelSpecification:

一般线形模型

一般线形模型：如果统计资料里包含包含自变量X和连续变化的反应变量Y，且它们的关系可以用：

则它们之间的关系为线性模型。

将该公式推广：上式称为一般线形模型。ModelSpecification:一般线形模型在计算机中，往往采用的是矩阵(数组)的方法来描述数学上的概念，因此就可以表示成：ModelSpecification:一般线形模型=+yN:扫描数,p:回归量数对比向量Contrasts我们设置对比向量来选择一项我们感兴趣的效应

对比向量

c

长度为p.

cTβ

是回归参数

β的线性组合.我们不关心全脑的β

参数。cTβ=1xb1+0xb2+0xb3+0xb4+0xb5+...cT=[10000…]cTβ=0xb1

+

-1xb2+1xb3+0xb4+0xb5+...cT=[0-1100…]cT=10000000T=contrastof

estimated

parametersvariance

estimateb1=cTb>0?b1

b2

b3

b4

b5...T-test

–

一维对比

–SPM{t}Question:Nullhypothesis:H0:cTb=0Teststatistic:T-contrastinSPMResMSimagecon_????imagebeta_????imagesspmT_????imageSPM{t}Foragivencontrastc:T-test:简单例子Q:在听词任务中激活区域在哪里

?cT=[10]Nullhypothesis:

听词任务versus静息SPMresults:HeightthresholdT=3.2057{p<0.001}Statistics:p-valuesadjustedforsearchvolumeset-levelcpcluster-levelpcorrectedpuncorrectedkEvoxel-levelpFWE-corrpFDR-corrpuncorrectedT(Zº)mmmm

mm0.000100.0005200.0000.0000.000

13.94

Inf0.000-63-27150.0000.000

12.04

Inf0.000-48-33120.0000.000

11.82

Inf0.000-66-2160.0004260.0000.0000.000

13.72

Inf0.000

57-21120.0000.000

12.29

Inf0.000

63-12-30.0000.000

9.89

7.830.000

57-3960.000350.0000.0000.000

7.39

6.360.000

36-30-150.00090.0000.0000.000

6.84

5.990.000

510480.00230.0240.0010.000

6.36

5.650.000-63-54-30.00080.0010.0010.000

6.19

5.530.000-30-33-180.00090.0000.0030.000

5.96

5.360.000

36-2790.00520.0580.0040.000

5.84

5.270.000-454290.01510.1660.0220.000

5.44

4.970.000

4827240.01510.1660.0360.000

5.32

4.870.000

36-2742Designmatrix0.511.522.510203040506070801voxel-levelpuncorrectedT(Zº)mmmm

mm

13.94

Inf0.000-63-2715

12.04

Inf0.000-48-3312

11.82

Inf0.000-66-216

13.72

Inf0.000

57-2112

12.29

Inf0.000

63-12-3

9.89

7.830.000

57-396

7.39

6.360.000

36-30-15

6.84

5.990.000

51048

6.36

5.650.000-63-54-3

6.19

5.530.000-30-33-18

5.96

5.360.000

36-279

5.84

5.270.000-45429

5.44

4.970.000

482724

5.32

4.870.000

36-2742Specify1stlevel

即以前版本的“fMRImodel”，spm5和spm8种的分析选项有所变化,使用Specify1stlevel做单个被试(singlesubject)分析；使用Specify2ndlevel做组分析(groupanalysis)。

选择Directory指定一个文件夹存放结果数据

。“Unitsfordesign”选择“Seconds”“Interscaninterval”输入“2”

选择“DataandDesign”后选择“NewSubject/Session”，再选择新出现的“Subject/Session”，选择“Scans”并用文件选择器选择相应任务的所有平滑后的功能图像（swra*.img）然后点击“down”.选择“Condition”后选择“”，然后选中新出现的“Condition”，“name”选项输入“nc”（表示数值一致刺激），“onset”输入“nc”（早已经做好的横向量，代表数值一致刺激启动的扫描数）选中“Durations”并输入“0”。再次选中“Condition”——“Newcondition”——“Condition”定义数值不一致刺激条件，名为“nic”，启动向量也为“nic”，同样定义数值中性刺激条件喝误差项。设置完毕后点击绿三角运行。

结束后将在选定目录中生成spm.mat文件。estimate很简单，我们只须选择刚生成的“spm.mat”文件点击“down”然后点击绿三角运行即可。结束后将生成一系列beta脑图。

results

选中“t—contrasts”，点击“definenewcontrast”,“contrast”定义为“1”，名称为“nc”,不校准，阈值取0.005。同样nic定义为“01”，中性任务nn定义为“001”，“nc-nn”定义为“10-1”，“nic-nn”定义为“01-1”，“nic-nc”定义为“-11”，其余设置一样。图为数值一致性任务“nc”的结果。

SPM结果显示注意我们指定或输入一系列的参数来进行每一步的处理，完成后我们可以通过Save按钮将每一步存为一个*.mat文件。以后我们可以通过Load按钮重新加载并使用这些*.mat文件，我们适当修改后（例如改变所运行的数据集）再选择Run按钮运行。

数据批处理（Batching）简介

BatchingSPMSPM99,SPM2[1]/wiki/SPM-Example_batch_script

[2]http://imaging.mrc-cbu.cam.ac.uk/imaging/SpmBatch5

CallSPMfunctionsdirectly[1,2]:

spm_realign

spm_slice_timing

spm_normalise

spm_write_sn

spm_smooth

spm_fmri_spm_ui

spm_spm

spm_FcUtil

spm_contrasts

Requiresknowledgeaboutinput/outputargumentsofmainSPMfunctionsandabouttheSPM.matstructure.SPM8BatchEditor:

“matlabbatch”structure.Dependencies.SPM5JobManager:

“jobs”structure.

Virtualfiles.Spm组分析组分析设计矩阵脑功能分析功能脑影像分析功能的空间定位功能交互作用功能连通有效连通脑连通的分析的分析方法：脑的结构组织所体现的功能特征主要表现在功能分区和功能整合。功能分区的形成是因为脑中功能相近的细胞和组织的物理位置倾向于聚集在一个小的区域，这使得某个认知功能的能量消耗及传输时间最小。分离的功能区之间的信息传递是通过有效连通进行的。脑功能连通的分类：脑功能连通功能连通：空间上分离的部位，在神经生理活动中的相互关系。有效连通：一个神经系统中一个基元对另一个基元的影响。有效连通分析：结构方程建模（SEM）

心理生理交互作用模型

（PPI）多元自回归（MAR）DCM心理生理交互作用模型将一个区域的激活水平与第二个区域的激活水平进行回归，那么回归的斜率将反映第二个区域对第一个区域的影响.如果在另一个实验参数条件下重复回归，这个斜率会产生变化，那么这种变化就是心理生理交互作用.基于一般线性模型(generallinearmodeling，GLM)

y=(V1A).b1+V1.b2+A.b3+e先设定一个感兴趣区域作为源区域，将成像图像的每一个体素都用实验条件、感兴趣区域的时间进程以及它们两个的交互作用进行建模.如果建模分析结果表明存在显著的交互作用，就提示激活区域和源区域之间存在受到实验条件调制的功能连通性.DCM模型简介节点状态方程：Z描述的是神经元的状态。u代表输入。θ代表我们要估计的参数。可得方程：

该方程说明了神经元状态的变化取决于神经系统内部神经结点自身和其他节点的状态，以及试验带来的扰动。DCM模型简介我们将上述方程进行双线性近似，得：其中：A代表内在耦合参数，与输入信号无关。B代表双线性耦合参数，由试验输入引起的输入的有效改变。C输入对神经元活动的外在影响。Photic:所有不同性质的光点,为第一类输入;Motion与Attention:分别为运动光点、加入注意成分的运动光点,为第二类输入。粗箭:第一类输入激活V1区;细箭:不同皮层区的内在连接;实线:第二类输入对内在连接的影响;长箭对应数字:皮层区间的内在连接强度;实线对应数字:第二类输入对有效连接的影响程度;粗箭对应数字:第一类输入对初级视觉皮层(V1)的影响程度。motion、attention施加的偶联效益:V1→V5内在连接本身不显著,motion增强了V1→V5的正向连接;attentiom

增强了IFG(额下回)→SPC(上顶叶皮层)、SPC→V5区的负反馈连接基于相关分析的脑功能连通性分析法

刺激函数分析法时间簇分析法基于内部条件的脑区间协方差分析法MR脑功能成像数据三维数组的存储形式，三个维度分别表示X、Y和Z轴。数组中每一个数称为一个体素，表示脑中一小块长方体的激活强度A

B

CDEFa1a2a3b2b1b3c3c2c1刺激函数分析法时间簇分析法基于感兴趣区的有效连接分析

基于感兴趣区的有效连接分析是通过信号值的相关分析建立脑回路的功能调节模型。其具体方法是通过对比同一任务中所有被测试者的不同感兴趣区的任务相对于控制的变化，找出这些感兴趣区之间的关系。任务过程中的时间协方差分析

时间依赖性连接的动态分析

选择两个感兴趣区，计算每个图像扫描的时刻，感兴趣区的相关系数，因此在整个实验过程中，可组成一组两个感兴趣区相关系数变化的序列，动态的分析整个过程中的脑区间的联系。设S(x0,y0,z0,n,i)为第i个被测试者某一感兴趣区t值最强点(x0,y0,z0)的S值，则对于第n幅脑图，由不同被测试者可组成如S(x0,y0,z0,n,1),S(x0,y0,z0,n,2),……S(x0,y0,z0,n,I)的样本，将它与另一感兴趣区t值最强点(x1,y1,z1)的样本S(x1,y1,z1,n,i)作相关分析，得到第n幅脑图时两个感兴趣区相关系数。将任务时间段内的相关系数取出，并以次绘出线图；将控制时间段内的相关系数平均，作为绘图的基线。VBMProcessingResting-statefMRI:

静息态数据处理Outline预处理

ALFF,ReHo,FC计算统计结果呈现预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08数据整理被试信息整理数据整理被试信息整理数据整理被试信息整理原始数据整理静息功能数据结构数据DTI数据数据整理被试信息整理原始数据整理静息功能数据结构数据DTI数据预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08DICOM->NIFTIMRIcroN’sdcm2niiguiSPM5’sDICOMImport预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08SliceTimingWhy?SliceTimingWhy?Huetteletal.,2004SliceTiming1:2:33,2:2:323322-(2/33)33SliceTiming预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08RealignWhy?RealignRealign检查头动?b=load(‘rp_name.txt’);c=max(abs(b));c(4:6)=c(4:6)*180/pi;预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08NormalizeWhy?Huetteletal.,2004mean_name.imgr*.imgEPI.nii-90-126-72;9090108333Normalize预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08SmoothWhy?减少被试间的空间结构差异；有利统计检验的有效性;w*.img平滑核Smooth预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08去线性漂移预处理数据整理

DICOM->NIFTI

去前10个时间点

SliceTimingRealign

NormalizeSmooth

去线性漂移

Filter:0.01-0.08滤波Why?关心低频的活动滤波Outline预处理

ALFF,ReHo,FC计算统计结果呈现文献管理ALFFZangetal.,2007ALFFReHo注意：预处理时不做平滑，计算完ReHo后再进行平滑。Zangetal.,2004ReHo功能连接提取协变量

头动参数:rp_name.txt

全脑信号、脑脊液、白质信号提取协变量提取协变量提取协变量提取协变量提取协变量提取协变量提取协变量功能连接提取协