一种基于功能模块的生物分子网络分析的方法

摘要

本发明属于生物信息技术领域。本发明提供了一种基于功能模块的生物分子网络如基因表达调控网络或蛋白质相互作用等网络间比较的方法，主要包括如下步骤：构建生物网络的邻接矩阵Madj，计算网络节点间的功能相似性矩阵Msim，计算网络边的功能权重矩阵，最小图熵算法挖掘网络模块，最后对网络模块进行功能富集分析，其中符号含义如说明书中所定义。

说明书

 

技术领域

本发明属于生物信息技术领域。更具体而言，本发明涉及生物分子网络如蛋白质-蛋白质相互作用网络或基因表达调控网络。

背景技术

过去的几十年，分子生物学的研究对象主要是生物体内部的单个组织、细胞或基因，利用的是经典还原论的思想，研究难免留于局部。实际上，生物体是一个复杂系统，生物分子之间存在层次性关联及相互作用，且生命现象并不是一个静态而是一个动态演化的过程。网络分析由于它的系统性及内部元素间的关联性已成为研究现代生物学及系统生物学的一个重要方法。

随着生物芯片、高通量测序、质谱等分子生物学实验技术的飞速发展，产生了海量的基因组学、转录组学、蛋白质组学及代谢组学等高通量数据，并伴随之产生了常见的生物网络如蛋白质相互作用网络、基因转录调控网络、代谢调控网络、非编码基因调控网络及各种跨组学调控网络等。这些网络的应用极大地促进了对基因表达调控机制、蛋白质功能机制、代谢物传递机制等的研究，并广泛应用于疾病机制研究，药物相互作用研究，育种遗传研究和考古进化研究等领域。

生物分子网络的复杂性，一方面体现在网络的分子组成，比如人类的基因数为20,000~30,000个，人类蛋白质及其修饰的种类更是指数级增长；另一方面是生物分子间关系的复杂性，比如基因表达转录调控、转录后调控、非编码基因调控及蛋白质翻译后修饰等。一般网络的拓扑性质可以从全局进行分析，但从整体上研究生物分子网络的功能可行性较低。此外，生物分子网络还具有无尺度性、小世界特征和模块化特征，即网络中大部分的节点度较小，而少部分节点与其它节点的连接度比较高；网络中的大部分节点都不是直接相连的，但可以通过有限几个节点进行连通；网络中的部分节点呈现高聚类性，模块内部连接紧密，模块间连接比较稀疏。模块是实现生物功能的主体，与功能之间的联系也更为紧密，承担着各个生物学功能的模块被定义为网络的功能模块。因此生物网络分析的策略是：首先从总体上分析网络的拓扑学性质，然后是挖掘全局网络中的模块，最后基于网络模块功能分析研究全局网络的功能机制。目前常用的网络模块识别算法和软件有层次聚类、k近邻、k均值、AutoSOME、MCL、MCODE、GLAY和MINE等，其中MCL马尔科夫聚类算法是一种基于随机流模拟的图聚类算法，MCODE是给定网络中每个节点一个权值(可以为节点的连接度)，根据贪婪算法分离出网络中稠密的区域，再对这些区域进行整合，添加或删除节点。上述的方法中，无一例外均是采用网络拓扑属性进行网络模块的识别，然后基于基因本体论(Gene Ontology, GO)方法进行功能富集得到模块的功能。这些方法存在的共同问题是在挖掘网络模块的时候，依据的是网络中节点的拓扑学性质（度、聚类系数、介数等）而没有考虑相邻节点间的功能相似性。

发明内容

本发明结合了网络中节点间的功能相似性和拓扑性质，使得生物分子网络由无权重的网络成为边具有功能相似性权重的网络，然后再基于权重模块识别方法来挖掘生物网络中的功能模块。

本发明提供了一种基于功能模块的生物分子网络分析方法，通过网络中节点间的功能相似性结合网络的拓扑性质，筛选出网络中功能关联的一些节点和边，组成功能模块。

在一个实施方案中，对于数学化的生物网络G(V，E)，V表示网络中的节点集合，E表示网络中边的集合，本发明是这样实现的，主要包括如下流程：

步骤1，基于所述生物网络G，构建网络节点间的邻接矩阵M_adj，矩阵的行和列分别表示网络中的节点集合V，矩阵由0和1构成，如果两个节点间存在边关系，则矩阵中此两个节点所在行和列相交的元素为1，其他情况则为0。

