CN104978498A - 生物分子网络拓扑结构比对的自适应方法 - Google Patents

Abstract

本发明涉及一种生物分子网络拓扑结构比对的自适应方法。该方法用于寻找两个生物分子网络生物学意义上的最优映射，其步骤如下：A、构建第一生物分子网络和第二生物分子网络的初始相似矩阵 。B、基于相似矩阵，获得第一生物分子网络中的节点和第二生物分子网络之间的比对映射矩阵。C、根据当前迭代步的匹配结果，自适应地更新相似矩阵，然后计算得到下一迭代步生物分子的相似矩阵。D、计算每一迭代步的映射矩阵的得分，然后判断是否结束计算。其特点在于针对节点的网络拓扑特点自适应地计算网络节点之间的相似性，有效降低仅用匈牙利算法而需要的大量计算时间，同时提高仅用贪心算法的准确率，从而找出比同类算法更好的映射。

Description

生物分子网络拓扑结构比对的自适应方法

技术领域

本发明涉及一种生物分子网络拓扑结构比对的自适应方法。

背景技术

生物分子网络是生物体内各种分子通过相互作用来完成各种复杂的生物功能的一个复杂网络。复杂网络的比对问题是一个图的全局比对问题，涉及到大量的计算，其隐含的子图同构问题已被证实是一个NP完全问题(Non-deterministic Polynomial的问题)，也即是多项式复杂程度的非确定性问题。目前，普遍采用图来表示生物分子网络，并以图论的方法来研究生物分子网络结构，图中的节点表示生物分子，图中的边表示生物分子之间的调控、相互作用各种关系。

目前，生物分子网络仅用图论的方法来研究生物分子网络的比对还面临着绪多的问题，主要包括：

(1)生物分子网络比对主要是基于相似矩阵进行最优匹配，而该最优匹配方法不适合生物网络的比对问题，该方法在匹配中计算复杂，耗时间长；

(2)在生物分子网络比对中，节点的拓扑结构是不同的，匹配节点对其邻居节点匹配的影响强弱不同，现有的技术没有进行自适应的计算。

发明内容

本发明的目的在于克服现有技术存在的不足，提供一种生物分子网络拓扑结构比对的自适应方法，该方法对两个不同物种的复杂生物分子网络拓扑结构进行比对，能找到不同生物分子网络之间的最优匹配，降低自适应的计算时间；该方法在网络变化后，仍可以找到原始生物分子网络的比对，得到更优的比对结果。

为达到上述目的，本发明的构思是：首先使用生物分子网络间的生物序列相似系数和拓扑结构相似系数，迭代计算两个生物分子网络中每一对生物分子的节点与节点之间的相似系数，使用匈牙利算法和贪心算法进行匹配，然后根据第一生物分子网络G_A中的节点和第二生物分子网络G_B中的节点之间的映射矩阵M^k，使用节点对的邻居节点之间的相似性，自适应地构建更新相似矩阵S^k；最后，根据迭代的连续相邻的生物分子的相似矩阵S^k+1和生物分子的S^k相似矩阵S^k匹配的比较，判断是否结束运算，结束时，获取结果最优的生物分子网络的相似矩阵S^k+1的相似参数。

根据上述发明构思，本发明采用下述技术方案：

生物分子网络拓扑结构比对的自适应方法，其特征在于该方法的具体操作步骤如下：

a.构建两个分子网络的初始相似矩阵S⁰，该两个分子网络分别记为第一生物分子网络G_A和第二生物分子网络G_B，所述的初始相似矩阵S⁰中的S⁰(a_i,b_j)表示节点a_i和节点b_i之间的相似系数，其中a_i∈G_A，b_j∈G_B a_i表示第一生物分子网络G_A中的节点，b_j表示第二生物分子网络G_B中的节点；

b.根据当前迭代步的相似矩阵S^k，使用匈牙利算法、贪心算法对第一生物分子网络G_A中的节点和第二生物分子网络G_B中的节点进行匹配，获得G_A中的节点和G_B中的节点之间的映射矩阵M^k，其中k表示迭代步数，初始时k＝0：其中每一个元素M^k(a_i,b_j)为1或0，1表示节点a_i和节点b_i匹配，0则表示不匹配；

