CN107832583B - 一种基于图匹配的跨物种生物通路发现方法 - Google Patents

Abstract

本发明公开了一种基于图匹配的跨物种生物通路发现方法，本发明是为了解决采用传统生物化学实验方法发现生物通路的低效率问题，和现有图匹配算法无法很好结合生物序列相似性和蛋白质交互网络结构相似性问题。通过本发明可以将生物序列相似性和网络结构相似性很好地融合，能够发现不同物种的蛋白质交互网络中共存的较大子结构，从而更有效地发现存在于不同物种中具有相似功能的生物通路，对生物学研究不同物种之间的联系有指导意义。

Description

一种基于图匹配的跨物种生物通路发现方法

技术领域

本发明属于图算法技术领域，涉及一种基于图匹配的跨物种生物通路发现方法，尤其涉及一种跨物种生物蛋白质交互网络中生物通路的发现方法。

背景技术

图是计算机科学中常用的一类抽象数据结构，以描述事物之间的复杂关系。图结构已广泛应用于多种领域，如万维网、公路网、社交网络、知识图谱、蛋白质交互网络等。随着以上应用领域的发展，图数据不断的快速产生和积累，如何对其进行有效的管理、查询和挖掘等已成为学术界和工业界的研究热点。其中一个重要的课题就是挖掘不同物种生物蛋白质交互网络中的具有相似功能的生物通路。生物通路(biological pathway)由蛋白质和蛋白质之间的交互组成，可视为执行特定生物功能的最小的生物系统，不同物种之间存在大量相似功能的生物通路。有效挖掘发现不同物种之间具有相似功能的生物通路，可以帮助人类在生物体系层次上深刻理解物种间的相似及差别，对于基因学和医药学的发展具有重要的指导意义。

传统生物通路发现方法需要大量的生物化学实验支撑，比较低效。通过图匹配方法，可以发现不同物种生物蛋白质交互网络(PPI网络)中相似结构和功能的子结构，根据这些匹配的子结构再根据生物化学方法验证它们是否是真正具有相似功能的生物通路，就比较有针对性和高效。图匹配问题目标是从两个由节点和边构成的图中得到节点一一映射的关系，实质上是图论中的子图同构问题，本身就是一个无法在多项式时间内解决的问题，随着生物蛋白质交互网络规模的扩大，蛋白质网络匹配问题面临着更加艰巨的挑战。

早期蛋白质网络匹配技术主要采用序列匹配技术，因为蛋白质是由基因序列转码的mRNA转译合成的，所以蛋白质也携带有基因序列信息，通过利用BLAST等序列匹配算法可以将序列信息相近的蛋白质进行匹配，然而只依赖蛋白质序列信息而忽略蛋白质网络的结构特性，导致匹配的准确性不高。后来PathBLAST、MaWISh、Graemlin等启发式算法开始采用计算局部网络相似性来进行蛋白质网络匹配，但这些局部匹配算法可能会造成误导，因为一个物种的蛋白质网络子结构可能会匹配到另一个物种的蛋白质网络中的多个子结构，这种一对多的关系给确定生物通路带来了困难。现在比较流行的是采用全局网络匹配算法，比如IsoRank、PATH、GA、GRAAL、L-GRAAL、NATALIE、GHOST、NETAL、MAGNA、SPINAL、HubAlign等。全局网络匹配算法强调两个蛋白质网络匹配的一对一映射关系，当且仅当两个蛋白质节点的邻居节点也匹配的情况下才会将两个蛋白质节点进行匹配。但现有全局网络匹配算法存在一个弊端，即进行网络匹配时不能很好地将蛋白质序列信息和蛋白质网络结构信息结合起来，造成序列信息匹配度和结构信息匹配度两者的对立，导致比较差的生物通路发现效果。

发明内容

为了解决上述技术问题，本发明提供了一种基于新的图匹配算法的跨物种生物通路发现方法。

本发明所采用的技术方案是：一种基于图匹配的跨物种生物通路发现方法，其特征在于，包括以下步骤：

