基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法

摘要

本发明公开了基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法。其步骤如下：A、根据生物分子的序列特征，构建网络T（Gt）和B（GB）中生物分子的初始相似矩阵。B、根据生物分子在网络中的拓扑的相似特征，计算生物分子的相似矩阵S。C、构建专家知识字典。D、基于专家知识，以邻居优先策略进行搜索，获得结果子网。E、计算结果子网R（Gr）与目标子网T（Gt）的相似得分。F、计算p值。G、结果子网可视化。与同类方法相比，该方法计算准确度更高，具有更高的稳定性，能更好地处理生物分子网络间的“空隙”（Gap），获得更多的保守边和节点，且能较好地反映生物分子网络之间的变化情况。

说明书

技术领域

本发明涉及计算机领域，提出了基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法。

背景技术

生物分子网络是复杂网络。在复杂网络中搜索与目标子网最相似的子网是一个局部网络比较问题，涉及到大量的计算，已被证实是一个NP完全问题（Non-deterministic Polynomial的问题，也即是多项式复杂程度的非确定性问题）。目前，研究人员普遍采用图来表示复杂网络，并以图论的方法来研究它们。对于生物分子网络而言，图中的节点表示生物分子，边表示生物分子之间的调控、相互作用等各种关系。

由于生物分子网络特有的生物学意义，仅用图论的方法来研究它们是不够的，其搜索还面临着更多的问题，主要包括：（1）每个生物分子都有其生物学意义，要明确一个网络中的某个生物分子和另外一个网络中的哪个生物分子最相似，不仅要考虑生物分子本身的序列，还要考虑它在网络中的拓扑位置；（2）无论是国际公开的数据库中的数据，还是自己通过生物实验获得的数据都存在假阳性和假阴性现象，目前只能通过这些不完全准确和不完整的数据研究生物分子网络；（3）对于要研究的不同的具体问题，网络中各个分子的地位并不是完全平等的，计算过程中要合理利用专家知识，以贴近生物学的实际应用背景。

目前，已有一些研究小组在进行这方面的研究，也开发了少量的工具。这些方法各有其优点，但也各有其局限性，无法满足系统生物学的需要。而这些局限性主要体现在对具有相对复杂的拓扑结构的网络搜索上，对于生物分子网络而言，为了能获得不同物种间的最相似网络，计算的准确度必须得到提高。同时，考虑到生物分子网络的进化和变异，不同物种的网络虽然不同，却有一定的保守性，算法应能在变异后的网络中找到原始的保守信息，能较好地体现网络拓扑的变化情况，且具有较高的稳定性。

发明内容

本发明的目的在于，为了解决上述问题而提供基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法，该方法能在复杂的生物分子网络中搜索到与目标子网最相似的结果子网，避免了以往算法不能合理利用专家知识而带来的计算的盲目性，并降低因原始信息缺失带来的误差，从而具有较高的稳定性。随着生物分子网络的进化和变异，该方法较少受到Gap的影响，从而可以搜索得到更多的保守边和节点，而变异的边数往往与未匹配边数一致，即方法能较好地体现拓扑的变化情况。

为达到上述目的，本发明的构思是：首先结合生物分子的序列特征及其在网络中的拓扑相似特征，计算各个生物分子之间的相似系数，其中拓扑相似特征主要考虑目标生物分子的邻居/非邻居分子之间的平均相似性，以降低原始信息缺失和不准确带来的误差，并提高其稳定性；然后根据具体问题和专家知识字典，将目标子网中的生物分子分类，确定K类分子的最相似分子；最后，根据生物分子之间的关系特点，如“与相似的蛋白质发生相互作用的那些蛋白质之间往往具有更高的相似度”，对N类分子采用邻居优先的策略进行搜索，获得结果子网。

根据上述发明构思，对于网络A（G_A）、网络B（G_B）及网络A中的目标子网T（G_t），本发明采用下述技术方案：