步骤2，构建网络节点间的功能相似性矩阵M_sim，矩阵的行和列分别表示网络中的节点集合V，矩阵由0.001-1之间的小数构成，如果两个节点间存在语义相似性，则此矩阵中此两个节点所在行和列相交的元素为这两个节点的功能语义相似性得分；功能语义相似性得分计算优选基因本体论中的分子功能(Molecular Function)或生物学过程(Biological Process)语义，其次是细胞学组成(Cellular Component)语义。如果两个分子节点间存在完全语义相似性则值为1，如果两个分子节点间不存在语义相似性则值为0.001，功能语义相似性方法优选GOSemSim算法(Yu G, Bioinformatics 2010)。

步骤3，计算网络边的功能权重矩阵M_E，                                                ，“”表示两个矩阵相同位置元素相乘。

步骤4，基于最小图熵算法（图熵，，，Ω表示子模块内的顶点集合，i属于Ω，n_i为节点与模块内的邻接节点间功能权重累加和，N_i为节点在全局网络中与其邻接节点间功能权重累加和）挖掘网络的子模块：首先选定种子节点（例如度不小于5的节点）及其邻居节点作为种子簇，然后逐个地删除种子簇中种子节点的邻居节点，尽可能地减小簇的图熵；接下来，再将种子簇的邻居节点中可以降低簇的图熵的节点添加进来，重复这个过程，直到图熵无法进一步降低，当图熵小于经验阈值（例如20）时循环停止。

在本发明的方法的一个实施方案中，所述方法还包括步骤5，对步骤4识别出的功能模块进行功能富集分析（例如在线功能富集分析工具DAVID或Perl平台的GO::TermFinder工具）。

基因本体论（Gene Ontology，简写为GO）包含了基因参与的生物过程、所处的细胞组成和发挥的分子功能三方面的功能信息，并将这些功能按概念范围的大小组织成“有向无环图”层次结构。两个分子节点间存在完全语义相似性是指着两个分子被注释到相同的Gene Ontology功能节点；两个分子节点间不存在语义相似性是指两个分子注释到的Gene Ontology功能节点在层次结构中的距离大于等于8，两个分子间的语义相似性定义为0.001。语义相似性计算可以通过现有技术中的一些算法进行计算，如R统计软件中的GOSemSim和DOSE插件等。

本发明的优势在于使用网络中邻近节点间的功能相似性作为节点的一个重要属性，使得模式识别算法挖掘到的功能模块内部具有较高的功能关联性，和传统的仅基于拓扑性质的模块发现方法相比更能分析出有真实生物学意义的功能模块。

本发明的创新点在于：

1) 将网络中的节点间相似性得分矩阵与节点邻接矩阵有效的整合，使得生物分子网络成为一个边具有功能权重的网络；

2) 最小图熵算法的改进，以往算法是基于节点的连接度进行模块识别，再对模块进行功能分析，而本发明的方法则是基于网络中节点间功能相似性直接进行功能模块的识别。

附图说明

图1示出本发明所述生物网络功能模块识别方法的流程图；

图2示出小鼠蛋白质相互作用网络；

图3示出小鼠蛋白质相互作用网络挖掘到的功能模块之一；

图4示出小鼠蛋白质相互作用网络挖掘到的功能模块之一；

图5示出本发明方法与MCODE、MCL和AutoSOME挖掘到的模块的newman模块度比较；

图6示出一种药物硝苯地平治疗机制相关的功能模块；

图7示出一种药物硝苯地平治疗机制相关的功能模块基因在钙信号通路的富集；

图8示出人类蛋白质相互作用网络中挖掘到的2个功能模块。

具体实施方式

AutoSOME：一种基于距离的网络聚类算法。

MCL：马尔可夫聚类算法。

GLAY：团体聚类(Community clustering)方法。

MINE：一种类似于MCODE改进的聚类算法。

MCODE：分子复合物检测(Molecule Complex Detection)算法。

Gene Ontology(GO)：基因本体论，基因功能注释数据库。

DAVID：注释、可视化和整合分析数据库。

GO::TermFinder：一个基于Perl语言的Gene Ontology富集分析工具。

Hub节点：与网络其它节点连接较多的节点被定义为Hub节点。

BIOGRID： 蛋白互作和遗传互作数据库。

INTACT：欧洲分子生物实验室分子相互作用数据库。

MINT：分子相互作用数据库(Molecular INTeraction database)。

Entrez Gene ID：由美国国立生物信息中心NCBI提供的一种基因的标识符。