步骤1：构建初始匹配(matching construction)阶段，通过结合各种相似性度量得到一个初始匹配方案，具体是：首先将两个物种的蛋白质交互网络(PPI网络)G₁和G₂中具有最高全局相似性、局部相似性、序列相似性和度相似性以及度数大于一定阈值的的节点进行匹配作为锚点A，然后从这些锚点A扩展，根据局部相似性和序列相似性匹配锚点的邻居节点，直至两个物种中节点总数较少的PPI网络中的所有节点均已匹配，得到初始匹配M。

具体实现包括以下子步骤；

步骤1.1：采用谱方法计算全局相似性S_g；

对于图G，它的邻接矩阵为A，对角度矩阵为D，拉普拉斯矩阵L＝D-A；对于两个物种的蛋白质交互网络G₁和G₂，分别计算它们的拉普拉斯矩阵的特征值，假设图G₁的拉普拉斯矩阵L₁的特征值为α₁≥α₂≥…≥α_n，图G₂的拉普拉斯矩阵L₂的特征值为β₁≥β₂≥…≥β_n，令Λ₁＝diag(α_i)，Λ₂＝diag(β_i)，L1和L2是对称半正定矩阵，则

其中U₁和U₂是正交矩阵；如果G₁和G₂是同构的，则存在一个列矩阵P使得PL₁P^T＝L₂，解得

则全局相似性

步骤1.2：计算局部相似性S_l；

假设图G中的顶点v的k步邻居子图为N_k(v)，并且

表示包括节点v的完整k步子图，V₁表示图G₁的顶点集合，V₂表示图G₂的顶点集合，节点u∈V₁和节点v∈V₂两者之间的局部相似性通过比较u和v的k步邻居子图来衡量，具体如下：假设d(u)、d(v)分别为节点u和节点v在G₁、G₂中的度，假设N_k(u)的所有节点的度按大小降序排列后分别为d_1,1,d_1,2,…，Nk(v)的所有节点的度按大小降序排列后分别为d_2,1,d_2,2,…；令节点u和节点v的k步子图的较小节点总数n_min＝min{|N_k(u)|,|N_k(v)|}，则G₁的节点u和G₂的节点v之间的局部相似性为

其中

和

分别表示包含节点u的k步子图的顶点数和边数，

和

分别表示包含节点v的k步子图的顶点数和边数，其中k步子图的最小度之和

步骤1.3：计算度相似性S_d和序列相似性S_seq；

假设G1的节点u的度为d(u)，G2的节点v的度为d(v)，那么u和v之间的度相似性为

通过BLAST计算出序列分数seq(u,v)，然后进行归一化得到序列相似性

步骤1.4：融合各种相似性，选择出锚点A；

首先由全局相似性矩阵S_g、局部相似性矩阵S_l得到拓扑相似性矩阵S_t＝S_g×S_l，再结合度相似性S_d得到网络结构相似性S_str(u,v)＝(1-θ)×S_t(u,v)+θ×S_d(u,v)，其中θ表示拓扑相似性和度相似性的平衡参数，值在[0,1]之间可调节。再结合结构相似性S_str和序列相似性S_seq得到最终的节点相似性S(u,v)＝(1-α)×S_str(u,v)+α×S_seq(u,v)，其中α表示平衡结构相似性和度列相似性的平衡参数，值在[0,1]之间可调节。

然后根据节点间的相似性分数S(u,v)将所有节点对进行降序排列，当节点对(u,v)满足以下两个条件时则将这两个节点进行匹配并加入锚点集合A；

条件1：

条件2：S(u,v)≥τ，其中τ是阈值；

步骤1.5：从选择出的锚点集合A出发将邻居节点进行扩展匹配，得到初始匹配M；

首先将锚点集合A中的匹配节点对都加入初始匹配M中，然后将锚点集合A中的所有匹配节点对(u,v)的邻居节点的笛卡尔积节点对(N(u)×N(v))加入优先队列Q中，并按照扩展相似性S_e(u,v)＝(1-α)×S_l(u,v)+α×S_seq(u,v)大小降序排列，然后逐个出列，如果扩展相似性最大的节点对(u,v)之前都没有与其他节点匹配过，那么将(u,v)匹配对加入M中，并将(u,v)的邻居节点笛卡尔积加入优先队列Q中，直至队列Q为空，得到初始匹配M。

