CN110379457A - 一种基于整数组合的网络模体子树枚举和统计方法 - Google Patents

Abstract

本发明公开了一种基于整数组合的网络模体子树枚举和统计方法，该方法首先利用整数的组合操作设计一种有效的非导出子树的枚举方法；其次，通过在子树枚举的过程中同时搜索一个根树的方式对所枚举的子树进行一个部分的分类来有效减少子树同构的判断数目；最后，通过把所查找子树的规范化标记存储在内存中的方法来去除规范化标记的查询过程，有效地提高了子树计数的方式。本发明实现简单、快速，相比于现存的子树统计方法，MTMO方法在运行时间上具有较明显的性能加速，且能用于查找更大规模的子树。

Description

一种基于整数组合的网络模体子树枚举和统计方法

技术领域

本发明涉及一种基于整数组合的网络模体子树枚举和统计方法

背景技术

近年来，随着高通量蛋白质组技术的快速发展，可获得的蛋白质相互作用数据迅速增长。蛋白质相互作用被构建为无向图(其中顶点对应蛋白质，边对应蛋白质之间的相互作用)有助于从网络水平上理解细胞机制的基本构件和组织。蛋白质相互作用网络类似于复杂网络，拥有小世界和无尺度等全局特性。而在2002年由Milo等首次提出的网络模体则是其中非常重要的一种局部性质，网络模体定义为一种在给定网络中频繁出现的连通子图，并且其出现的次数要多于在相应的随机网络中的次数，被认为是复杂网络的基本构件块。

然而，网络模体发现涉及到图的同构判断，其运行时间随着模体规模的增大呈指数增长，且当前的模体发现方法仅能有效查找中等规模的模体。因此，如何在网络模体发现的过程中减少或者去除子图同构判断及查找出更大规模的模体将是一个挑战。

当前大部分网络模体发现方法都是以导出子图为基础，例如ESA、ESU和Kavosh。然而，值得注意的是，由于一个非树型的导出子图包含了多个相同规模大小的非导出子图，因此，在一个网络中存在更多的非导出子图同构到一个给定的拓扑结构中，使得非导出子图的发现是更具有挑战性的。此外，由于目前每一物种的蛋白质相互作用数据是不完整的，且包含了大量的假阳性数据，从而，一个网络中的某一个特定模体的一个子图出现可能在另外一个网络中其子图出现包含了附加的边，反之亦然。而非导出子图的定义考虑了子图在网络中所有可能的出现，因此，对于蛋白质网络而言，统计子图的非导出出现是一种更加健壮的方式。

鉴于任何非树型连通子图可以通过相应的树型子图进行边的扩展而得到，因此仅需要考虑子树模式的提取。然而，子树枚举和统计仍然是一个计算上具有挑战性的问题。Omidi等提出的MODA方法通过在输入网络中采用子树映射策略来查找该子树的出现频率。然而，MODA方法不能扩展到上千个顶点的大规模网络，且其性能瓶颈主要在于树型子图的映射。另外，一些研究者提出了多个以颜色编码技术为基础的串行和并行实现方法来估计网络中非导出子树和有界树宽子图出现数目。为了查找某一规模为k的特定子树t，颜色编码技术首先通过使用k种颜色来给网络中的每个顶点随机赋予一种颜色，然后通过动态规划方法来近似估计网络中与查询子树t同构的子树t′数目，其中子树t′中的每个顶点在该子树中是具有唯一的颜色。2008年，Alon等描述了如何采用颜色编码技术来统计树型子图的非导出出现数目，并使用树型子图的分布揭示了单细胞生物体与多细胞生物体之间的显著差异。随后，Zhao等结合基于流的划分和颜色编码技术实现了一个并行的子树统计方法PARSE。最近，Slota等又提出了一种基于颜色编码的多核并行子树统计方法FASCIA。

上述这些非导出子树统计方法都是属于以子树为中心，以子树为中心的方法需要预先产生一个给定大小的所有非同构子树，然后每一次在目标网络中只能映射一种特定类型的子树，导致网络的同一部分被重复映射多次，没有利用先前的子图查找信息，潜在地需要更多的时间进行重复的搜索。