A、 计算G_t和G_B的初始相似矩阵                                                ：根据生物分子的序列特征，构建网络T（G_t）和B（G_B）中生物分子的初始相似矩阵，其中的每个元素表示节点和节点之间的序列相似系数；

B、 计算G_t和G_b的相似矩阵S：根据生物分子在各自网络中的拓扑相似特征，计算生物分子的相似矩阵S，矩阵中的每个元素为节点和节点之间的相似系数；

C、  构建专家知识字典：字典中包含了网络T（G_t）和B（G_B）中由专家确定的最相似的生物分子对；

D、  采用邻居节点优先策略进行网络搜索：利用专家知识，基于相似矩阵S，以邻居优先策略进行搜索，获得结果子网；

E、  计算结果子网（G_r）与目标子网（G_t）的相似得分；

F、  计算p值，分析目标子网的统计学意义；

G、 结果子网（G_r）可视化。

本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法，与现有技术相比较，具有以下突出的实质性特点和显著优点：

1.   该方法建立专家知识字典，避免了以往算法不能合理利用专家知识而带来的计算的盲目性。

2.   该方法结合生物分子本身的序列特征及其在网络中的拓扑特征计算生物分子的相似系数，实现了图论方法和生物学应用背景的有机结合。

3.   该方法在计算生物分子拓扑结构相似的时候，强调生物分子在网络拓扑结构上的平均相似性，而弱化它们的不相似性，有效降低了因为原始数据的不准确和不完整带来的误差。提高了算法的稳定性，且算法较少受到Gap的影响，能较好地体现网络的变化情况。

4.   该方法采用邻居节点优先进行网络搜索，符合生物分子网络的生物学意义，并降低了计算复杂度。

附图说明

图1是本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法的流程图。

图2是图1中步骤B所述的根据生物分子在各自网络中的拓扑相似特征，对无向图计算生物分子的相似矩阵的具体流程图，对有向图的计算与此类似。

图3是图1中步骤D所述的基于专家知识进行搜索的流程图。

图4是图3中步骤D3所述的对N类生物分子根据邻居优先的策略进行搜索配对的流程图。

图5是图1中步骤F所述的计算p值的流程图。

图6是图5中步骤F1所述的生成随机网络的流程图。

图7是本发明与同类方法对经典示例的计算结果对比。

图8是为了不失一般性，对图7算例的1~7条边各进行100次拓扑变换后，本发明与同类方法的计算结果对比图。

图9是为了不失一般性，各以最多100种方式删除图7算例的1~6个节点后，本发明与同类方法的计算结果对比图。

图10是为了不失一般性，对图7算例增加节点后，本发明与同类方法的计算结果对比图。

图11是对果蝇和人类网络搜索比对时，采用或不采用专家知识的结果对比。

具体实施方式

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

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

对于网络A（G_A）、网络B（G_B）及网络A（G_A）中的目标子网T（G_t），本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法，如图1－图6所示，包括以下步骤：

A、 根据生物分子的序列特征，构建网络T（G_t）和B（G_B）中生物分子的初始相似矩阵，其中的每个元素表示节点和节点之间的序列相似系数。具体步骤包括：

A1、取，其余参数取缺省值，用BLAST计算中所有分子在的序列相似分子。

A2、按以下公式计算这些生物分子之间的相似系数：

B、 根据生物分子在各自网络中的拓扑相似特征，计算生物分子的相似矩阵S，矩阵中的每个元素为节点和节点之间的相似系数，具体步骤包括：

B1、计算生物分子和在网络拓扑结构上的平均相似性，其相似性的各个方面在有向网络中分别由 、 、 和 表示，在无向网络中则由和表示；