c.根据步骤b所得的映射矩阵M^k，使用节点的邻居节点之间的相似性，自适应地更新相似矩阵S^k，然后结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，计算得到下一生物分子的相似矩阵S^k+1；

d.计算每一迭代步的映射矩阵M^k的得分SS^k，然后计算生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值，判断是否结束迭代计算，若生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值小于阈值λ，λ体现允许的计算误差，则结束计算，取第m(0＜＝m＜＝k)步的映射矩阵M^m为最终映射结果，其中m需满足该步的映射矩阵M^m的得分SS^m最大；否则，若生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值不小于阈值λ，则不结束计算，返回步骤B到步骤F继续进行计算，直到前后两次生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值小于设定的阈值λ，则结束计算；

根据权利要求1所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述的步骤b的具体方法为：

b-1、遍历G_A的每一个节点，如果这个节点在网络G_B中的相似节点数不小于一个阈值该阈值为一经验值，就把这个节点在相似矩阵S^k中对应的行归入到匈牙利矩阵H^k中，否则归 入贪心矩阵G^k中；

b-2、使用匈牙利算法处理匈牙利矩阵H^k，并在G^k中标记那些已经匹配的列；如果在有效的时间内获得结果，则进行下一步，否则结束处理，把匈牙利矩阵H^k合并到贪心矩阵G^k中；

b-3、使用贪心算法处理贪心矩阵G^k中由那些没有被标记的列构成的矩阵，和b-2中匈牙利矩阵H^k的结果合并，获得最终映射M^k。

上述步骤c的迭代计算G_A和G_B的相似矩阵S^k其具体计算方式采用下述公式：

$S^k(a_i,b_j)=S^k(a_i,b_j)+\frac{\underset{a_u\↔a_i,b_v\↔b_j}{\Σ}{S}^k(a_u,b_v)}{\mathrm{deg}\left(a_i\right)},M^k(a_u,b_v)=1$

其中，表示a_u和a_i之间存在边，表示b_v和b_j之间存在边，a_u和b_v是在当前迭代步完成后相匹配的节点，deg(a_i)为节点a_i的度；结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，迭代计算得到生物分子的相似矩阵S^k+1。

上述的迭代计算得到两个分子网络的相似矩阵S^k+1的具体步骤如下：

C-1.计算两个分子网络节点a_i和节点b_i在网络拓扑结构上的平均相似性，其相似性的各个方面由N₁(a_i,b_j)和N₂(a_i,b_j)表示；

C-2.在两个分子网络节点节点a_i和节点b_i的初始相似系数S⁰(a_i,b_j)基础上集成它们在网络拓扑结构上的平均相似性，计算S^k+1，其中上标k和k+1代表迭代次数，具体公式如下：

$S^{k+1}(a_i,b_j)=\frac{N_1(a_i,b_j)+N_2(a_i,b_j)}{2}\×(1-\∂)+S^0(a_i,b_j)\×\∂$

归一化，具体公式为：

$S^{k+1}=\frac{S^{k+1}}{max\left(S^{k+1}\right)}$

N₁(a_i,b_j)表示a_i和b_j的邻居节点之间的平均相似性，N₂(a_i,b_j)表示非邻居节点之间的平均相似性，max(S^k+1)表示矩阵S^k+1的最大值，表示一个0到1之间的值，由用户输入确 定，用以控制节点生物相似性和拓扑相似性的权重。

上述的计算生物分子a_i和b_j在网络拓扑结构上的平均相似性的具体方法为：

其中，N₁(a_i,b_j)表示节点a_i和b_j的邻居节点，即这些节点和a_i或b_j之间存在边，之间的平均相似性，N₂(a_i,b_j)表示节点a_i和b_j的非邻居节点，即这些节点和a_i或b_j之间不存在边，之间的平均相似性，由deg(a_i)表示节点a_i的度，表示a_u和a_i之间存在边，表示a_u和a_i之间不存在边，a_u∈G_A表示a_u是G_A中的一个节点，b_v∈G_B表示b_v是G_B中的一个节点，S(a_u,b_v)表示a_u和b_v的相似系数，n₁为G_A的节点数，n₂为G_B的节点数，上标k代表迭代次数。