目前对于树型子图的查找方法研究还处于刚起步阶段，Ferreira等在2011年提出了第一个输出敏感的方法，即该方法在规模为n的图G中枚举所有大小为k的子树所需的时间为O(sk)，其中s表示图G中所有的k-子树的数目，数量上正比于n^k。但该方法仅是理论研究，没有可用于比较的实验结果，同时，编程实现较为困难。

发明内容

本发明所要解决的技术问题是提供一种基于整数组合的网络模体子树枚举和统计方法，该方法实现简单，在运行时间上具有较明显的性能加速，且能用于查找更大规模的子树。

一种基于整数组合的网络模体子树枚举和统计方法，采用邻接边扩展的整数组合原则依次枚举网络模体对应的给定无向图G中的所有子树，同时使用根树数据结构记录每次枚举得到的子树的拓扑信息，得到所述给定无向图G的标记根树，完成对所枚举子树的部分同构分类；最后利用AHU方法得到所述标记根树中-每一个叶子结点对应子树的二进制正则标记，完成各类同构子树的数目统计。

所述标记根树中的叶子结点是指标记根树中不存在孩子结点的结点；

一个已知的网络，设定k取值后，会对应一个具体的标记根树，尽管标记根树的每一条路径会对应一个子树同构类，但是不同的路径对应的一个类子树可能是同构的，使用标记根树只是确保到达相同的路径的子树是同构的；对应不同路径的子树最后使用AHU方法产生正则标记来进行同构判断，同构子树的正则标记相同。

每一类同构子树根据中心结点的个数(1个或2个)产生1个或2个正则标记，即使是两个也没有关系，可以用来区分不同的同构子树的。

完成了子树统计后，进一步可以扩展查找非子树类型的子图，从而完成模体的查找。

进一步地，所述采用邻接边扩展的整数组合原则依次枚举网络模体对应的给定无向图G中的所有子树的过程如下：

首先，对所述给定无向图G中所有顶点从1开始进行连续整数标号；

其次，按照设定的k-子树，依次以每个顶点作为子树的第一个顶点，建立以第一个顶点为根的基于邻接关系产生的最大深度为k的隐含树；

其中，所述隐含树中每个顶点的孩子由邻接边来产生，一条邻接边的顶点能被选择作为一个顶点的孩子仅当对应的邻接顶点还没有被包含在当前所枚举的子树中；

在建立隐含树时，所有的孩子顶点的整数标号必须大于根节点的标号。

该限制实现了去除图G中先前考虑的顶点的过程。

进一步地，在所述枚举子树的过程中，先确定子树模式，然后依据子树模式进行子树枚举；

所述子树模式是指子树的层数以及每一层的节点数，子树模式的数量为所有正整数累加等于k-1形成的正整数序列的数量；

正整数序列中的第i个数的数值表示子树的第i层顶点的孩子顶点的数量。

进一步地，采用循环门排序方法进行所述正整数序列组合的选择。

进一步地，所述标记根树是在子树搜寻过程中，按照搜寻顺序对顶点进行标记编号，记录每个标记顶点的双亲顶点的标记编号，每搜寻到一个新的顶点，则产生一条比较根树的路径连接边，且以新的顶点的双亲顶点的标记编号作为根树路径中一条连接边的权值；

所述标记根树中每条路径均包含k-1条连接边。

进一步地，对所述各子树的二进制正则标记转换为十进制数，以最大十进制数和最小十进制数之差加1作为统计数据数组长度，以统计数据数组依次记录各类同构子树的同构频率，实现网络模体子树统计；

所述各类同构子树的同构频率是指各类同构子树中子树数量占所述给定无向图G中的所有子树数量的比值，从属于同一同构子树的子树二进制正则标记相同。

有益效果

本发明提供了一种基于整数组合的网络模体子树枚举和统计方法，从网络中心的角度结合整数组合的思想来枚举子树的非导出出现，同时在子树枚举的过程中通过搜索一个根树数据结构的方式对所枚举的子树进行一个部分的分类以减少子树同构判断的数目，而且使用了一种基于数组的索引方案来简化子树的计数方式。

