基于图像和基因数据整合分析的大脑功能特异性基因特征识别

王天杰

摘要：在脑科学领域，包含基因和多类图像数据的多模态数据的结合分析在揭示大脑的功能和机理的过程中越来越发挥着重要的作用。本文通过运用统计学中的偏最小二乘回归方法，结合基因数据和两类图像数据，对大脑的功能发育进行了相关分析，并查找文献进行了结果的比较，分析了与大脑功能发育相关的基因，并对功能网络进行了相关分析，得到了可靠的结果。

关键词：偏最小二乘回归（PLS） 脑图像 基因 网络

中图分类号：TP391 文献标识码：A 文章编号：1009-3044（2018）28-0200-06

1 绪论

人体各个器官中都有人体全部的基因，但是不同器官中基因的表达不同，而且不同的表达与器官的不同功能有关，基因突变会导致基因表达的变化，进而导致相关器官的功能异常。本研究的目的是研究人的不同年龄阶段以及不同脑区的大脑发育与人的全基因表达之间的关系，揭示显著表达的基因在该年龄段或该脑区大脑发育中的作用，即找到与大脑发育显著相关的基因，并认为这些基因的突变更有可能在疾病中扮演重要角色，进而发现如阿兹海默症，精神分裂症，多动症等大脑相关疾病的致病基因与这些基因的高度相关性，为致病基因的发现和疾病的治疗提供参考。为了达到这一目的，本研究需要的数据有涵盖多个年龄段的大脑的结构数据以及人体所有基因在相应年龄段大脑中的表达数据，结构数据主要是大脑皮层的生物数据，如皮层厚度，脑区体积等。

fMRI 基于血氧水平依赖（Blood Oxygen Level Dependent， BOLD）效应，通过测量局部脑区中脱氧血红蛋白浓度的改变来间接表明大脑神经元的功能活动，主要分为任务态（Task-related）和静息态（Resting-state）[1]。任务态fMRI通过事先设定好的特定任务刺激被试然后检测大脑的反应，但对于一些认知受损或认知尚未发育完成的被试来说完成任务有一定困难。而rs-fMRI 不需要执行特殊的任务，只需被试闭上眼睛，大脑清醒但不刻意做任何系统思考，克服了上述缺点[11]。

1.1 布罗德曼分区系统

布罗德曼分区最早由德国神经科医生科比尼安·布洛德曼（Korbinian Brodmann）提出。他的分区系统包括每个半球的52个区域。其中一些区域今天已经被细分，例如23区被分为23a和23b区等。

1.2 Destrieux分区系统

Destrieux分区系统是基于脑沟深度的在freesurfer中广泛应用的分割方案。这种划分方案能够自动而精确的定义脑沟和脑回。这种划分方案把大脑的每个半球划分为74个脑区，即共148个脑区。

图1 结构分区模板Destrieux

2 方法

2.1 腦图像数据

本研究所使用的脑图像数据是人脑连接组计划（HCP）的WU-Minn HCP Lifespan Pilot Data数据集中的从不同年龄段的人采集的结构磁共振数据（structure MRI）和静息态功能核磁共振影像数据（rfMRI）[10]。探索人脑是21世纪伟大的科学挑战之一。 人类连接体项目（HCP）正在通过阐明脑功能和行为的神经通路来解决这一挑战的一个关键方面。 解密这个令人惊讶的复杂接线图将揭示什么使我们成为独特的人类，什么使每个人与所有其他人不同。HCP项目由华盛顿大学，明尼苏达大学和牛津大学（WU-Minn HCP联盟）领导，联盟全面使用最先进的无创神经影像学方法，将目标总数1200名健康成年人的人脑电路进行全面绘制。 它将产生并贡献关于脑连通性，大脑连通性与行为的关系，以及遗传和环境因素对脑电路和行为的影响的宝贵信息[2]。该项目对于其数据有专门的处理流程，经数据预处理后用于分析。本数据处理的先决条件是：1，64位Linux操作系统；2）安装FSL （FMRIB Software Library）v5.0.6；3）安装FreeSurfer version 5.3.0-HCP；4）安装Connectome Workbench version 1.0；5）安装HCP version of gradunwarp version 1.0.2。6）安装HCP-pipelines。安装完成后，设置好工作路径和环境变量，然后开始处理数据。