Cytoscape：一种开源的生物分子网络可视化工具平台。

ClusterMaker：一个基于Cytoscape平台开发的整合了多个聚类算法的网络模块识别工具。

本发明的目的是从复杂的生物网络中去除不真实的网络关系，保留最可靠的网络关系，这些被保留的网络关系代表了真实的生物学意义，从而能从系统层面提供生命现象的线索。以下的实施例证明了本发明方法的有效性。这些实施例是非限制性的，本发明的方法可以用于其他生物网络。

实施例一、以高血压小鼠模型硝苯地平给药治疗后的一组基因芯片表达数据构成的基因和蛋白质相互作用网络为例，具体地描述本发明方法：

整合BIOGRID、INTACT、MINT、NIA Mouse Protein-Protein Interaction Database等蛋白质相互作用数据库中的蛋白质相互作用数据，在剔除重复数据和自相互作用之后，得到了一个全局的小鼠基因和蛋白质相互作用网络。由于这类数据比较有限，因此再利用小鼠蛋白在其它模式生物中的直系同源和旁系同源蛋白质的相互作用数据预测出这些蛋白质在小鼠中的相互作用，共收集了65,850条小鼠蛋白质相互作用的数据。

通过对基因芯片数据进行均值中心标准化(mean centered normalization)，将标准化后表达值大于1的基因，视为实验中有表达的基因。将这部分基因映射到上述小鼠蛋白质相互作用网络上，最终得到一个对应的蛋白质相互作用网络（定义为G）。网络G如图2所示，共包含2,183个基因及6,656条边。

步骤1：计算网络G的邻接矩阵M_adj，M_adj为一个2,183行2,183列的矩阵，每行每列各表示唯一的一个基因，如果邻接矩阵6,656条边中存在某两个基因的互作关系则矩阵中相应元素为1，否则为0，M_adj的部分元素如下：

步骤2：采用R的GOSemSim软件包，基于Gene Ontology的分子参与的生物学过程，计算2,183个基因间两两间的语义相似性得分，M_sim的部分元素如下：

Entrez Gene ID74596110911663996865374596 Cds1110.4380.089110911 Cds2110.4380.08966399 Tsfm0.4380.43810.07968653 Samm500.0890.0890.0791

 步骤3，计算网络边的功能权重矩阵M_E，，“”表示两个矩阵相同位置元素相乘；计算得到M_E的部分元素相应如下：

Entrez Gene ID74596110911663996865374596 Cds1000.4380.089110911 Cds2000.4380.08966399 Tsfm0.4380.43800.07968653 Samm500.0890.0890.0790

步骤4：基于最小图熵算法（图熵，，，Ω表示子模块内的顶点集合，i属于Ω，n_i为节点与模块内的邻接节点间功能权重累加和，N_i为节点在大网络中与其邻接节点间功能权重累加和）第3步得到的权重网络共挖掘得到26个功能子模块，图3和图4是挖掘的功能模块的示例，从图中可以明显看出本发明方法挖掘到的功能模块内部分子具有较高的功能一致性，比如图3中POLR2G、POLR2J、POLR2F均属于RNA聚合酶二家族，YWHAH和YWHAG属于酪氨酸/色氨酸激活蛋白；图4中存在着大量琥珀酸脱氢酶家族成员蛋白(Sdh)，线粒体核糖体蛋白(MRP)。以图3为例，假设图3就是一个大网络，POLR2F、POLR2G、POLR2J、TCEA1和MED6为一个种子模块，则POLR2F节点的n值为POLR2F与POLR2G功能权重1、POLR2F与POLR2J功能权重1、POLR2F与TCEA1功能权重0.90及POLR2F与MED6功能权重0.54之和3.44，N为所有与POLR2F连接基因间的功能权重累加和4.44，则p为0.775(3.44/4.44)。

步骤5：基于在线功能富集分析工具DAVID对步骤4挖掘的26个功能模块分别进行功能分析，以p<0.1为富集分析显著性阈值，26个模块显著富集到的生物学过程数目最少3个，最多16个，模块内节点富集的功能具有紧密的关联。比如图3所示的模块富集分析结果如下：

