基于生物分子网络分析的精神分裂症功能模块挖掘

廖苑君,陈应坚,赵小蕾,孙胜南,林 帆,覃继恒,饶绍奇*

(广东医科大学a.公共卫生学院;b.医学系统生物学研究所,中国广东东莞523808)

精神分裂症(schizophrenia,SCZ)是一种严重影响人类生命健康和生存质量的复杂性精神疾病,患病率约为1%,遗传贡献高达83%[1]。患者可出现幻觉、妄想、情感淡漠、社交退缩等典型精神症状,同时自知力丧失或不完整,造成严重的思维、情感、行为等多方面的障碍[2]。SCZ病症具有较高的复发率和精神症状残留,且具有不可预知性,给家庭和社会带来巨大的经济负担,并占据了大量的医疗资源,导致资源利用价值下降。在此背景下,研究SCZ病因与病理机制对其预防及治疗具有重要的意义。目前,大规模组学数据的发展为研究遗传学机制提供了便利,然而分析方法主要采用的是单位点、单基因分析策略[3]。由于SCZ病因复杂,其发生、发展过程涉及多个生物分子(如基因、miRNA等),这些生物分子构成了一个复杂的网络,因此与单位点、单基因分析方法相比,生物分子网络分析能够更加全面揭示SCZ的分子机制[4]。

在生物分子网络分析中,节点代表生命体内的各种分子,边(网络连接)代表分子间的互作关系。根据对象的不同,生物分子网络可分为基因调控网络、蛋白质互作网络、信号转导网络等[5]。迄今为止,生物分子网络分析在生物医学领域已被广泛应用,具有成熟的分析方法,可用来筛选和挖掘网络中关键的核心节点及功能模块,从而提示具有重要生物学意义的生物联系或功能结构。Xu等[6]发现基于功能注释定义的功能模块不仅可以用于癌症的分型或预测,甚至能够挖掘癌症中潜在的亚型。Zhao等[4]通过开展冠心病风险通路内、通路间网络研究,挖掘了冠心病与阿尔茨海默病共享的功能模块。本研究借助网络分析理论,深入挖掘与SCZ相关的功能模块,以揭示潜在的分子致病机制,为SCZ的基础研究提供新思路。

1 材料与方法

1.1 数据来源