结构预处理：Prefreesurfer：为每个人生成自己的未失真的本地结构体积空间，对齐T1w和T2w图像，进行偏置场矫正，把每个人自己的本地结构体积空间融合到标准体积空间中。本地结构体积空间有每个人大脑各个结构的体积和面积等数据，是对每个人的大脑最逼真的体现，而标准体积空间是为了同时研究不同人的大脑，使他们具有可比性，把他们都对齐到标准体积空间中。该部分的第一步是MR梯度非线性引起的畸变的校正。在结构处理中使用的所有图像（T1w，T2w和场图幅度和相位）必须校正梯度非线性失真。为了校正失真，由每个梯度线圈产生的磁场由球面谐波膨胀（特定于3T Connectome扫描仪中的SC72梯度特性）建模[8]。 然后，使用FreeSurfer（Jovicich等人，2006）中提供的gradient_nonlin_unwarp软件包的定制版本完成校正。Freesurfer：这一部分基于FreeSurfer version 5.2，且有许多改进之处。这一部分的目标是把整个大脑划分为已知的结构，重建大脑白质/灰质分界面以及软膜表面（隔离大脑灰质与脑脊液的膜），并进行freesurfer的折叠式表面配准。Postfreesurfer：生成所有的体积和表面文件，以便在connectome workbench中可视化观察。下采样融合后的表面以进行连接分析。HCP的数据中有三种表面，每个人的本地表面（～136k个顶点），高分辨率的Conte69标准表面（～164k个顶点），以及低分辨率的Conte69标准表面（～32k个顶点，适用于低分辨率数据如fMRI或弥散张量数据的跨目标分析）。

在结构预处理后，得到了每个人全脑的结构数据，为了本研究的目的，从中选出四项，即每个脑区的厚度，表面积，体积，曲度，进行后续的分析。基因数据只有16个脑区，每个脑区的左右半球的基因表达值取均值，为了保持图像数据与基因数据的一致，在进行脑区匹配后，对于一个脑区的基因表达值，在图像数据中，左右脑同一脑区的体积相加，表面积相加，厚度和曲度取均值。

功能预处理： 功能预处理在结构预处理后进行[13]。

fMRIVolume：去除空间失真，重新调整体积以弥补头部动作带来的误差，将fMRI数据与结构数据融合，减少偏置场误差，将4D图像归一化为全局平均值，并用最终的脑表面来掩盖数据。这一部分的结果输出可以用来进行基于整个大脑的fMRI数据分析。fMRIVolume预处理部分在很多方面像PreFreeSurfer预处理部分，它们解决类似的问题[6]。 需要在完成HCP结构预处理（PreFreeSurfer，FreeSurfer和PostFreeSurfer）的基础上进行。与PreFreeSurfer预处理部分一样，第一步是使用与PreFreeSurfer部分相同的方法校正梯度非线性引起的失真。接下来是重新对准时间序列，纠正接受扫描的人的头部动作引起的误差。fMRIVolume预处理部分的最后一步是把所有的融合变换和失真矫正步骤联系起来，形成单一的非线性变换，以应用在单个的重采样步骤中。

fMRISurface：这部分在fMRIVolume的基础上进行，其任务是把每个人的时间序列数据从各自的大脑转换到标准空间中。这部分的输出结果可以用来进行静息态或任务态的fMRI分析。第一步是一种新颖的部分体积加权带状约束体积对表面映射算法，其中白质/灰质交界面和软膜表面用于定义哪些fMRI体素在灰質带内。为了获得最大的准确性，在本步骤中使用本地网格表面（在MNI体积空间中），因为它们没有被重新采样，因此能够最精确地重现这个大脑中的组织轮廓。一些体素将部分在灰质色带内，部分在其外部。分配给每个表面顶点的强度值是完全或部分在灰质带内的体素的加权平均值，部分体素根据其在色带内的部分体积进行加权[9]。在每个顶点，算法找到包含白质/灰质分界面和软膜表面包含该顶点的三角形（两个表面的顶点必须具有对应关系），并将这些三角形与通过连接两个表面对应三角形的边形成的四边形相连接，构建多面体。然后，它在每个体素中创建点的网格，并测试每一个点，看它是否在多面体内。如果多面体是扭曲的而导致测试结果模糊，算法判定该点是半在内的。用于体素的权重是内部网格点的数量，加上半在多面体内网格点数量的一半。每个个体受试者的时间序列被采样到一组左右半球表面顶点的标准集，以及皮质下灰质体素的标准集中，形成标准的CIFTI灰度空间。两个半球表面的时间序列以及每个皮质下结构的时间序列被组合在尺度为灰度×时间的单个数据矩阵中。CIFTI密度时间序列是fMRI预处理部分的最终输出，代表着灰质部分时间序列数据最简洁而仍然完整的表示[12]。