B11、在有向网络中，表示节点a和b的入邻居（即有边从这些邻居节点指向a或b）之间的平均相似性，表示节点a和b的出邻居（即有边从a或b指向这些邻居节点）之间的平均相似性，表示节点a和b的非入邻居节点（指那些没有边进入a和b的节点）之间的平均相似性，表示a和b的非出邻居节点（指那些没有边从a和b指向它的节点）之间的平均相似性。由表示进入节点a的边的条数，即a的入度，表示从节点a引出的边的条数，即a的出度，表示存在一条从到a的边，表示从到a不存在边，表示存在一条从a到的边，表示从a到不存在边，表示是图中的一个节点，表示是中的一个节点，表示的相似系数，为图G_t的节点数，为图G_B的节点数，上标k代表迭代次数，且

B12、在无向网络中，表示节点a和b的邻居节点（即这些节点和a或b之间存在边）之间的平均相似性，表示节点a和b的非邻居节点（即这些节点和a或b之间不存在边）之间的平均相似性，由表示节点a的边的度数，表示和a之间存在边，表示和a之间不存在边，表示是中的一个节点，表示是中的一个节点，表示的相似系数，为图G_t的节点数，为图G_B的节点数，上标k代表迭代次数，且

B2、在生物分子和的序列相似系数基础上增加它们在网络拓扑结构上的平均相似性，迭代，直至一致收敛或振荡收敛；具体公式如下，其中上标k和k+1代表迭代次数，代表和的相似系数，、 、 和 表示有向网络中和出/入邻居节点和非邻居节点之间的平均相似性，和表示无向网络和的邻居节点和非邻居节点之间的平均相似性，表示和的序列相似系数：

在有向网络中

在无向网络中 

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

C、  构建专家知识字典，字典中包含了网络T（G_t）和B（G_B）中由专家确定的最相似的生物分子对。

D、 基于专家知识，以邻居优先策略进行搜索，获得结果子网，并可视化。具体步骤包括：

D1、将目标子网T（G_t）中的生物分子分成两类，一类为K类，其中的生物分子是与具体问题相关的重要分子且存在于专家知识字典中；另一类为N类，包含了网络T（G_t）中除K类以外的其它生物分子。

D2、根据专家知识字典，在网络B（G_B）中找到与K类生物分子具有最高相似度的生物分子，使它们配对。

D3、对N类生物分子根据邻居优先的策略进行搜索配对，直到全部配对完成，获得结果子网。具体步骤包括：

D31、优先队列PQ根据节点对的相似系数维护节点对：对子网T（G_t）中的每个N类节点，都在网络B（G_B）中找出与它最相似节点，并把这些节点对加入PQ。

D32、在每一次搜索配对过程中，选择在优先队列中还没有匹配的节点中最匹配的节点对，将它们标记为“已匹配”。

D33、同时，对邻近的未匹配节点对，在其原有相似系数的基础上增加系数，从而增大它们被选择匹配的机会。

D34、这个搜索过程在子网T（G_t）中的所有节点都被匹配时结束。

E、  计算结果子网R (G_r) 与目标子网T (G_t) 的相似得分，其相似得分定义如下：

设目标子网为，结果子网为，其中，，分别代表网络，的节点集合，且₁，₂，即网络中有₁个节点，网络₂个节点；表示节点存在于网络中，，分别表示结果子网中与对应的节点；、分别代表网络、的边集合，表示边的两个端点是节点，表示边是网络的一条边；表示边的权重；表示网络的节点和结果子网中与其对应的节点的序列相似系数。则结果子网相对于目标子网的得分为

在无向图中：

其中 

在有向图中：

其中

F、  计算p值，分析目标子网的统计学意义。具体步骤如下：

F1、生成网络B（G_B）的n个随机网络。具体方法是：

F11、随机选择一对边A-B和C-D；

F12、将这两条边重新组合，使A-D相连且B-C相连；

F13、如果这些新形成的边在网络中已经有一条或几条存在，则终止这一步并重新选择一对新边进行变换。这样就可以防止两个相同的节点之间有多条边出现；

F14、重复F11-F13，形成一个与网络B（G_B）对应的随机网络。