上述的映射矩阵的得分SS^k侧重于生物意义的打分PES^k和表示拓扑意义的边正确率EC^k计算得到，其具体计算方式如下：

SS^k＝PES^k+100×EC^k

其中PES^k和EC^k的计算方式如下：设网络A为G_A＝(V₁,E₁)，网络B为G_B＝(V₂,E₂)，其中，V₁，V₂分别代表网络G_A，G_B的节点集合，且|V₁|＝n₁，|V₂|＝n₂，即G_A中有n₁个节点，G_B有n₂个节点；a_i,a_u∈V₁表示节点a_i，a_u存在于G_A中，分别表示在第k次迭代 时G_B中与a_i，a_u对应的匹配节点；E₁，E₂分别代表G _A，G_B的边集合，(a_i,a_u)＝e_iu表示边e_iu的两个端点是节点a_i，a_u，e_iu∈E₁表示边e_iu是G_A的一条边；表示G_A的节点a_i和G_B中与其对应的节点的初始相似系数。G_A与G_B的PES^k为：

${\mathrm{PES}}^k=\frac{1}{2}\underset{i,u=1}{\overset{n_1}{\Σ}}{\mathrm{ES}}^k\left(e_{iu}\right)+\underset{i=1}{\overset{n_1}{\Σ}}{\mathrm{PS}}^k\left(a_i\right)$

其中

EC^k为：

上述的结束计算的判断，有两个依据条件，只要满足一个，即可结束计算：

D1、S^k+1与S^k或S^k与S^k-1对应元素差值绝对值的最大值小于一个限定的阈值(0.01)：

max(|S^k+1-S^k|)＜＝0.01

或

max(|S^k-S^k-1|)＜＝0.01

D2、映射矩阵的得分连续三步不变化：

SS^k+1(G_A,G_B)＝SS^k(G_A,G_B)＝SS^k-1(G_A,G_B)。

上述步骤B中所述的匹配，结合了匈牙利算法和贪心算法各自的优点。匈牙利算法在处理此类问题中是公认的优秀算法，可以得到较优的结果，但是由于生物分子网络具有生物背景，两个网络之间的序列相似节点非常少，即G_A中的一个节点，在G_B中可能只有少数的几个节点和其具有序列相似性，在相似矩阵中表现为很多行都只有一个或几个有效值，其它位 置都是0，对于这种稀疏矩阵，匈牙利算法很难在有效的时间内获得结果；贪心算法没有时间问题，但是结果又不理想。为了能获得相似性更高的匹配，并有效降低计算时间，本发明在匹配时对相似矩阵进行了处理。

所谓邻居优先策略是指在生物分子网络比对中，已匹配节点的邻居节点更相似，即如果一对节点相匹配，那么它们的邻居节点间的相似性应该增加，以增加这些邻居节点匹配的可能性。对于一对节点a_i∈G_A和b_j∈G_B，它们所有邻居节点的匹配情况都应该对这对节点的匹配产生影响，但是每一个节点在网络中的拓扑位置是不同的，应该根据节点自身的拓扑位置及其邻居节点的匹配情况而进行自适应的计算。

根据当前迭代步的匹配结果，使用自适应邻居优先策略，更新相似矩阵。。为了避免由于节点度过大使得这样影响太过明显以削弱其本身相似性的作用，应该进行控制，

本发明的生物分子网络比对自适应混合方法，与现有技术相比，具有以下突出的实质性特点和显著优点：

1、该方法结合匈牙利算法和贪心算法中各自的优点，得到第一生物分子网络G_A中的节点和第二生物分子网络之间的比对映射及相应的比对映射矩阵M^k，能在有效时间内找到较好的匹配结果，降低计算时间。

2、该方法能根据生物分子在网络中的拓扑位置的不同，自适应地结合生物分子属性和拓扑特点，计算生物分子的相似性，并更能符合生物分子网络的生物学意义。

3、该方法可以找到专家知识词典中给出的最优匹配，且匹配结果和有专家知识词典的匹配结果基本一致，其匹配结果的匹配蛋白质对中拥有更多的GO共同项，能够进行更好的生物功能预测。