在具体的研究过程中需要根据真实数据的产生方式，调节处理过程中的相关参数，以得到正确的结果。

2.2 脑基因数据

该数据集从收集自57个死后人脑的一个或两个半球的1，340个组织样本产生，从胚胎发育到成年后期，代表多种族的男性和女性[7]。该数据发现大约86%的蛋白质编码基因使用严格标准进行表达，其中超过90%的蛋白质编码基因在整个转录本或外显子水平跨区域和/或时间差异调控。 这些时空差异中的大多数发生在出生前，随后在出生后寿命期间区域转录组之间的相似性增加。这是一个一个关于不同时空人脑转录组的全面，公开的数据集。从57个人的1340个样本中提取的基因数据，涉及16个脑区，包括小脑皮层，丘脑，纹状体，扁桃体，海马和11个新皮层区域。每个样本是16个脑区之一的表达数据。这57个人分为15个年龄段，为了与图像数据一致，取最后五个年龄段，6岁～12岁，12岁～20岁，20岁～40岁，40岁～60岁，60岁以上，分别对应图像数据的8岁～9岁，14岁～15岁，25岁～35岁，，45岁～55岁，65岁～75岁。先按照年龄段划分基因数据样本，得到5个年龄段的样本数分别是42，87，220，91，93。每个年龄段所有同一脑区的左右脑基因表达数据取平均值，作为该年龄段该脑区的基因表达值。最后得到每个年龄段所有脑区的基因表达值。

2.3 脑图像数据与脑基因数据的匹配

脑图像数据是整个大脑的完整图像数据，按照脑图像数据使用的大脑模板，左右大脑半球各被划分为74个脑区，图像数据包括每个脑区的结构信息，如脑区的体积，平均厚度，表面积，曲度[5]。由于脑图像数据与脑基因数据划分大脑的模板不同，需要找到他们之间的等价关系，再把基因数据和图像数据匹配起来。以布罗德曼分区系统作为中介，分别找到图像数据和基因数据中脑区和布罗德曼脑区的对应关系，即可实现二者之间的匹配。从而得到一个脑区的基因数据和其结构数据。最后找到9个脑区的对应。

3分析

3.1 偏最小二乘回归（partial least square regression）[4]

偏最小二乘回归（又名潜在结构投影），结合并扩展了主成分分析（PCA）和多元线性回归的特点，它的目标是从一组自变量预测一组因变量。这种预测从自变量中抽取一组称为潜在变量的正交因子，而且这些潜在变量的预测能力最强。

在本研究中，定义PLS组分为基因表达权重的线性组合，该组合与MRI数据高度相关。具体的，在本研究中，自变量为基因，因变量为MRI数据，PLS的目的是抽取每个时间段与大脑结构数据高度相关的基因的权重组合即每个时间段的特定基因，以及每个脑区与此脑区结构数据高度相关的基因的权重组合，即该脑区发育过程中不同于其他脑区的特定基因。

3.2 静息态fMRI数据处理与功能脑网络的构建

在结构预处理的基础上进行rfMRI数据的预处理[3]，即功能预处理。经过两步功能预处理的步骤后，用最后得到的数据进行分析。得到标准空间中27个人的91292个位置点的时间序列数据，时间点之间的间隔为0.72秒，共采集了420个时间点。然后获取每个脑区的时间序列，具体做法是，采用HCP提供的多模态脑区划分模板（这个模板把每个半球划分为180个脑区，共360个脑区），抽取每个脑区的时间序列，以每个脑区作为网络的节点，然后求每两个脑区之间的时间序列相关系数，得到功能网络的连接矩阵。