F2、在每个随机网络中用同样的方法搜索同一个目标子网的相似子网，得到n个结果子网。

F3、用T检验计算p值。p值反映了计算结果有多大概率是由两个无关网络随机计算的结果，p值越接近于0，说明所得到的结果越显著越不可能是随机出现的结果，因此越可能具有生物学意义；反之，p值越接近于1，则所对应的结果就越不显著，越可能是由于无意义的随机计算得到的。

G、 结果子网可视化。

参照图7，示出了本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法与同类代表性方法，即NBM和MNAligner，对于网络搜索例的计算结果对比。 图7A和图7B分别是示例的两个网络及其初始相似矩阵。计算结果表明，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法和NBM方法都能找到两个网络之间较多的相似边，即，使其保持更相似的拓扑结构，而本发明的总得分高于NBM，说明在找到相似拓扑结构的同时，本发明找到的匹配节点具有更高的相似系数。可见本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法能比同类方法搜索到与目标更相似的结果子网，说明在准确性方面本方法更优。

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

参照图8，示出了对图7算例的1~7条边（约占总边数的6.7%~46.7%）各进行100次拓扑变换获得的700个不同的网络，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法与同类代表性方法NBM的计算结果对比。取G_t和G_B均为图7A中的第一个网络，即两个网络完全相同且其节点为｛A, B, C, D, E, F, G, H, I, J, K, L｝，保持G_t节点数目不变，修改边以改变网络的拓扑结构，具体为分别修改1，2，……，7条边各100次，获得700个不同拓扑结构的网络；同源表中对角线元素取该矩阵元素的最大值0.8，其余元素取自图7中的同源表；邻居优先参数ω以步长0.1取值0.1~1。图8示出了ω=0.1时，本发明和NBM对这700个G_t相对于G_B的计算结果，当ω取其它值时具有类似特性。其中：

图8A为变换G_t的1~7条边各100次后的总分平均分，实验表明，对于这700个网络，本发明的计算结果的总分平均分普遍高于NBM方法，说明本发明计算的G_r与G_t普遍更相似。

图8B为这700个网络的总分方差，图中本发明的总分方差显著低于NBM方法，说明本发明的方法具有更高稳定性。

图8C为变换G_t的1~7条边时，变换后的G_t与G_B（与原始G_t相同）边的匹配情况。实验表明，本发明的计算结果能较好地符合边的变换情况，即当变换n条边的时候，则有41%~96%的结果为n条边未匹配，而NBM的最高比例为23%~89%，且所占比例最高的网络不一定是n条边未匹配（即其峰值与变换的边数不一致，不能体现拓扑的变化情况）。

图8D为700个网络的平均边正确率（Edge Correctness，EC），即G_t与G_B中匹配的边在Gt总边数中的比例。

图8E为700个网络的EC方差，图中本发明的EC方差显著低于NBM方法，说明本发明的方法具有更高稳定性。

总之，图8的实验结果表明，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法比同类代表性方法NBM相比，具有更高的准确性和稳定性；随着生物分子网络的进化，本发明更能反应网络拓扑结构的变化情况。

参照图9，示出了各以最多100种方式删除图7算例的1~6个节点后，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法与同类代表性方法NBM的计算结果对比。取G_t和G_B均为图7A中的网络，以最多100种方式各删除G_t中的1~6个节点（占G_t总节点数的8.3%~50%），并连通删除节点的邻居节点使得网络的连通特性不变，具体为每次删除1个节点可获得12个网络，每次删除2个节点可获得66个网络，每次删除3~6个节点分别获得100个网络，如此共获得478个不同拓扑结构的网络；同源表取自图7中的同源表；邻居优先参数ω以步长0.1取值0.1~1。图9示出了ω=0.1时，本发明和NBM对这478个G_t相对于G_B的计算结果，当ω取其它值时具有类似特性。其中：