步骤2：优化匹配M得到最优匹配M^*；

匹配优化(matching refinement)阶段，步骤1中根据启发式算法得到的初始匹配M不能保证是最优的，采用一种新提出的算法进行匹配优化，具体是：首先从两个物种其中之一的蛋白质交互网络G₁中随机选择一个顶点覆盖集C，然后保留顶点覆盖集C和初始匹配M的节点交集F₁，以及保留F₁在M中对应的G₂的节点集合F₂，对G₁中不属于F₁集合的节点和G₂中不属于F₂集合的节点进行匹配优化，得到优化后匹配M^*，如果M^*的匹配效果比M好，则更新M，如此多次迭代优化直至M不再更新，就得到接近最优的最终匹配结果。

其中对G₁中不属于F₁集合的节点和G₂中不属于F₂集合的节点进行匹配优化，得到优化后匹配M^*，具体实现包括以下子步骤：

步骤2.1：构建一个带权二分匹配图G_b，一边是包含V₁-F₁的节点结合，另一边是包含V₂-F₂的节点集合，对于所有u∈V₁-F₁和v∈V₂-F₂，在G_b中增加一条边，边(u,v)的权重为w(u,v)＝|M[N(u)∩F₁]∩(N(v)∩F₂)|；

步骤2.2：利用匈牙利树算法计算出G_b的最大带权二分匹配M_b，使得二分图中所有边的权重之和最大化，那么优化后的匹配M^*＝(M∩(F₁×F₂))∪M_b。

判断匹配是否更好的标准是：

其中

如果u₁和u₂之间存在边，则

为1，否则

为0；如果v₁和v₂之间存在边，则

为1，否则

为0。

步骤3：利用匹配M^*发现生物通道；

利用步骤2生成的两个物种的蛋白质交互网络匹配结果，结合生物数据库挖掘两个物种间共存的生物通路。KEGG PATHWAY数据库中存储了现有实验证实的生物通路，其中一个生物通路(例如hsa03010)的表示是由物种代号(hsa表示人类)和一个数字组成(03010)，相同数字的生物通路表示具有相似的生物功能。具体是：首先在KEGG PATHWAY数据库中找到两个物种的蛋白质交互网络(PPI网络)中分别涉及到的所有生物通路，因为相同数字的生物通路表示具有相似的生物功能，所以可以得到两个物种中具有相似生物功能的生物通路集合，然后结合步骤2中生成的匹配结果，可以得到两个物种的具有相似生物功能的生物通路之间的结构映射关系，这个结构映射关系对于生命科学家研究不同物种之间生物通路之间的关联性具有指导意义。

本发明具有以下优点：图匹配算法将生物序列相似性和网络结构相似性很好地融合，能够发现两个物种的蛋白质交互网络中共存的较大子结构，从而更有效地发现存在于不同物种中具有相似功能的生物通路。

附图说明：

图1是本发明实施例的流程图；

图2是本发明实施例在样本图上匹配之后发现的共存的生物通路子结构，其中(1)hsa03010生物通路，(2)sce03010生物通路；

图3是本发明实施例在样本图上发现的共存生物通路被APID dataserver证实的子结构，其中(1)hsa03010生物通路，(2)sce03010生物通路。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明主要基于一种新提出的图匹配算法，将不同物种的蛋白质交互网络(PPI网络)进行最优匹配，根据匹配的结果在KEGG PATHWAY数据库中进行查询，找到物种间共存的生物通路子结构。通过本发明，我们提供了一种新的跨物种生物通路发现方法，比传统生物化学方法更加高效，且算法效果比现有匹配算法更好。

本发明提供的方法能够用计算机软件技术实现流程。参见图1，实施例以人类(human，生物代号HSapiens)的蛋白质交互网络(PPI网络)和酵母菌(yeast，生物代号SCerevisiae)的蛋白质交互网络(PPI网络)为例，样本图的属性信息参见表1，