通过利用提出的标记根树，使得到达相同叶子结点的所有子树都是同构的，即到达相同叶子结点的所有子树只需要调用一次AHU方法计算该同构类子树的正则标记。AHU方法的时间复杂度为O(k²)，而搜索标记根树的子过程的最大时间复杂度是O(k)。也就是说，对大部分枚举的子树，其判断同构的时间成本从O(k²)降到了O(k)，而对另外的一些子树(至多等于标记根树中的叶子结点数)，则需要增加额外的搜索成本O(k)。相比于大规模复杂网络中的子树的数目，标记根树中叶子结点的数目相对较少，这意味着很大一部分的子树将会到达相同的叶子结点，而每个叶子结点仅当首次被创建时才需要调用AHU方法。显然，通过搜索标记根树将会显著减少子树统计时间。此外，当子树的规模不大时，标记根树中所有的结点数目相对较少，因此构建标记根树所需要的内存开销也是可以忽略的。

附图说明

图1为枚举和统计3-子树示意图，其中，(a)为图G₁，(b)为图G₁中包含顶点1的所有3-子树，(c)为图G₁中去除顶点1包含顶点2的所有3-子树，(d)为图G₁中去除顶点1和顶点2包含顶点3的所有3-子树；

图2为本发明MTMO的流程图；

图3为一个网络实例；

图4为MTMO枚举子树实例图，其中(a)按照(1,1,1)模式建立的隐含树，(b)按照(1,2)模式建立的隐含树，(c)按照(2,1)模式建立的隐含树，(d)按照(3)模式建立的隐含树，所有图中所选择的顶点和边均通过图形加粗来标记；

图5为枚举4-子树的标记根树实例图；

图6为MTMO搜索标记根树实例图；

图7为所有5-子树的非同构类别及相应的正则标记示意图，其中，(a)度序列为(1,1,1,1,4)的子树，(b)度序列为(1,1,1,2,3)的子树，(c)度序列为(1,2,2,2,1)的子树。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明：

一种基于整数组合的网络模体子树枚举和统计方法，采用邻接边扩展的整数组合原则依次枚举网络模体对应的给定无向图G中的所有子树，同时使用根树数据结构记录每次枚举得到的子树的拓扑信息，得到所述给定无向图G的标记根树，完成对所枚举子树的部分同构分类；最后利用AHU方法得到所述标记根树中每一个叶子结点对应子树的二进制正则标记，完成各类同构子树的数目统计。

实施例1：

一种基于整数组合的网络模体子树枚举和统计方法，在后续描述中，该方法简称为MTMO；

本发明将子树统计问题定义为：给定一个规模为n的无向图G＝(V,E)和一个正整数k，查找出图G中所有可能的k-子树的频率。

一个子树T的两个出现T₁与T₂之间可能共享顶点，事实上，T₁与T₂的顶点集合可能是相同的。仅当T₁与T₂之间至少有一条不同的边时，才认为是子树T的两个不同出现。例如，如图1所示，图G₁中3-子树的频率为9。

基于整数组合的子树枚举和统计方法MTMO的整个流程如图2所示。首先输入一无向图G和一个正整数k，方法MTMO可以划分为2个子过程：子树枚举和子树分类。

1)子树枚举：

为了枚举图G中的所有非导出k-子树，首先枚举出包含一个特定顶点的所有k-子树，随后，从图G中去除该顶点，并对图G中剩下的其他顶点重复上述过程。为了去除已经考虑过的顶点和便于实现，对图G中所有顶点给予从1开始的连续整数标号。

为了查找出一个特定顶点在图G中参与的所有k-子树，以该顶点为根的基于邻接关系产生的最大深度为k的隐含树被建立。隐含树中每个顶点的孩子由邻接边来产生，一条邻接边能被选择作为一个顶点的孩子仅当对应的邻接顶点还没有被包含在当前所枚举的子树中。此外，在建立隐含树时，所有孩子顶点的整数标号必须大于根节点的标号。该限制实现了去除图G中先前考虑的顶点的过程。