4、该方法将生物分子网络的生物学相似性和拓扑结构相似性进行了有效的结合，使比对结果在生物学意义和拓扑结构上都具有优势。

附图说明

图1是本发明的生物分子网络比对自适应混合方法的流程图。

图2是图1步骤B所述的结合匈牙利算法和贪心算法进行匹配的具体流程图。

图3是图1中步骤C所述的基于当前步的匹配结果，使用自适应邻居优先方法进行相似矩阵更新的具体流程图。

图4是是本发明与同类型算法结果对比。

图5是为了不失一般性，对图5算例的G_A边数的5％、10％、15％、20％、25％、30％、35％、 40％各进行50次拓扑变换后，本发明与同类方法的结算结果对比。

图6是本发明与只使用匈牙利算法和只使用贪心算法的结果对比。

图7是本发明与INM算法在专家知识词典和共同GO项方面的对比。

具体实施方式

以下结合附图对本发明的优选实施例进一步详细说明。

本实施例中，本发明的生物分子网络比对自适应混合方法的实验在上海大学系统生物技术研究所的集群计算机上完成，该集群由14台IBM HS21刀片服务器和2台x3650服务器组成计算和管理节点，网络连接采用千兆以太网和infiniband 2.5G网。每个节点配置两个双核CPU和4GB内存，每个CPU为intel xeon 5150 2.66GMhz主频，两台图形工作站作为前端机，可以进行科学数据可视化。

本实施实例在酵母和人类网络上进行。本发明的生物分子网络比对自适应混合方法，如图1－图3所示，包括以下步骤：

A、构建第一生物分子网络和第二生物分子网络的初始相似矩阵S⁰：假设第一生物分子网络记为G_A，第二生物分子网络记为G_B，构建第一生物分子网络G_A和第二生物分子网络G_B中生物分子的初始相似矩阵，记为S⁰，其中，S⁰(a_i,b_j)表示节点a_i∈G_A和节点b_j∈G_B之间的相似系数，a_i表示第一生物分子网络G_A中的节点，b_j表示第二生物分子网络G_B中的节点；

B、构建第一生物分子网络G_A和第二生物分子网络G_B节点间的映射矩阵M^k，其中k表示迭代步数，初始时k＝0：根据当前迭代步的相似矩阵S^k，使用匈牙利算法、贪心算法对第一生物分子网络G_A中的节点和第二生物分子网络G_B中的节点进行匹配，获得G_A中的节点和G_B中的节点之间的映射矩阵M^k，其中每一个元素M^k(a_i,b_j)为1表示节点a_i∈G_A和节点b_j∈G_B匹配，为0则不匹配：

B1、遍历G_A的每一个节点，如果这个节点在网络G_B中的相似节点数不小于一个阈值(我们在进行蛋白质相互作用网络比对的时候，经验值为5)，就把这个节点在相似矩阵S^k中对应的行归入到匈牙利矩阵H^k，否则归入贪心矩阵G^k。

B2、使用匈牙利算法处理匈牙利矩阵H^k，并在G^k中标记那些已经匹配的列。如果在有 效的时间内获得结果，则进行下一步，否则结束处理，把匈牙利矩阵H^k合并到贪心矩阵G^k中。

B3、使用贪心算法处理贪心矩阵G^k中由那些没有被标记的列构成的矩阵，和之前的匈牙利矩阵结果合并，获得最终映射M^k。

C、迭代计算第一生物分子网络G_A和第二生物分子网络G_B的相似矩阵S^k：首先根据当前迭代步的匹配结果M^k，使用节点的邻居节点之间的相似性，自适应地更新相似矩阵S^k，然后结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，计算得到下一迭代步生物分子的相似矩阵S^k+1：

C1、根据当前迭代步的匹配结果，使用自适应邻居优先策略，更新相似矩阵。所谓邻居优先策略是指在生物分子网络比对中，已匹配节点的邻居节点更相似，即如果一对节点相匹配，那么它们的邻居节点间的相似性应该增加，以增加这些邻居节点匹配的可能性。对于一对节点a_i∈G_A和b_j∈G_B，它们所有邻居节点的匹配情况都应该对这对节点的匹配产生影响，但是每一个节点在网络中的拓扑位置是不同的，应该根据节点自身的拓扑位置及其邻居节点的匹配情况而进行自适应的计算。为了避免由于节点度过大使得这样影响太过明显以削弱其本身相似性的作用，应该进行控制，所以本发明使用邻居节点的相似性及节点a_i本身的度进行自适应的计算，其具体计算方法如下：