表1样本图数据(人类和酵母菌的PPI网络)

PPI network	#nodes	#edges	#average degree
				HSapiens(human)	13276	110528	16.651
SCerevisiae(yeast)	5831	77149	26.462

本发明的具体实现包括以下步骤：

步骤1：首先通过计算融合全局相似性、局部相似性、度相似性和序列相似性，得到人类和酵母菌的PPI网络中所有节点间的节点相似性，然后根据节点相似性将所有节点对降序排列，将满足要求的节点对加入锚点集合。然后从锚点集合中的节点对出发，根据局部相似性和序列相似性，扩展匹配它们的邻居节点，直至较小PPI网络中所有节点都被匹配，得到初始匹配。

实施例中构建初始匹配的具体实施过程如下：

首先，计算全局相似性S_g、局部相似性S_l、度相似性S_d和序列相似性S_seq，并融合为节点相似性S。

计算全局相似性S_g的方法为：先计算得到人类PPI网络G₁和酵母菌PPI网络G₂的拉普拉斯矩阵L₁和L₂，然后计算本征值，得到对角矩阵Λ₁和Λ₂，从而根据

和

得到U₁和U₂，则

计算局部相似性S_l的方法为：先计算得到G₁中节点u和G₂中节点v的k步邻居子图N_k(u)和N_k(v)，然后根据发明内容中提到的

计算式计算得到。

计算度相似性S_d的方法为：先计算得到G₁中节点u和G₂中节点v的节点度d(u)和d(v)，然后由

计算得到。

计算序列相似性S_seq的方法为：将BLAST算法计算得到的序列相似性归一化即可。

融合各个相似性为节点相似性S的方法为：

S(u,v)＝(1-α)×S_str(u,v)+α×S_seq(u,v)，

其中S_str(u,v)＝(1-θ)×S_t(u,v)+θ×S_d(u,v)，且S_t＝S_g×S_l，α和β的值均为[0,1]之间可调节，这里均选取为0.5。

然后，将所有节点对根据融合后的节点相似性S大小降序排列，将满足约束条件的节点对加入锚点集合A。这里约束条件为：

(1)

(2)S(u,v)≥τ，其中τ是满足不小于0.5的阈值，这里选取为0.5。

最终，从锚点集合A中的每个节点对出发，扩展匹配它们的邻居节点，直至G₁和G₂中较小图中的所有节点均被匹配，即得到初始匹配M。扩展匹配的具体过程如下：将锚点集合A中的所有匹配节点对(u,v)的邻居节点的笛卡尔积节点对(N(u)×N(v))加入优先队列Q中，并按照扩展相似性S_e(u,v)＝(1-α)×S_l(u,v)+α×S_seq(u,v)大小降序排列，然后逐个出列，如果扩展相似性最大的节点对(u,v)之前都没有与其他节点匹配过，那么将(u,v)匹配对加入M中，并将(u,v)的邻居节点笛卡尔积加入优先队列Q中，直至队列Q为空，即得到初始匹配M。

步骤2：上述启发式算法得到的初始匹配M不一定是最优的，继续将初始匹配M采用覆盖集和匈牙利树二分匹配算法进行优化得到最优匹配M^*。

实施例中优化初始匹配M得到最有匹配M^*的具体实施过程如下：

首先，随机选取G₁的一个顶点覆盖集C，假设初始匹配M中所有G₁节点集合为P₁，令F₁＝C∩P₁，F₁中所有节点由初始匹配M映射的在G₂中的节点结合为F₂。保留F₁和F₂的匹配节点对直接加入到匹配M^*中。

然后，构建一个带权二分图G_b，一边是所有V₁-F₁的节点集合，另一边是所有V₂-F₂的节点集合，二分图所有节点之间加一条边，边的权重为w(u,v)＝|M[N(u)∩F₁]∩(N(v)∩F₂)|。然后利用匈牙利树算法(Hungarian Algorithm)求解最大二分匹配问题，将计算的解，即所有剩余节点的匹配对加入到匹配M^*中。