图9A为删除节点以改变拓扑结构后的总分平均分。实验表明，对于这478个网络，本发明的计算结果的总分平均分普遍高于NBM方法，说明本发明计算的G_r与G_t普遍更相似。

图9B为删除节点而改变拓扑结构后的总分平均变化率。实验表明，本发明的计算结果能较好地体现网络拓扑的变化，即总分变化趋势与网络节点的变化趋势基本一致；而NBM的计算结果则与网络拓扑的变化有较大的差异。

图9C为删除节点而改变拓扑结构后的总分方差。从图示可以看出，本发明的总分方差显著低于NBM的总分方差，说明本发明具有更高的稳定性。

图9D为删除节点而改变拓扑结构后的平均EC。从图中可以看出，本发明的平均EC普遍高于NBM，说明本发明能得到更多的保守边，即具有更高的准确性。

图9E为删除节点而改变拓扑结构后的EC方差。从图中可以看出，本发明的EC方差显著低于NBM，说明本发明具有更高的稳定性。

参照图10，示出了图7算例增加节点后，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法与同类代表性方法NBM的计算结果对比。取G_t和G_B均为图7A中的第一个网络，即两个网络完全相同且其节点为｛A, B, C, D, E, F, G, H, I, J, K, L｝，在G_t中随机增加1~2个节点（占G_t总节点数的8.3%~16.7%）形成Gap，重复7次实验；同源表中对角线元素取该矩阵元素的最大值0.8，其余元素取自图7中的同源表，新增加的节点与原来节点的同源系数均取0.8；邻居优先参数ω以步长0.1取值0.1~1。图10示出了ω=0.1时，本发明和NBM的计算结果，当ω取其它值时具有类似特性。其中：

图10A为增加节点后的总分变化；图10B为增加节点后的边保守率（即在增加节点导致网络拓扑发生变化后，原有的边还有多少仍能与自己匹配）；图10C为增加节点后的节点保守率（即在增加节点导致网络拓扑发生变化后，原有的节点还有多少仍能与自己匹配）。实验结果表明，在G_t和G_B完全一致的情况下，如果G_t增加节点发生拓扑结构的变化，本发明计算结果的总分变化比较平稳，而边保守率和节点保守率都相对较高。这说明对于不同物种生物网络存在的Gap，本发明较不容易受其影响，能较好地反映网络的保守信息。

参照图11，它是对果蝇和人类网络搜索比对中，采用或不采用专家知识的结果对比。基于上海大学生命学院分子生物学课题组的果蝇模型实验和相关文献记载，我们试图基于三个与帕金森病有密切关系的果蝇蛋白质CG7176, CG9277, CG17870来研究人类与帕金森病相关的蛋白质相互作用网络（Protein Interaction Network, PIN）。其中，果蝇蛋白质PIN数据来自于DIP数据库，包含7038个蛋白质和20720条相互作用；人类PIN数据来自于HPRD数据库，包含6340个蛋白质和23591条相互作用；目标子网即取自果蝇PIN，是由CG7176, CG9277, CG17870及与其有直接相互作用的蛋白质共同构成的网络。

图11A为本次实验的目标子网。

图11B为本发明在不采用专家知识字典的情况下的计算结果，表中列出了其中仅有的两对具有相似功能的蛋白质。

图11C为专家知识字典的内容。本实验中CG7176, CG9277, CG17870为K类蛋白，此表中列出了这些蛋白的匹配蛋白。

图11D为本发明在采用专家知识字典的情况下的计算结果（p-value=10^-19），表中列出了其中九对具有相似功能的蛋白质。

图11的实验说明，对于特定的生物问题，专家知识的引入将有效地提高结果的准确性。

综上所述，图7～图11表明，本发明的基于专家知识与拓扑相似的邻居优先生物分子子网搜索方法，与同类代表性方法相比，其总体计算准确度更高，具有更高的稳定性，能更好地处理生物分子网络间的Gap，更能反映生物分子网络之间的保守信息和差异情况。

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