$S^k(a_i,b_j)=S^k(a_i,b_j)+\frac{\underset{a_u\↔a_i,b_v\↔b_j}{\Σ}{S}^k(a_u,b_v)}{\mathrm{deg}\left(a_i\right)},M^k(a_u,b_v)=1$

其中，表示a_u和a_i之间存在边，表示b_v和b_j之间存在边，a_u和b_v是在当前迭代步完成后相匹配的节点，deg(a_i)为节点a_i的度；

C2、结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，迭代计算得到生物分子的相似矩阵S^k+1，具体步骤如下：

C21、计算生物分子a_i∈G_A和b_j∈G_B在网络拓扑结构上的平均相似性，其相似性的各个方面由N₁(a_i,b_j)和N₂(a_i,b_j)表示；

C22、在生物分子a_i∈G_A和b_j∈G_B的初始相似系数S⁰(a_i,b_j)基础上集成它们在网络拓扑 结构上的平均相似性，计算S^k+1；具体公式如下：

$S^{k+1}(a_i,b_j)=\frac{N_1(a_i,b_j)+N_2(a_i,b_j)}{2}\×(1-\∂)+S^0(a_i,b_j)\×\∂$

其中上标k和k+1代表迭代次数，N₁(a_i,b_j)表示节点a_i和b_j的邻居节点(即这些节点和a_i或b_j之间存在边)之间的平均相似性，N₂(a_i,b_j)表示节点a_i和b_j的非邻居节点(即这些节点和a_i或b_j之间不存在边)之间的平均相似性，max(S^k+1)表示矩阵S^k+1的最大值，表示一个0到1之间的值，由用户输入确定，用以控制节点生物相似性和拓扑相似性的权重，deg(a_i)表示节点a_i的度，表示a_u和a_i之间存在边，表示a_u和a_i之间不存在边，a_u∈G_A表示a_u是G_A中的一个节点，b_v∈G_B表示b_v是G_B中的一个节点，S(a_u,b_v)表示a_u和b_v的相似系数，n₁为G_A的节点数，n₂为G_B的节点数，上标k代表迭代次数，N₁(a_i,b_j)及N₂(a_i,b_j)的数学定义如下：

C23、归一化，具体公式为：

$S^{k+1}=\frac{S^{k+1}}{max\left(S^{k+1}\right)}$

D、判断计算是否结束：首先计算每一迭代步的映射矩阵M^k的得分SS^k，然后计算迭代运算中的两个连续相邻的生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值最大值，判断是否结束生物分子网络的相似矩阵S^k+1的计算。若生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值最大值小于阈值λ，则结束计算，取第m(0＜＝m＜＝k)步的映射矩阵M^m为最终映射结果，其中m需满足该步的映射矩阵M^m的得分SS^m最大；否则，若生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值最大值不小于阈值λ，此处λ取0.01，则不结束计算，返回步骤B继续进行计算，直到前后两次生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值最大值小于设定的阈值λ，则结束计算。或者，SS^k连续三次不变化，则结束计算。以上结束计算的条件，只要满足一个，即可结束计算：

D1、S^k+1与S^k或S^k与S^k-1对应元素差值绝对值的最大值小于一个限定的阈值(0.01)：

max(|S^k+1-S^k|)＜＝0.01

或

max(|S^k-S^k-1|)＜＝0.01

D2、映射结果的得分连续三步不变化：

SS^k+1(G_A,G_B)＝SS^k(G_A,G_B)＝SS^k-1(G_A,G_B) 