最终，如果M^*是比M更好的匹配，那么更新M＝M^*，然后重复优化匹配M的所有步骤，直至M^*达到稳定。这里判断匹配是否更好的标准是：

其中

如果u₁和u₂之间存在边则

为1，否则为0，

同理。

步骤3：在KEGG PATHWAY数据库中查询人类PPI网络G₁和酵母菌PPI网络G₂各自涉及到的所有生物通路(pathway)，然后根据最优匹配M^*找到其中具有最大公共子结构的生物通路。

图2给出了KEGG PATHWAY中查询到的人类代号为hsa03010的生物通路，共包含132个蛋白质节点和1924条边，还有查询到的酵母菌代号为sce03010的生物通路，共包含175个蛋白质节点和2311条边，其中最优匹配M^*发现了hsa03010和sce03010的最大公共子结构包含63个蛋白质分子和1406条边，该公共子结构在图2中分别均以深色标出。hsa03010和sce03010具有相同的数字标号03010，表明人类和酵母菌的这两个生物通路具有相似的生物功能。为了进一步证明人类生物通路hsa03010和酵母菌生物通路sce03010之间的关联性，我们采用APID dataserver查询这两个生物通路中已被生物实验证实确实存在的蛋白质及蛋白质之间的联系。APID查询结果表明hsa03010生物通路中总共有126个蛋白质节点和1748条边是被实验证实确实存在的，sce03010生物通路中总共有165个蛋白质节点和192条边是被实验证实确实存在的，它们之间存在一个由26个蛋白质节点和32条边构成的较为完整的公共子结构，该结构展示在图3中。该结果说明本发明提出的图匹配方法可以有效地发现人类和酵母菌中具有相似功能的生物通路公共子结构，该结果对于生物学研究物种之间的联系有指导意义。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (3)

1.一种基于图匹配的跨物种生物通路发现方法，其特征在于，包括以下步骤：

步骤1：构建初始匹配M；

首先将两个物种的蛋白质交互网络G₁和G₂中具有最高全局相似性、局部相似性、序列相似性和度相似性以及度数大于一定阈值的节点进行匹配作为锚点集合A，然后从这些锚点集合A扩展，根据局部相似性和序列相似性匹配锚点的邻居节点，直至两个物种中节点总数较少的蛋白质交互网络中的所有节点均已匹配，得到初始匹配M；

步骤1的具体实现包括以下子步骤；

步骤1.1：采用谱方法计算全局相似性S_g；

对于图G，它的邻接矩阵为G_A，对角度矩阵为D，拉普拉斯矩阵L＝D-G_A；对于两个物种的蛋白质交互网络G₁和G₂，分别计算它们的拉普拉斯矩阵的特征值，假设图G₁的拉普拉斯矩阵L₁的特征值为α₁≥α₂≥…≥α_n，图G₂的拉普拉斯矩阵L₂的特征值为β₁≥β₂≥…≥β_n，令Λ₁＝diag(α_i)，Λ₂＝diag(β_i)，L₁和L₂是对称半正定矩阵，则

其中U₁和U₂是正交矩阵；如果G₁和G₂是同构的，则存在一个列矩阵P使得PL₁P^T=L₂，解得

则全局相似性

步骤1.2：计算局部相似性S_l；

假设图G中的顶点v的k步邻居子图为N_k(v)，并且

表示包括节点v的完整k步子图，V₁表示图G₁的顶点集合，V₂表示图G₂的顶点集合，节点u∈V₁和节点v∈V₂两者之间的局部相似性通过比较u和v的k步邻居子图来衡量，具体如下：假设d(u)、d(v)分别为节点u和节点v在G₁、G₂中的度，假设N_k(u)的所有节点的度按大小降序排列后分别为d_1,1,d_1,2,…，Nk(v)的所有节点的度按大小降序排列后分别为d_2,1,d_2,2,…；令节点u和节点v的k步子图的较小节点总数n_min＝min{|N_k(u)|,|N_k(v)|}，则G₁的节点u和G₂的节点v之间的局部相似性为

其中

和

分别表示包含节点u的k步子图的顶点数和边数，

和