图2为HCP的多模态分区模板，理解复杂的人类大脑皮层需要一个它的主要子区域（即脑区）的地图（或分割）。 制作准确的地图已经是神经科学的一个世纪以来的目标。HCP使用自己产生的高精度的多模态磁共振图像和客观的半自动神经解剖学方法，把每个脑区划分为180个区域，区域之间的分解线上有着皮质结构，功能，连接，拓扑上的急剧变化。除了之前用解剖显微技术和其他专业性的具体方法得到的83个脑区，还发现了97个新脑区。

4 结果

PLS回归的结果中，有每个基因的与大脑结构相关的权重，依据权重系数对基因排序后，再运用超几何分布检测致病基因的显著性。在超几何分布中，N=16938，是总的基因个数，M是致病基因的总个数，n是排序靠前的若干个基因，k是这n个基因中包含致病基因的个数。P值的阈值设为0.05，小于0.05即表现出显著性。

在对阿兹海默症的研究过程中，取前100个基因，发现了致病基因在VFC，M1C两个脑区的显著性，取前200个基因，没有发现新的脑区，前300个基因，又发现了致病基因在A1C，MFC，ITC三个脑区的显著性，取前400、500个基因，这几个脑区仍然具有显著性，在查閱文献后，这几个脑区都得到了验证，是受阿兹海默症影响的脑区。证明了本研究所用方法的有效性，为致病基因的发现提供了一定程度的参考。在对精神分裂症的研究过程中，发现了其致病基因在8岁～9岁，65岁～75岁两个年龄段的显著性，查阅文献后得到了确认。在对帕金森症的研究过程中，发现了其致病基因在A1C，M1C，IPC，V1C 4个脑区的显著性，查阅文献后得到了确认。

功能脑网络：

功能性连接描述节点（ 可表示神经元、神经集群、功能脑区等不同尺度上的脑功能单元） 之间的功能性信号在某一时段内统计意义上的关系，但不反映节点之间的因果关系。功能性脑网络是基于脑/神经的功能性信号（ 电信号、磁信号、反映血液动力学或代谢的信号等） 构建而成，在微观尺度上可由单个神经元之间的电位发放关系构建，在中尺度（ meso-scale） 上可由反映神经元集群活动的局部场电位来构建，在宏观尺度上，可由特定功能脑区之间的EEG/MEG/fMRI 等来构建。

图4是通过R-fMRI构建人脑功能性脑网络的流程图。（1）提取时间序列数据C. 提取的是rfMRI数据，B.每个解剖单元（即网络节点的数据都要提取）B.每个解剖单元来自自现有的脑区划分模板 A. 每个解剖单元也可以是体素。（2）计算功能连接相关系数矩阵（即网络的边）D. 计算任意两个节点时间序列数据的相关系数。（3）将相关矩阵阈值化为二元连接矩阵（即关联矩阵E）。（4）将关联矩阵可视化为图（F）。

GRETNA是一个运用图论的网络分析工具箱，它整合了目前神经科学领域研究的大多数网络指标，使研究人员能够对脑连接体的拓扑结构进行综合分析。本研究运用该工具计算了每个人脑功能网络的各项指标，包括对相关系数取绝对值后二值化网络，二值化的方法是采用稀疏阈值的方式，即以网络的实际边数与该网络理论上能有的最大边数的比值作为阈值，这种方法能够保证，对节点数相同的网络取阈值后，有相同个数的边被保留下来，以便进行同一阈值下不同网络间的比较。由于对阈值的定义没有明确的标准，阈值从相关系数0.05到0.4，等间隔取36个阈值，得到不同的网络。生成100个有相同节点数，边数以及度分布的随机网络，在GRETNA中，随机网络用马尔科夫链算法生成。脑网络和随机网络进行比较是为了检验其特征与随机网络是否有显著不同。计算的全局特征包括：小世界（small world），全局效率（Global Efficiency），局部效率（Local Efficiency），富人俱乐部系数（Rich Club），模块（Modularity），生成模块的算法为贪婪优化算法，层次性（hierarchy），同步性（synchronization）。节点的局部特征包括：度，节点效率，近中心性（betweenness）。