参照图4，示出了本发明的生物分子网络比对自适应混合方法与同类代表性方法，即INM、NBM和MI-GRAAL，对于网络比对例的计算结果对比。图4A是两个网络的来源以及网络的边数节点数，图4B是示例的网络G_A，网络G_B由于规模太大，很难进行可视化展示而没有给出。图4C给出了四个算法比对结果的侧重生物意义的打分PES^k和表示拓扑意义的边正确率EC^k，以及表示拓扑意义的最大共同连通子网(LCCS)的节点数和边数。计算结果表明，本发明的生物分子网络比对自适应混合方法，其比对结果不管是在侧重生物意义的评判标准PES得分上，还是表示拓扑结构意义的EC值上，比起INM算法、NBM算法和 MI-GRAAL算法都有很大的优势。在比对结果的最大连通子网中，本发明的生物分子网络比对自适应混合方法得到的节点数和INM相同，都优于NBM，而在边数上大大高于INM和NBM算法，MI-GRAAL算法的结果最大连通子网虽然节点数多了一个，但是边数非常少，也就是说结果的连通性非常差，这样的比对结果也是不可取的。可见本发明的生物分子网络比对自适应混合方法能比同类方法得到生物学意义和拓扑结构意义都更相似的比对，在准确性方面本方法更优。

本发明与NBM采用了不同的邻居优先的方法，下面图5的实验主要针对NBM进行。

参照图5，示出了对图4算例的7～56条边(约占总边数的5％～40％)，各进行50次拓扑变换获得的400个不同的网络，本发明的生物分子网络比对自适应混合方法与同类代表性方法NBM的计算结果对比。取G_A和G_B均为图4A中的两个网络，以改变网络的拓扑结构，代表物种的进化，具体为以5％(7条边)为增量，分别修改G_A的7，14，……，56条边各50次，获得400个不同拓扑结构的网络；同源表为酵母和人类蛋白质间的序列相似性计算获得；生物相似性权重参数α取0.1～1。图7示出了α＝0.4时，本发明和NBM对这400个不同网络相对于G_B的计算结果，当α取其它值时具有类似特性。其中：

图5A为变换G_A的7～56条边各50次后的总分平均分，其中横轴表示变换G_A的边数，纵轴表示每次变换边后得到的PES平均值。实验表明，对于这400个网络，本发明的比对结果的总分平均分普遍高于NBM方法，说明本发明的比对结果在生物学意义上普遍更相似。

图5B为变换G_A的7～56条边各50次后的EC平均值，其中横轴表示变换G_A的边数，纵轴表示每次变换边后得到的EC平均值。实验表明，对于这400个网络，本发明的比对结果的EC平均值普遍高于NBM方法，说明本发明的比对结果在拓扑意义上也普遍更相似。

图5C和5D为变换G_A的7～56条边各50次后的比对结果最大连通子网的节点数和边数的平均值。其中横轴表示变换G_A的边数，纵轴表示每次变换边后得到的最大共同连通子网(LCCS)边数和节点数的平均值。实验表明，对于这400个网络，本发明的比对结果的LCCS节点数和边数普遍高于NBM方法，说明本发明的比对结果在拓扑意义上普遍更相似。

图5E为变换G_A的7～56条边各50次后比对结果中的保守匹配对(和初始G_A的比对结果相同的匹配节点对数目)的平均值，其中横轴表示变换G_A的边数，纵轴表示每次变换边后得到的保守匹配对的平均值。实验表明，对于这400个网络，本发明的比对结果找到的保守匹配对普遍高于NBM方法，说明本发明在物种进化时更能找到保守的匹配，具有更高的稳定性。

总之，图5的实验结果表明，本发明的生物分子网络比对自适应混合方法与同类代表性方法NBM相比，比对结果在生物学意义和拓扑结构意义上都更相似；随着生物分子网络的进化，本发明更能找到保守的匹配对，具有更高的稳定性。

参照图6，示出了对图4算例的网络以及另外两个同类网络使用本发明和仅使用贪心算法或仅使用匈牙利算法结果比较。图6A为用来比对的网络的来源及网络的边数节点数。6B给出了各个算法在不同数据下比对结果的侧重生物意义的打分PES^k和表示拓扑意义的边正确率EC^k，以及表示拓扑意义的最大共同连通子网(LCCS)的节点数和边数。计算结果表明，本发明的生物分子网络比对自适应混合方法，和仅使用贪心算法时，对所有数据都可以在有效的时间内获得结果，但是本发明的比对结果明显优于贪心算法的结果，匈牙利算法更多时候不能在有效的时间内获得结果，可以算出结果的时候，本发明也可以获得和匈牙利算法一样好的结果。说明本发明可以在有效的时间内获得较优的匹配结果。