分别表示包含节点v的k步子图的顶点数和边数，其中k步子图的最小度之和

步骤1.3：计算度相似性S_d和序列相似性S_seq；

假设G₁的节点u的度为d(u)，G₂的节点v的度为d(v)，那么u和v之间的度相似性为

通过BLAST计算出序列分数seq(u,v)，然后进行归一化得到序列相似性

步骤1.4：融合各种相似性，选择出锚点集合A；

首先由全局相似性S_g、局部相似性S_l得到拓扑相似性矩阵S_t＝S_g×S_l，再结合度相似性S_d得到网络结构相似性S_str(u,v)＝(1-θ)×S_t(u,v)+θ×S_d(u,v)，其中θ表示拓扑相似性和度相似性的平衡参数，值在[0,1]之间可调节；再结合结构相似性S_str和序列相似性S_seq得到最终的节点相似性S(u,v)＝(1-α)×S_str(u,v)+α×S_seq(u,v)，其中α表示平衡结构相似性和序列相似性的平衡参数，值在[0,1]之间可调节；

然后根据节点间的相似性分数S(u,v)将所有节点对进行降序排列，当节点对(u,v)满足以下两个条件时则将这两个节点进行匹配并加入锚点集合A；

条件1：

条件2：S(u,v)≥τ，其中τ是阈值；

步骤1.5：从选择出的锚点集合A出发将邻居节点进行扩展匹配，得到初始匹配M；

首先将锚点集合A中的匹配节点对都加入初始匹配M中，然后将锚点集合A中的所有匹配节点对(u,v)的邻居节点的笛卡尔积节点对(N(u)×N(v))加入优先队列Q中，并按照扩展相似性S_e(u,v)＝(1-β)×S_l(u,v)+β×S_seq(u,v)大小降序排列，其中β表示平衡S_l(u,v)和S_seq(u,v)的平衡参数；然后逐个出列，如果扩展相似性最大的节点对(u,v)之前都没有与其他节点匹配过，那么将(u,v)匹配对加入M中，并将(u,v)的邻居节点笛卡尔积加入优先队列Q中，直至优先队列Q为空，得到初始匹配M；

步骤2：优化匹配M得到最优匹配M^*；

首先从蛋白质交互网络G₁中随机选择一个顶点覆盖集C，然后保留顶点覆盖集C和初始匹配M的节点交集F₁，以及保留F₁在M中对应的G₂的节点集合F₂，对G₁中不属于F₁集合的节点和G₂中不属于F₂集合的节点进行匹配优化，得到优化后匹配M^*，如果M^*的匹配效果比M好，则更新M，如此多次迭代优化直至M不再更新，就得到接近最优的最终匹配结果；

步骤3：利用匹配M^*发现生物通道；

首先在KEGG PATHWAY数据库中找到两个物种的蛋白质交互网络中所有蛋白质节点涉及的生物通路集合；

然后利用匹配M^*的节点映射关系找到具有最大公共子结构的生物通道。

2.根据权利要求1所述的基于图匹配的跨物种生物通路发现方法，其特征在于：步骤2中所述对G₁中不属于F₁集合的节点和G₂中不属于F₂集合的节点进行匹配优化，得到优化后匹配M^*，具体实现包括以下子步骤：

步骤2.1：构建一个带权二分匹配图G_b，一边是包含V₁-F₁的节点结合，另一边是包含V₂-F₂的节点集合，对于所有u∈V₁-F₁和v∈V₂-F₂，在G_b中增加一条边，边(u,v)的权重为w(u,v)＝|M[N(u)∩F₁]∩(N(v)∩F₂)|；

步骤2.2：利用匈牙利树算法计算出G_b的最大带权二分匹配M_b，使得二分图中所有边的权重之和最大化，那么优化后的匹配M^*＝(M∩(F₁×F₂))∪M_b。

3.根据权利要求1所述的基于图匹配的跨物种生物通路发现方法，其特征在于：步骤2中判断匹配是否更好的标准是：

其中

如果u₁和u₂之间存在边，则

为1，否则

为0；如果v₁和v₂之间存在边，则

为1，否则

为0。

Images