小世界（small world）：小世界系数定义为给定网络的聚类系数或特征路径长度与有相同度分布的随机网络的聚类系数和特征路径长度的比值。与此随机网络相比，小世界网络有大的聚类系数，小的路径长度。

效率（Efficiency）：指通过网络交换信息的效率。小世界网络是全局和局部都高效的网络。

全局效率（global efficiency）：网络所有节点之间最短路径的倒数之和除以N*（N-1），N是网络的节点个数。

局部效率（local efficiency）：是子网络全局效率的平均值。

富人俱乐部系数（Rich Club）：网络中的富人俱乐部现象是指，网络的主要节点（hub）之间的联系比其他节点之间的联系要密切。富人俱乐部现象的存在能提供一个网络高级结构的重要信息，如弹性，层次性，特异性。

模块（Modularity）：即网络中的社区结构，是网络的顶点被分成组的趋势，组内的连接密集而组间的连接稀疏。

层次性（hierarchy）：真实网络是无标度网络，即节点的度分布是幂律分布。同时又有高的聚类度。小的模块之间以层次结构彼此组合，从而产生了网络的层次性。

度：與每个节点相连的边数。

中介中心性（betweenness）：中介中心性指的是一个结点担任其他两个结点之间最短路的桥梁的次数。一个结点充当“中介”的次数越高，它的中介中心度就越大。中介中心性主要是由美国社会学家林顿·弗里曼（Freeman，1979）教授提出来的一个概念，它测量的是一个点在多大程度上位于图中其他“点对”的“中间”。他认为，如果一个行动者处于多对行动者之间，那么他的度数一般较低，这个相对来说度数比较低的点可能起到重要的“中介”作用，因而处于网络的中心。中介中心性的另一个重要作用就是它能够分辨出谁是“跨界者”（boundary spanners）——那些在两个或多个团体中扮演着不可或缺的桥梁作用的个体。

人类大脑功能网络展示了高效的小世界，分类，分层和模块化组织特性，并拥有高度关联的中心节点，并且这些发现对不同的分析策略是有鲁棒性的。

在构建的静息态功能网络中，发现了脑功能网络的模块化结构，长尾度分布，以及富人俱乐部系数曲线。与有相同度分布的随机网络相比，功能网络的全局特征都表现出了非平凡性，如同配性，聚类系数，全局效率，最短路径长度，小世界特性以及模块化程度。随着年龄增长，这些特性在减弱，趋向于随机网络，如分别取两个阈值时，模块数由第一个年龄段的7个变为第五个年龄段的6个，及由第一个年龄段的3个变为第五个年龄段的2个。这是有序到无序的过程，表现了生物体衰老过程中的熵增。

参考文献：

[1] Glasser M F， Coalson T S， Robinson E C， et al. A multi-modal parcellation of human cerebral cortex[J]. Nature， 2016.

[2] Bellivier F， Golmard J L， Rietschel M， et al. Age at onset in bipolar I affective disorder： further evidence for three subgroups[J]. American Journal of Psychiatry， 2003， 160（5）： 999-1001.

[3] Tzourio-Mazoyer N， Landeau B， Papathanassiou D， et al. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain[J]. Neuroimage， 2002， 15（1）： 273-289.

[4] Salimi-Khorshidi G， Douaud G， Beckmann C F， et al. Automatic denoising of functional MRI data： combining independent component analysis and hierarchical fusion of classifiers[J]. Neuroimage， 2014， 90： 449-468.

[5] Bullmore E， Sporns O. Complex brain networks： graph theoretical analysis of structural and functional systems[J]. Nature Reviews Neuroscience， 2009， 10（3）： 186-198.

[6] Cohen J D， Daw N， Engelhardt B， et al. Computational approaches to fMRI analysis[J]. Nature Neuroscience， 2017， 20（3）： 304-313.