作为SCZ遗传公共数据库之一,SZDB数据库[7](http：//www.szdb.org/)收录了大量的数据集和候选基因。对SZDB数据库进行挖掘,纳入标准：1)基因表达数据集;2)数据集的患者病例数超过10例,经筛选得到5套完整数据集[8～12]。利用该数据库提供的在线处理分析工具筛选差异表达基因。首先对数据进行标准化(归一化)预处理,随后利用R和Bioconductor平台中的Limma包进行差异表达分析。考虑到基因表达丰度低的mRNA的差异表达倍数不一定足够大,本研究不设置差异倍数阈值。为了更好地控制差异表达结果和假阳性率,将在3套或以上数据集中具有显著性差异表达(P＜0.05)的基因纳入后续分析。最终共挖掘了944个与SCZ显著相关的基因,并定义为初始种子基因集。

HPRD 数据库[13](http：//www.hprd.org/)是人类蛋白质参考数据库,含有人类蛋白质互作注释信息。迄今为止,该数据库储存了41 327对蛋白质互作(protein-protein interaction,PPI)对子数据,数据可信度高,证据可靠。本研究提取该数据库的PPI对子数据,构建致病基因网络。

1.2 致病基因网络的构建

本研究构建的网络节点集合包含了初始种子基因集及其在PPI的一级邻居基因,对于集合中的任意两个不同节点,假设其产物在PPI网络中存在互作关系即视为网络中连接。因此,构建的网络可视为SCZ特异性的风险基因网络,其中节点代表基因,节点间的边代表基因间相互作用。

1.3 网络功能模块的挖掘

生物分子倾向于协同执行一定的生物过程或发挥特定生物功能。同样地,在生物分子网络中,其对应的节点连接高度紧密,聚合为节点簇,构成网络的模块化(modularity)。因此,与同样规模的随机网络相比,模块化网络连接更为紧密。从功能模块的角度分析网络,能够将生物分子进行功能分区的划归,对预测其功能、理解生命活动的层级和结构具有关键提示作用。功能模块的结构划分评价可看作是网络结构分解的优化问题。基于研究假设,本研究通过Newman算法[14～15]对构建的网络进行分解,从而识别网络中高度模块化的功能模块。该算法将网络分解问题转变为二次型优化问题,并将之简化为特征向量提取问题。Newman算法拟优化的目标函数为其中,当Aij=1时说明节点i和节点j之间存在互作关系,反之取值为0;m代表网络中边的总数;ki和kj分别代表节点i和j的连通度;s代表有n个节点指示变量的列向量,向量内元素的取值为1或-1,si和sj分别代表节点i和j所属的功能模块,即si=即两节点位于同一功能模块时,Q才能取得最大值。为了方便,将等式转换为二次型的形式：Q=为s的倒置,由列向量变为行向量;B代表着一个内部元素为Bij=Aij-kikj/2m的对称矩阵,当s取B的最大特征值对应的特征向量时,Q值最大。因为s只可取值为±1,所以s只保留矩阵特征向量对应的正负号,通过符号方向将网络分解为两部分。通过以上方法重复网络一分为二的过程,直到网络不能分解为止。最终将这些稠密、不可分割的子网提取出来,并定义为功能模块。本研究设定的功能模块节点数大于50,上述算法通过R语言“igraph”软件包完成。

1.4 网络拓扑属性分析

对网络的拓扑属性进行评价时,评价指标包括以下五方面的内容。

1)网络直径(network diameter)：即网络中任意两个连通节点间最短路径的最大值,可反映网络的大小。

2)连通度(degree)：某个节点的连通度是指网络中直接与该节点相连的边数,一般用符号k表示,能够衡量该节点在维持网络完整性中的作用。在随机网络中,连通度的分布符合泊松分布[16],通过检验某个节点的连通度是否显著大于一般节点,可进一步筛选和挖掘核心节点。基于理论假设,某个节点的连通度大于或等于t的概率为：其中 λ=NP,λ 代表节点的期望连通度,N代表节点数目,P代表任意两个节点发生连接的概率。通过Bonferroni校正,以P＜0.01为差异有统计学意义。当P值小于预设检验水准时,该节点即定义为核心节点。

3)介数(betweenness)：描述网络中某节点的路径占其他节点间最短路径的比例。介数越高,表示其在维持网络紧密连接性方面起更加重要的作用。

4)聚类系数(clustering coefficient)：衡量网络中某节点与周围节点的密集连接程度。

5)无标度网络(scale-free network)：指由少数高连通度节点及多数低连通度节点组成的网络。一般认为,无标度网络中连通度的分布符合幂律分布,即P(k)∝k-α,其中k代表接通度,α代表指数参数,使用极大似然估计法[17]。本研究采用KS(K-olmogorov-Smirnov)拟合优度检验评判无标度特性。

1.5 功能富集分析

通过R语言(版本3.5.0)“clusterProfiler”软件包[18]对每个功能模块进行KEGG(Kyoto Encyclopedia of Genes and Genomes)功能富集分析。通过FDR校正,以P＜0.001为差异有统计学意义,并针对每个功能模块最显著的前5条富集通路开展后续分析。根据富集到的KEGG通路类别评估功能模块之间的相互作用,利用Cytoscape软件(版本3.6.1)构建功能模块与通路分类的关系。

本研究遵循上述流程开展,技术路线见图1。

2 结果

2.1 网络的构建与拓扑属性

通过PPI引导,构建的SCZ特异性基因网络含有3 207个节点和5 207个互作对子。由于网络中的节点之间并非全部紧密连接,故将孤立作用的节点(TRIM63、PDK4、DLAT 等)去除掉,只对最大子网进行后续分析,最终生成的网络包含3 062个节点和5 120个互作对子。通过网络拓扑属性分析得到网络直径为12,聚类系数为0.034,网络平均路径长度为4.972。

2.2 网络功能模块的挖掘与拓扑属性分析