子树枚举的原理充分利用了一个整数的组合操作，一个规模为k的子树包含了k-1条边，因此，为了提取k-子树，需要考虑整数k-1的所有可能的组合。在数学上，整数k-1的一个组合是由正整数组成的一个序列，该序列之和为k-1。具有相同和的两个序列若其元素顺序不同，也被看作是该和的两个不同组合。整数k-1的一个组合可以表示成：k₁,k₂,...k_m,其中k₁+k₂+...k_m＝k-1。为了基于该组合来查找k-子树，隐含树中的第i层的k_i条边被选中作为k-子树的边集。同时需要注意的是，为了确保当前枚举的子树加入一条边后其规模也加1，隐含树的每一层所选择的多条边中不能对应到相同的邻接顶点。最终，所选择的包含了k个顶点的k-1条边就定义了图G中的一个非导出k-子树。

对于一个特定的层i，可能有k_i<n_i，其中n_i是隐含树中的第i层的候选边数。那么在层i，需要考虑C(n_i,k_i)(C(n,k)表示从n个元素里面选择k个元素的不同组合数)种不同组合的边选择。本发明通过采用循环门排序(revolving door ordering)方法来实现边的所有组合选择，循环门排序方法被认为是生成组合数的最快方法，且是一个具有常量摊销时间的方法。

以图3给出的网络为例，根据子树枚举的原理，此处给出了一个简单实例，包含顶点1的所有4-子树被枚举，如图4所示。如前所述，整数k-1＝3的所有组合需要被考虑来作为不同的选择模式。整数3的所有组合有4种，分别是：(1,1,1)，(1,2)，(2,1)和(3)。以(1,1,1)组合模式为例，其中顶点1作为隐含树的根结点，按照模式的要求此时需要选择顶点1的一条邻接边，而顶点1有三条有效的候选邻接边{(1,2),(1,3),(1,4)}，通过使用循环门排序方法，边(1,2)是第一条被选择的边。同样按照模式的要求此时需要选择顶点2的一条邻接边，尽管顶点2有三条邻接边{(2,1),(2,3),(2,5)}，但是顶点1已经被选择作为当前枚举子树中的顶点，所以顶点2只有两条有效的候选邻接边{(2,3),(2,5)}，继续使用循环门排序方法，边(2,3)是第二条被选择的边。同样按照模式的要求此时需要选择顶点3的一条邻接边，尽管顶点3有三条邻接边{(3,1),(3,2),(3,6)}，但是顶点1和顶点2已经被选择作为当前枚举子树中的顶点，所以顶点3只有一条有效的候选邻接边{(3,6)}，所以边(3,6)是最后一条被选择的边，此时边集{(1,2),(2,3),(3,6)}就构成了一个大小为4的非导出子树。图4中其他的子树可以参照类似的方式得到。

子树枚举是统计子树频率的关键过程，为了枚举一个顶点参与的所有k-子树，需要考虑整数k-1的所有可能的组合，整数k-1总共有2^k-2种不同的组合，而每一种组合需要根据邻接关系选择k-1条边，且考虑到图G中每个顶点的度不会超过D(用D表示图G中所有顶点的最大度，有D≤n)，那么整个枚举过程的时间复杂度为O(2^k-2nD^k-1)≈O(2^k-2n^k)。

2)子树分类：

枚举出一个子树之后，为了统计图G中每一种非同构类的子树出现的数目，需要对该子树的同构类进行判别。传统的以网络为中心的方法需要对每一个出现的子图分别进行同构测试，然而实际的非同构类的数目是相对比较少的，特别是对于树型拓扑结构的子图。为了避免对每一个查找的子树分别进行同构测试，在子树枚举的过程中，使用一个根树数据结构来记录枚举子树的拓扑信息，从而对枚举子树进行一个部分的分类。

由于子树的枚举过程建立一个隐含树，那么，在所枚举的子树顶点集中，除了第一顶点，其他扩展的顶点都具有唯一的双亲。因此，使用标记根树来存储每一个新增顶点的双亲信息。图5给出了枚举4-子树时一个标记根树的举例，其中根树中边上的数字标记表示新增顶点连接到当前枚举子树中的第i个顶点。一般地，对于k-子树的枚举，标记的根树数据结构具有以下几个属性：

(1)从根结点到叶子结点的路径长度为k-1。

(2)假设根结点所在的层次为0。

(3)第i(i＝0,1,…k-2)层的每一个内部结点至多有i+1个孩子结点。