[7] Destrieux C， Fischl B， Dale A， et al. Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature[J]. Neuroimage， 2010， 53（1）： 1-15.

[8] Burgess G C， Kandala S， Nolan D， et al. Evaluation of denoising strategies to address motion-correlated artifacts in resting-state functional magnetic resonance imaging data from the Human Connectome Project[J]. Brain Connectivity， 2016， 6（9）： 669-680.

[9] Gordon E M， Laumann T O， Adeyemo B， et al. Generation and evaluation of a cortical area parcellation from resting-state correlations[J]. Cerebral cortex， 2014： bhu239.

[10] Wang J， Zuo X， He Y. Graph-based network analysis of resting-state functional MRI[J]. Frontiers in systems neuroscience， 2010， 4： 16.

[11] Wang J， Wang X， Xia M， et al. GRETNA： a graph theoretical network analysis toolbox for imaging connectomics[J]. Frontiers in human neuroscience， 2015， 9： 386.

[12] Satterthwaite T D， Baker J T. How can studies of resting-state functional connectivity help us understand psychosis as a disorder of brain development？[J]. Current opinion in neurobiology， 2015， 30： 85-91.

[13] McKeown M J， Hansen L K， Sejnowsk T J. Independent component analysis of functional MRI： what is signal and what is noise？[J]. Current opinion in neurobiology， 2003， 13（5）： 620-629.

[14] Smith S M， Miller K L， Salimi-Khorshidi G， et al. Network modelling methods for FMRI[J]. Neuroimage， 2011， 54（2）： 875-891.

[15] Marrelec G， Krainik A， Duffau H， et al. Partial correlation for functional brain interactivity investigation in functional MRI[J]. Neuroimage， 2006， 32（1）： 228-237.

[16] Abdi H. Partial least squares regression and projection on latent structure regression （PLS Regression）[J]. Wiley Interdisciplinary Reviews： Computational Statistics， 2010， 2（1）： 97-106.

[17] Whitaker K J， Vértes P E， Romero-Garcia R， et al. Adolescence is associated with genomically patterned consolidation of the hubs of the human brain connectome[J]. Proceedings of the National Academy of Sciences， 2016， 113（32）： 9105-9110.

[18] Guerra-Carrillo B， Mackey A P， Bunge S A. Resting-state fMRI： a window into human brain plasticity[J]. The Neuroscientist， 2014， 20（5）： 522-533.

[19] Short-term test–retest reliability of resting state fMRI metrics inchildren with and without attention-deficit/hyperactivity disorder

[20] Kang H J， Kawasawa Y I， Cheng F， et al. Spatio-temporal transcriptome of the human brain[J]. Nature， 2011， 478（7370）： 483-489.

[21] Yi L， Wang J， Jia L， et al. Structural and functional changes in subcortical vascular mild cognitive impairment： a combined voxel-based morphometry and resting-state fMRI study[J]. PloS one， 2012， 7（9）： e44758.

[22] Yao Y， Lu W L， Xu B， et al. The increase of the functional entropy of the human brain with age[J]. Scientific reports， 2013， 3： 2853.

[23] Glasser M F， Sotiropoulos S N， Wilson J A， et al. The minimal preprocessing pipelines for the Human Connectome Project[J]. Neuroimage， 2013， 80： 105-124.

[24] Lindquist M A. The statistical analysis of fMRI data[J]. Statistical Science， 2008： 439-464.

[25] Cao M， Wang J H， Dai Z J， et al. Topological organization of the human brain functional connectome across the lifespan[J]. Developmental cognitive neuroscience， 2014， 7： 76-93.

[26] Murphy K， Fox M D. Towards a consensus regarding global signal regression for resting state functional connectivity MRI[J]. NeuroImage， 2016.

[27] Beckmann C F， Smith S M. Probabilistic independent component analysis for functional magnetic resonance imaging[J]. IEEE transactions on medical imaging， 2004， 23（2）： 137-152.

[28] Beckmann C F， DeLuca M， Devlin J T， et al. Investigations into resting-state connectivity using independent component analysis[J]. Philosophical Transactions of the Royal Society of London B： Biological Sciences， 2005， 360（1457）： 1001-1013.

【通聯编辑：唐一东】