通过Newman分解算法对最终生成的网络进行分解,得到14个节点数大于50的功能模块。14个功能模块的拓扑属性和无标度检测结果如表1所示。从表1可知,功能模块内的节点数差别很大。例如：最大的功能模块(M9)由522个节点和860条边组成,每个节点平均含有1.6条边,而最小的功能模块(M14)由51个节点和52条边组成,每个节点平均含有1条边。从直径的分布结果来看,功能模块尺寸越大,其直径反而呈下降的趋势,这表明随着功能模块尺寸的增加,节点间聚集的程度会越发紧密。这些功能模块的拓扑性质指标接近,直径数值在8～15波动,特征路径长度均小于6,呈现“小世界”的特性,符合六度分离理论,支持了生物学上无标度网络的常见表现[19]。无标度检验结果显示,各功能模块指数参数估计值均为2～3,连通度均近似符合幂律分布(KS检验,P＞0.05),支持了功能模块的无标度性质,可视为无标度网络。

图1 数据获取及分析流程图Fig.1 Flow chart of data acquisition and analysis

2.3 核心基因的识别

利用泊松分布检验筛选出102个核心基因(P＜0.01)。通过检索PubMed数据库,得到57个核心基因在已有文献中被报道与SCZ相关,其余的45个核心基因暂无文献支持其与该疾病相关(表2)。在14个功能模块中,M9含有的核心基因最多,多达21个;M7含有的核心基因最少,仅有1个。图2是我们利用R语言软件展示的部分功能模块与核心基因分布图。

表1 各功能模块的基本特征和拓扑属性Table 1 Basic characteristics and topological properties of each functional module

根据连通度对基因进行排序,列出连通度较大的前10个基因的拓扑参数(表3),结果显示EGFR基因具有最高的连通度,维系着网络中15.29%的基因间最短路径的连通路径,提示其在网络中发挥重要作用。

2.4 模块的功能富集分析

KEGG通路富集分析结果显示,模块内的基因显著富集于细胞凋亡(hsa04210)、ErbB信号通路(hsa04012)、细胞周期(hsa04110)、磷脂酶D信号通路(hsa04072)、PI3K-Akt信号通路(hsa04151)等多个有统计学意义的信号通路(表4)。根据富集到的KEGG通路类别评估这14个功能模块之间的相互作用,利用Cytoscape软件构建功能模块与通路分类的关系,结果如图3所示。图3信息表明,提取的SCZ功能模块主要参与折叠、分类和降解(folding,sorting and degradation)、细胞通讯(cellular communication)、免 疫 系 统(immune system)、传染病(infectious diseases)、癌症(cancers)及内分泌系统(endocrine system)等多个功能类别。

3 讨论

表2 每个功能模块的核心基因Table 2 Hub genes of each functional module

表3 连通度排名前10的网络节点的拓扑属性Table 3 Topological properties of some network nodes

近年来随着基因组学和蛋白质组学的发展,众多研究机构利用组学实验手段产生了大量与疾病相关的实验数据。利用公共数据库资源或整合实验数据,能够帮助研究者抽取有重要生物学意义的生物信息,重构不同层次的生物分子网络。从功能模块的角度分析网络的功能学特征,能够挖掘有意义的生物分子或功能结构。为此,本研究基于网络模型,利用SZDB和HPRD公共数据库构建SCZ特异性的风险基因网络,并进一步通过Newman算法对构建的网络进行分解,挖掘了紧密联系的功能模块和核心基因,分析了功能模块行使的生物学功能。

图2 部分功能模块与核心基因分布图含有标签的节点代表对应功能模块的核心节点,节点的大小代表该基因的连通度大小。Fig.2 Part of functional modules and hub gene distributionThe node containing the tag represents the hub node of the corresponding functional module,and the size of the node represents the degree of the gene.