(4)第i(i＝0,1,…k-1)层上的结点数目至多为C(2i,i)/(i+1)，这恰好是卡特兰数列的第i项。

(5)每条边被标记为一个数字，指出当前枚举子树的新增顶点的双亲信息。

伴随着当前枚举子树被边扩展的同时，标记根树可以被搜索。由于枚举子树扩展一条边的同时会新增一个顶点，那么可以使用这个新增顶点的双亲信息去搜索标记根树。需要注意的是，初始时标记根树的指针设置为根结点，在每一次搜索的过程中，标记根树要么创建一个新的分枝，要么沿着一条已经存在的路径。图6给出了一个子树在枚举过程中搜索标记根树的实例。

当枚举子树的规模到达了k，标记根树的指针将会移到合适的叶子结点上。标记根树中从根结点到任何叶子结点的一条路径对应于一个特定子树类型的一种顶点连接模式，因此，到达相同叶子结点的所有子树都是同构的，然而，两个同构的子树可能到达两个不同的叶子结点，即，可能存在多个不同的叶子结点对应于相同的子树同构类。表1给出了不同k-子树非同构类数目以及相应标记根树的最大叶子结点数目(为第k-1个卡特兰数)。因此，对于每一个叶子结点，需要计算一次其对应子树的正则标记以确定该子树的同构类别。

表1 不同k-子树非同构类数目以及相应标记根树的最大叶子结点数目

对于一个普通的k-树，使用AHU方法可以得到一个或两个长度为2k的二进制正则标记。图7给出了所有5-子树的非同构类别及相应的正则标记。有了子树的正则标记后，还需要使用该正则标记搜索一个二叉树并返回相应的叶子结点，该叶子结点就对应于一个特定的非同构类别。

本发明为了简化子树的统计方法，将二进制的正则标记转换为十进制数字，并将十进制数字对应于一个数组的索引。从而，每一种非同构的k-子树的频率被存储在相应的数组元素中。然而，该数组中大部分的元素没有被使用，导致内存空间的浪费。考虑到对于一个给定规模的所有非同构子树类别，星型树的正则标记是最小的，而路径树的正则标记最大，且最大值和最小值之间的差异不大。因此，为了节约内存开销，使用一个更小规模的数组，其长度为最大值与最小值之差加1。其中，数组的第一元素存放的是星型树的频率，而数组的最后一个元素存放的是路径树的频率。从而，图G中每一种非同构的k-子树的频率被存储在相应的数组元素中，整个子树枚举和统计过程结束。表2给出了不同k-子树正则标记的最小值、最大值以及相应的数组长度。

表2 不同k-子树正则标记的最小值、最大值以及相应的数组长度

基于整数组合的子树枚举和统计方法有效性验证

为了验证方法MTMO的有效性，发明人将MTMO方法应用于三个有代表性的网络，其中一个是生物网络：来自芽殖酵母的蛋白质相互作用网络，另外两个是非生物网络：一个电子网络和一个海豚的社会网络。网络中所有顶点的自身环被移除，表3给出了这三个网络的基本拓扑属性。将方法MTMO与当前经典的子树统计方法MODA进行了比较和分析；此外，在本实施例中，还测试了标记根树对MTMO方法性能的影响。

表3 实验数据集

除非另有说明，本实施例使用的计算平台的配置为Intel Xeon X5670CPU，主频2.93GHz，48GB内存，操作系统为Ubuntu 12.04。MTMO方法采用C++语言编程实现，编译源代码的编译器版本为GNU gcc/g++4.6.3，并使用“-O3”选项。

1.方法MTMO与MODA方法的性能比较

MODA方法是一个最近开源的网络模体发现工具并用于非导出子图的计数。该软件运行于Windows平台下的Microsoft Visual Studio。因此，为了获得更为直接的比较，发明人同样在Windows平台下采用C++语言编程实现了MTMO方法，两个方法均运行于一台服务器，配置为Intel Xeon X5670CPU，主频2.93GHz，48GB内存，操作系统为Windows 7。考虑到MODA方法无法适用于上千个顶点的大规模网络，此处仅使用其他两个非生物网络进行比较。表4给出了两方法的实验比较结果，从表中可以看出，相比于MODA，MTMO获得了显著性的加速，且随着子树规模的增大，优势更加明显。此外，MTMO具有更小的内存开销，能够统计出更大规模的子树。