参照图7，它是果蝇和人类的蛋白质相互作用网络比对中，本发明和INM算法比对结果得到的和专家知识词典中一致的最优匹配数以及比对结果的共同GO项。

图7A是本次试验的G_A。

图7B是专家知识词典给出的三对最优匹配。

图7C是在不是用专家知识词典进行强制匹配的时候，本发明的生物分子网络比对自适应混合方法及INM算法找到的和专家知识词典中相同的最优匹配数。

图7D是比对结果中拥有至少1个，2个，3个，4个共同GO项的匹配对数。

图7的实验说明，本发明的生物分子网络比对自适应混合方法在不使用专家知识词典进行强制匹配的时候，就可以找到专家知识词典中给出的最优匹配，且比对结果中更多的匹配对拥有共同的GO项。

综上所述，图4～图7表明，本发明的生物分子网络比对自适应混合方法，与同类代表性方法相比，可以在找到较优的匹配，并有效降低计算时间，能找到在生物学意义和拓扑结构意义上都更相似的比对结果，更好地应对物种的进化，在物种进化网络拓扑变化后，可以找到更多的保守匹配节点对，具有更高的稳定性。

本文结合说明书附图和具体实施例进行阐述只是用于帮助理解本发明的方法和核心思想。本发明所述的方法并不限于具体实施方式中所述的实施例，本领域技术人员依据本发明的方法和思想得出的其它实施方式，同样属于本发明的技术创新范围。本说明书内容不应理解为对本发明的限制。

Claims (7)

1.生物分子网络拓扑结构比对的自适应方法，其特征在于该方法的具体操作步骤如下：

a.构建两个分子网络的初始相似矩阵S⁰，该两个分子网络分别记为第一生物分子网络G_A和第二生物分子网络G_B，所述的初始相似矩阵S⁰中的S⁰(a_i,b_j)表示节点a_i和节点b_i之间的相似系数，其中a_i∈G_A，b_j∈G_B a_i表示第一生物分子网络G_A中的节点，b_j表示第二生物分子网络G_B中的节点；

b.根据当前迭代步的相似矩阵S^k，使用匈牙利算法、贪心算法对第一生物分子网络G_A中的节点和第二生物分子网络G_B中的节点进行匹配，获得G_A中的节点和G_B中的节点之间的映射矩阵M^k，其中k表示迭代步数，初始时k＝0：其中每一个元素M^k(a_i,b_j)为1或0，1表示节点a_i和节点b_i匹配，0则表示不匹配；

c.根据步骤b所得的映射矩阵M^k，使用节点的邻居节点之间的相似性，自适应地更新相似矩阵S^k，然后结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，计算得到下一生物分子的相似矩阵S^k+1；

d.计算每一迭代步的映射矩阵M^k的得分SS^k，然后计算生物分子的相似矩阵S^k+1和S^k之间对应元素差值的绝对值，判断是否结束迭代计算，若生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值小于阈值λ，λ为允许的计算误差，则结束计算，取第m，0≤m≤k步的映射矩阵M^m为最终映射结果，其中m需满足该步的映射矩阵M^m的得分SS^m最大；否则，若生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值不小于阈值λ，则不结束计算，返回步骤b到步骤d继续进行计算，直到前后两次生物分子的相似矩阵S^k+1和S^k之间对应元素的差值的绝对值的最大值小于阈值λ，则结束计算。

2.根据权利要求1所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述的步骤b的具体方法为：

b-1、遍历G_A的每一个节点，如果这个节点在网络G_B中的相似节点数不小于一个阈值(该阈值为一经验值)，就把这个节点在相似矩阵S^k中对应的行归入到匈牙利矩阵H^k中，否则归 入贪心矩阵G^k中；

b-2、使用匈牙利算法处理匈牙利矩阵H^k，并在G^k中标记那些已经匹配的列；如果在有效的时间内获得结果，则进行下一步，否则结束处理，把匈牙利矩阵H^k合并到贪心矩阵G^k中；