本研究共挖掘了14个与SCZ相关的功能模块,并通过拓扑学指标评价各模块的拓扑性质,结果显示所有功能模块的拓扑属性接近,符合六度分离理论,具有小世界性,支持了生物学上无标度网络的常见表现。KS拟合优度检验结果显示,所得功能模块的连通度均近似服从幂律分布,符合无标度性。无标度性是指功能模块包含了极少数连通度较高的核心基因,它们位于功能模块的中枢,发挥“信息枢纽”的作用,能够维持其功能的稳定性,而这些核心基因一旦受到破坏将会影响功能模块的功能,进而导致疾病的发生发展。利用泊松分布检验筛选得到102个核心基因,通过检索PubMed数据库,发现共有57个核心基因已被报道与SCZ相关,例如EGFR、HAX1、IL1-R1、RALGDS、NSF等基因。其中,EGFR基因在网络中具有最高的连通度,维系着网络中15.29%的基因间最短路径的连通路径。EGFR是具有酪氨酸激酶活性的表皮生长因子受体家族(epithelial growth factor receptor family)成员。相关研究发现,EGFR能够抑制多巴胺D3样受体的信号传导,是SCZ潜在的治疗靶点[20]。除此之外,一些尚未有报道证实与SCZ相关的核心基因也在研究过程中被发现,例如 SVIL、DNAJA1、RABAC1、STX6 等基因,虽然这些基因仍未经过完善的研究和验证,但从本研究的结果来看,它们可能会成为治疗SCZ的重要靶标。

表4 各功能模块的KEGG通路分析结果Table 4 KEGG pathways of each functional module

通过对所得功能模块进行KEGG信号通路的富集分析,发现大部分通路已有文献支持与SCZ的发病机制相关,例如：TGF-β1(hsa04350)和Hippo信号通路(hsa04390)在诱导Th17细胞的分化(hsa04659)过程中起着积极作用,进而介导组织的炎症反应[21～23];血管紧张素转换酶是肾素-血管紧张素系统(hsa04614)的关键酶,可以催化神经肽的降解并调节多巴胺能和5-羟色胺能神经传递[24];长时程抑制(hsa04730)作为突触可塑性的主要形式之一,在学习和记忆形成的机制中发挥重要作用[25];FcγR介导的吞噬(hsa04666)过程能够诱导肌动蛋白细胞骨架的组装(hsa04810),从而影响细胞的运动、分裂、黏附及吞噬过程,目前研究发现精神疾病的发病与神经突触被过度吞噬密切相关[26];肿瘤抑制因子PTEN对PI3K-Akt信号通路(hsa04151)起到负性调节作用,导致SCZ患者的认知功能受到影响[27～28]。同时,我们也得到了未明确证实与SCZ发病机制相关的通路,例如Ras信号通路(hsa04014)。据报道,Ras基因的失调能够导致5-羟色胺系统和多巴胺系统受损[29];Ras信号通路(hsa04014)在神经细胞的生长、分化等多个生物学过程中具有重要作用[25],推测其与精神疾病的发病密切相关。

本研究根据富集到的KEGG通路类别探寻这14个功能模块之间的相互作用,得到了更为深入的功能联系。结果显示,得到的功能模块主要参与折叠、分类和降解、细胞通讯、免疫、传染病、癌症及内分泌系统等多个功能类别。其中,信号转导(M2,M3,M6,M9,M12)、细胞通讯(M1,M9,M11,M12)以及折叠、分类和降解(M2,M6,M13,M14)这3种功能类别的连通度较高,提示这些功能与SCZ的病理机制关系更为密切。除此之外,功能模块M4、M7、M11共同参与免疫炎症反应;功能模块M3、M5、M7共同调控细胞周期;功能模块M10主要富集于冠心病通路,提示功能模块M10可能成为解释SCZ和冠心病共享的关键功能模块。功能模块分析显示,大部分功能模块不是单独发挥作用的,它们共同影响SCZ的发生发展,具有共享的病理机制。

图3 功能模块与通路分类的关系绿色圆圈表示功能模块,红色圆圈表示生物学通路对应的功能类别。Fig.3 Relationship between functional modules and pathway classificationThe green circle represents the functional modules,and the red circle represents the functional category corresponding to the biological pathway.

总的来讲,本研究从功能模块的角度分析了SCZ网络的功能学特征,为疾病的基础研究提供了新的思路。此外,研究过程采用的方法具有普适性和推广性,能够进一步推广到其他疾病研究。需要指出的是,研究的预测结果在一定程度上受到数据库的信息量完整性和正确性的影响,存在一定的局限性,在今后的研究中仍需要整合更多的数据信息来加以完善。