表4 MTMO与MODA实验结果的比较

MEM：表示内存不足

2.标记根树对性能的影响

为了验证标记根树的加速效果，另外发明人实现了未使用标记根树搜索的方法版本，记为NTMO，即NTMO方法对于每一个枚举子树都要进行一次同构测试。表5比较了MTMO和NTMO两个方法在上述三个网络中进行k-子树(3≤k≤12)统计所需要的时间，由于计算资源的限制，仅给出了能在一周时间内运行完成的实验数据。

表5 MTMO与NTMO实验结果的比较

-：表示一周时间内未运行完

从表5可以看出，MTMO方法在任一网络和任何子树规模下都要优于NTMO方法，并且达到了一个数量级的加速比。而且，随着子树规模的增大，加速比也在显著增加。更加准确地说，加速比与网络中出现子树的数目和这些子树的非同构类的数目的比例有关。例如，当在酵母蛋白质网络中统计所有可能的7-子树时，加速比达到了70多，但是在电子网络中其加速比的大小仅为34。这主要是由于在子树非同构类的数目相同的情况下，酵母网络中出现的子树数目远大于电子网络中出现的子树数目的缘故。

实验证明，本发明相比其他方法在运行时间上具有较明显的性能加速，且能用于查找更大规模的子树，从而使得在大规模的网络中统计更大规模的子树成为可能。

Claims (6)

1.一种基于整数组合的网络模体子树枚举和统计方法，其特征在于，采用邻接边扩展的整数组合原则依次枚举网络模体对应的给定无向图G中的所有子树，同时使用根树数据结构记录每次枚举得到的子树的拓扑信息，得到所述给定无向图G的标记根树，完成对所枚举子树的部分同构分类；最后利用AHU方法得到所述标记根树中-每一个叶子结点对应子树的二进制正则标记，完成各类同构子树的数目统计。

2.根据权利要求1所述的方法，其特征在于，所述采用邻接边扩展的整数组合原则依次枚举网络模体对应的给定无向图G中的所有子树的过程如下：

首先，对所述给定无向图G中所有顶点从1开始进行连续整数标号；

其次，按照设定的k-子树，依次以每个顶点作为子树的第一个顶点，建立以第一个顶点为根的基于邻接关系产生的最大深度为k的隐含树；

其中，所述隐含树中每个顶点的孩子由邻接边来产生，一条邻接边的顶点能被选择作为一个顶点的孩子仅当对应的邻接顶点还没有被包含在当前所枚举的子树中；

在建立隐含树时，所有的孩子顶点的整数标号必须大于根节点的标号。

3.根据权利要求1或2所述的方法，其特征在于，在所述枚举子树的过程中，先确定子树模式，然后依据子树模式进行子树枚举；

所述子树模式是指子树的层数以及每一层的节点数，子树模式的数量为所有正整数累加等于k-1形成的正整数序列的数量；

正整数序列中的第i个数的数值表示子树的第i层顶点的孩子顶点的数量。

4.根据权利要求3所述的方法，其特征在于，采用循环门排序方法进行所述正整数序列组合的选择。

5.根据权利要求3所述的方法，其特征在于，所述标记根树是在子树搜寻过程中，按照搜寻顺序对顶点进行标记编号，记录每个标记顶点的双亲顶点的标记编号，每搜寻到一个新的顶点，则产生一条比较根树的路径连接边，且以新的顶点的双亲顶点的标记编号作为根树路径中一条连接边的权值；

所述标记根树中每条路径均包含k-1条连接边。

6.根据权利要求3所述的方法，其特征在于，对所述各子树的二进制正则标记转换为十进制数，以最大十进制数和最小十进制数之差加1作为统计数据数组长度，以统计数据数组依次记录各类同构子树的同构频率，实现网络模体子树统计；

所述各类同构子树的同构频率是指各类同构子树中子树数量占所述给定无向图G中的所有子树数量的比值，从属于同一同构子树的子树二进制正则标记相同。