GOIDTermCount%p.valueGO:0006350基因转录947.3684210.0010689GO:0046907细胞内转运526.3157890.0021446GO:0008104蛋白定位631.5789470.0024133GO:0006886细胞内蛋白转运421.0526320.0053959GO:0006913细胞核与细胞质中蛋白转运315.7894740.0070212GO:0051169细胞核转运315.7894740.0073069GO:0032774RNA合成315.7894740.0118443GO:0006605蛋白靶向315.7894740.0131305GO:0045449调控转录842.1052630.0188419GO:0006357RNA聚类酶II结合启动子的调控421.0526320.0455435GO:0015031蛋白转运421.0526320.0522363GO:0051170细胞核进口210.5263160.074152GO:0034504细胞核内蛋白定位210.5263160.0778426

图3所示功能模块显著富集分子功能主要有蛋白转运、蛋白定位、转录、调控转录等，这些生物学过程全部与RNA转录密切相关。

对于初始的硝苯地平给药后小鼠蛋白质相互作用网路G，用Cytoscape的ClusterMaker工具包选用AutoSOME、MCL和MCODE三种算法分别挖掘网络G中的模块，AutoSOME、MCL和MCODE挖掘到包含5个蛋白节点以上的模块分别为182个、112个和49个。采用Newman模块度（Newman M．E．J．and M．Girvan．Finding and evaluating community structure in networks，模块度是一个复杂网络模块性的度量指标，它不仅考虑了模块的内聚性，也考虑模块之间的耦合性）评估AutoSOME、MCl、MCODE和本发明方法，如图5所示，本发明方法挖掘到模块的模块度（中值0.0460）与MCODE（中值0.052）相近较高，明显高于MCL（中值0.0169）和AutoSOME（中值0.00515）算法，说明本发明方法挖掘到的模块内部紧密程度与MCODE相当，且都具有较强的内部关联。

硝苯地平通过1）阻止钙内流2）阻碍细胞壁内钙的释放；3）阻断血管膜上α-肾上腺素受体4）抑制磷酸二酯酶的活性5）与钙调节素相作用6）激活Na+，K+ -ATP酶7）激活钙离子泵等几种方式实现抑制血管、支气管和子宫平滑肌的兴奋-收缩偶联，扩张全身血管，达到降低血压的效果。通过DAVID通路富集分析发现图6所示功能模块中的基因显著富集在MAPK信号通路、钙信号通路、心肌收缩及平滑肌收缩等通路，与硝苯地平的部分药理机制相近相符。图7为图6中功能模块基因在钙信号通路的富集。

通过小鼠蛋白质互作网络的分析实例，可以看出本发明方法可以有效地挖掘网络中具有密切功能关联的子模块。

实施例二、以人类Human Protein Reference Database (HPRD)数据库中提供的实验验证的39,240个蛋白质相互作用为基础，从中随机筛选3000条蛋白质相互作用关系，去除蛋白质自身相互作用，得到由2,095条边1,478个蛋白节点的网络G。

步骤1：计算网络G的邻接矩阵M_adj，M_adj为一个1,478行1,478列的矩阵，每行每列各表示唯一的一个基因，如果随机网络G中存在某两个蛋白的互作边则矩阵中相应元素为1，否则为0。

步骤2：采用R的GOSemSim软件包，基于基因本体(GO)中的生物学过程，计算1,478个蛋白质两两间的语义相似性得分M_sim，M_sim的中值为0.307。

步骤3，计算网络G边的功能权重矩阵M_E，，“”表示两个矩阵相同位置元素相乘，M_E矩阵的中值为0.289。

步骤4：基于最小图熵算法，人类蛋白质互作网络G挖掘得到20个功能子模块（模块内包含的蛋白节点大于10）。图8为2个挖掘到的功能子模块示例，图5左侧模块中可以直观看出ORC6L和ORC2L均为DNA复制起始相关的复合物，SMARCA4和SMARCB1均属于一种染色质调节蛋白；右侧模块中可以直观看出GTF2H1和GTF2H5均属于转录因子II家族蛋白，ACTL6A和ACTR6则分别是肌动蛋白及其受体。

步骤5：基于在线功能富集分析工具DAVID对步骤4挖掘的20个功能模块分别进行功能分析，以p<0.1为富集分析显著性阈值，其中图8右侧模块显著富集的生物学过程为组蛋白修饰、组蛋白乙酰化、转录、RNA生物合成、DNA修复和DNA代谢等紧密关联的功能。

通过上述分析实例，可以得出虽然本发明仅仅是针对全部人类蛋白质相互作用网络的一部分进行模块挖掘，但加入互作节点间的功能相似性参数，依然可以有效地挖掘出功能密切关联的模块。

以上是对本发明的描述，而非限定，基于本发明的思想的其他实施方式，均在本发明的保护之中。