b-3、使用贪心算法处理贪心矩阵G^k中由那些没有被标记的列构成的矩阵，和b-2中匈牙利矩阵H^k的结果合并，获得最终映射M^k。

3.根据权利要求1所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述步骤c的迭代计算G_A和G_B的相似矩阵S^k其具体计算方式采用下述公式：

其中，表示a_u和a_i之间存在边，表示b_v和b_j之间存在边，a_u和b_v是在当前迭代步完成后相匹配的节点，deg(a_i)为节点a_i的度；结合生物分子的初始生物相似性和生物分子在各自网络中的拓扑相似特征，迭代计算得到生物分子的相似矩阵S^k+1。

4.根据权利要求3所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述的迭代计算得到两个分子网络的相似矩阵S^k+1的具体步骤如下：

C-1.计算两个分子网络节点a_i和节点b_i在网络拓扑结构上的平均相似性，其相似性的各个方面由N₁(a_i,b_j)和N₂(a_i,b_j)表示；

C-2.在两个分子网络节点节点a_i和节点b_i的初始相似系数S⁰(a_i,b_j)基础上集成它们在网络拓扑结构上的平均相似性，计算S^k+1，其中上标k和k+1代表迭代次数，具体公式如下：

归一化，具体公式为：

N₁(a_i,b_j)表示a_i和b_j的邻居节点之间的平均相似性，N₂(a_i,b_j)表示非邻居节点之间的 平均相似性，max(S^k+1)表示矩阵S^k+1的最大值，表示一个0到1之间的值，由用户输入确定，用以控制节点生物相似性和拓扑相似性的权重。

5.根据权利要求4所述的生物分子网络拓扑结构比对的自适应方法，其特征在于中所述的计算生物分子a_i和b_j在网络拓扑结构上的平均相似性的具体方法为：

其中，N₁(a_i,b_j)表示节点a_i和b_j的邻居节点，即这些节点和a_i或b_j之间存在边，之间的平均相似性，N₂(a_i,b_j)表示节点a_i和b_j的非邻居节点，即这些节点和a_i或b_j之间不存在边，之间的平均相似性，由deg(a_i)表示节点a_i的度，表示a_u和a_i之间存在边，表示a_u和a_i之间不存在边，a_u∈G_A表示a_u是G_A中的一个节点，b_v∈G_B表示b_v是G_B中的一个节点，S(a_u,b_v)表示a_u和b_v的相似系数，n₁为G_A的节点数，n₂为G_B的节点数，上标k代表迭代次数。

6.根据权利要求1所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述的映射矩阵的得分SS^k侧重于生物意义的打分PES^k和表示拓扑意义的边正确率EC^k计算得到，其具体计算方式如下：

SS^k＝PES^k+100×EC^k

其中PES^k和EC^k的计算方式如下：设网络A为G_A＝(V₁,E₁)，网络B为G_B＝(V₂,E₂)，其中，V₁，V₂分别代表网络G_A，G_B的节点集合，且|V₁|＝n₁，|V₂|＝n₂，即G_A中有n₁个节点，G_B有n₂个节点；a_i,a_u∈V₁表示节点a_i，a_u存在于G_A中，分别表示在第k次迭代时G_B中与a_i，a_u对应的匹配节点；E₁，E₂分别代表G_A，G_B的边集合，(a_i,a_u)＝e_iu表示边e_iu的两个端点是节点a_i，a_u，e_iu∈E₁表示边e_iu是G_A的一条边；表示G_A的节点a_i和G_B中与其对应的节点的初始相似系数。G_A与G_B的PES^k为：

其中

EC^k为：

7.根据权利要求1所述的生物分子网络拓扑结构比对的自适应方法，其特征在于所述的结束计算的判断，有两个依据条件，只要满足一个，即可结束计算：

D1、S^k+1与S^k或S^k与S^k-1对应元素差值绝对值的最大值小于一个限定的阈值(0.01)：

max(|S^k+1-S^k|)＜＝0.01

或

max(|S^k-S^k-1|)＜＝0.01

D2、映射矩阵的得分连续三步不变化：

SS^k+1(G_A,G_B)＝SS^k(G_A,G_B)＝SS^k-1(G_A,G_B)。