CN102779241B

CN102779241B - 基于人工蜂群繁殖机制的ppi网络聚类方法

Info

Publication number: CN102779241B
Application number: CN201210232874.0A
Authority: CN
Inventors: 雷秀娟; 李永明; 田建芳; 裘国永; 吴爽; 尤梦丽
Original assignee: Shaanxi Normal University
Current assignee: Shaanxi Normal University
Priority date: 2012-07-06
Filing date: 2012-07-06
Publication date: 2015-04-22
Anticipated expiration: 2032-07-06
Also published as: CN102779241A

Abstract

本发明公开了一种基于人工蜂群繁殖机制的PPI网络聚类方法，具体包括如下步骤：将PPI网络转化为无向加权图；参数设置；对PPI网络的每个结点和边预处理；计算所有结点的加权网络综合特征值；初始化蜂后；婚飞过程；幼蜂的局部搜索；蜂后的选优；计算当前适应度并比较，直到输出全局最优聚类结果。该方法不需要预先设定聚类个数，在聚类过程中能够自动得到，避免了人为设定聚类个数的主观性，且时间复杂度降低明显。采用MIPS数据库做实验仿真，结果比较接近标准数据库且在正确率、查全率和运行时间等指标上性能较优。与其他的聚类方法相比，该方法采用基于繁殖机制人工蜂群方法机理，能自动确定聚类个数，实现聚类过程，有效提高了聚类效果和计算效率。

Description

基于人工蜂群繁殖机制的PPI网络聚类方法

技术领域

本发明属于生物信息计算领域，具体涉及蛋白质交互作用（Protein-proteininteraction，PPI）网络的自动获得聚类个数的聚类方法，PPI网络具有小世界、无尺度特性，本发明可推广并适用于其他小世界、无尺度网络的聚类问题。

背景技术

现有的聚类方法有很多种，主要包括：基于划分的方法、基于密度的方法、基于网络的方法、基于模型的方法、基于层次的方法、模糊聚类方法、谱聚类方法、功能流模拟方法、整体聚类方法等等。但是这些方法要么对应用领域和数据特性有特殊要求，要么方法本身存在一些缺陷，有些不适用于蛋白质相互作用（PPI）网络，有些对PPI网络的聚类效果不理想。

基于划分的方法以k_means和k_mediod为代表，对于大数据集该方法具有相对可伸缩性和高效性，但是该方法最大缺点在于必须事先知道要划分的目标类k的确定个数，且大都根据对象之间的距离进行聚类，只能发现球状的簇，不适合发现非凸面形状的簇，对于发现任意形状的簇遇到了困难，并且对孤立点是极其敏感的，但PPI网络恰恰形状不规则，有大量孤立点存在。

基于密度的方法将簇看作是数据空间中被低密度区域分隔开的稠密对象区域。该方法的主要思想是：只要邻域中的密度（对象或数据点的数目）超过某个阈值，就继续聚类，直至密度边缘。但是该方法的缺陷在于不能对存在大量稀疏结点的网络进行分类，因此对于存在有大量稀疏点PPI网络不适用。

基于网络的方法将空间划分为有矩形网格单元组成的网络，然后在网格单元上进行聚类。但是该方法的缺陷在于输入参数对聚类结果的影响很大，且这些参数比较难设置。当数据中有噪音时，如果不加特殊处理，方法的聚类质量会很差，并且该方法对于数据维度的可伸缩性较差。

基于模型的方法基本思想是：为每个聚类假设一个模型，再去发现符合模型的数据集对象，试图将给定的数据与某个数学模型达成最佳拟合，该方法可以通过构建反映数据点空间分布的密度函数来定位聚类，但是该方法的主要缺点在于对于要聚类的海量数据很难建立合适的、通用的模型。

基于层次的方法是根据给定的簇间距离度量准则，构造和维护一颗由簇和子簇形成的聚类树，直到满足某个终止条件为止。层次聚类方法的缺陷在于它对非球形簇和尺寸非均匀簇分裂以及孤立点非常敏感，该方法对噪声数据的鲁棒性比较差。

模糊聚类方法的主要思想是把需要识别的实物与模板进行模糊比对，从而得到所属的类别，其缺陷在于：聚类个数需要预先给定而且没有准则可遵循；只能识别团状的簇，不能识别不规则簇和带状簇，很多情况下对噪声点很敏感；对初始聚类中心敏感，容易陷入局部最优，难以取得全局最优或者整个聚类过程需要很长时间才能收敛到全局最优，从而影响聚类效果。

谱聚类方法本质的思想是将聚类问题转化成图的最优划分问题。谱聚类方法的大致流程是用将要聚类的数据样本构造一个相似度矩阵，计算此矩阵的特征值和特征向量，根据选用不同的特征向量，找出数据之间的内在联系，以达到聚类不同的数据目的。对特征向量采用不同的处理方法，就会产生不同的聚类方法。谱聚类方法虽然能在任意形状的样本空间上聚类，但对于海量数据聚类其邻域矩阵的选取和构造始终是一难题。

功能流模拟方法（Flow方法）比较符合蛋白质相互作用的实际效果，是目前国际上较为有效的一种PPI网络聚类方法。该方法采用了网络流的思想，这种方法与图割的概念有些相似。其思想是将每一个已功能注释过的蛋白质作为功能流方法中的一个“源点”，将自己的信息通过网络中的边“流”向其相邻的结点。这种流动方式有一些简单规则进行控制。该方法的缺陷在于时间复杂度较高，需要提前确定聚类个数，没有考虑到结点间距离的作用效果，而且还需要人为设定合并阈值。

整体聚类方法使用网络拓扑结构的距离度量对蛋白质网络进行初始聚类，然后使用基本聚类方法获得一系列的聚类结果，最后在该聚类集合上使用整体聚类方法得到最终的聚类结果。整体聚类方法的目标就是把多个独立不同的聚类融合成为单一的全面聚类，从而提高对无尺度网络聚类的质量。但此方法缺乏全局目标函数，每一步都需要局部地确定需要合并的聚类。此外，该方法的时间复杂度较高。

综上，上述聚类方法大都需要事先确定聚类个数，没有准则可遵循，因此是非常主观的，且这些方法的适应性普遍较差，对大型蛋白质相互作用网络都不能取得良好的聚类效果。

发明内容

针对上述现有聚类方法的缺陷或不足，本发明要解决的技术问题在于克服小世界、无尺度网络功能模块聚类问题时事先人为设定聚类个数的主观性，该方法采用基于繁殖机制人工蜂群方法机理，能自动确定聚类个数，实现聚类过程，有效提高了聚类效果和计算效率。

为了达到上述目的，本发明采用如下的技术解决方案：

一种基于人工蜂群繁殖机制的PPI网络聚类方法，具体包括如下步骤：

1）将PPI网络转化为一个无向加权图：

将PPI网络转化成一个无向加权图G（V，E），其中，V={v_i,i=1,2,…,n}为结点v_i的集合，E为边e的集合，结点v_i表示蛋白质，边e表示蛋白质之间的相互作用，w_ij表示结点v_i和结点v_j之间相互作用的大小，也就是结点v_i和结点v_j之间边e_ij的权值，若v_i和v_j之间没有边则w_ij=0，i＝1,2,…,n，j＝1,2,…,n。

2）参数设置：

令count，maxcount分别表示控制外循环的当前迭代次数和外循环对应的最大迭代次数，maxcount∈[10,1000]，并令count=1；iter，maxiter分别表示控制内循环的当前迭代次数和内循环对应的最大迭代次数，maxiter∈[10,200]；N和S分别表示蜂后婚飞时的能量和速度，N∈[50,1000]，S∈[10,500]；fval，gfval分别表示当前适应度和全局最优适应度，令gfval＝∞；cluster，gcluster分别表示当前的聚类结果和全局最优聚类结果；visited是结点被访问标记；n表示蛋白质结点个数，T表示蜂后与雄蜂交配成功的计数器；

3）对PPI网络的每个结点v_i和每个边e进行预处理：

所述对结点v_i预处理是计算结点v_i的加权网络综合特征值com-value_i；所述对边e预处理是计算改进的边聚集系数CC_i，j；

4）计算所有结点的加权网络综合特征值com-value_i的代数平均值Av-com-value，将大于Av-com-value*W的结点保存，其中W∈[0.5,2]；

5）初始化蜂后：

令iter=1，蜂后代表聚类中心，初始化蜂后就是确定第一个聚类中心，从大于Av-com-value*W的结点中随机选取一个结点作为第一个蜂后，并令该蜂后结点的visited＝1；

6）婚飞过程：

给蜂后的能量和速度赋初值，令N＝100，S＝50，令蜂后与雄蜂交配成功的计数器T＝0，将与蜂后结点的改进的边的聚集系数CC_i,j大于零的结点作为要与蜂后交配的雄蜂，将所有雄蜂结点按照该雄蜂结点与蜂后结点的改进的边的聚集系数CC_i,j降序排列，排序后的雄蜂依次与蜂后交配，每交配成功一次，蜂后的速度S和能量N就以式10和式11衰减一次；通过式9计算每个雄蜂结点与蜂后结点的交配成功概率P(i)，同时生成一个[0,1]之间的随机数rand，若P(i)>rand，则交配成功，将该雄蜂的精子加入蜂后的受精囊中，并将该雄蜂结点的访问标记visited修改为1，计数器T＝T+1；交配失败则继续与下一个雄蜂交配，直至蜂后的能量N小于能量阈值Thred∈[0.01,0.99]或者蜂后的受精囊中的精子数量大于M∈[20,200]，一次婚飞过程结束；进入幼蜂的局部搜索过程；

P(i)＝exp[-Δ(f)/S(t)]

式9

其中，P(i)表示雄蜂结点v_i和蜂后交配成功的概率；-Δ(f)表示雄蜂结点v_i和蜂后结点加权网络综合特征值com-value_i的差值；S(t)是蜂后在时刻t的速度S；

蜂后的速度S和能量E以式10和式11的方式衰减：

S(t+1)＝α*S(t) 式10

N(t+1)＝N(t)-β*N(t)/M 式11

其中，β*N(t)/M是每次转移后能量的消耗量；β∈[0.5,2]为能量衰减因子；E(t)表示当前蜂后的能量，E(t+1)表示与一个雄蜂交配后的蜂后的能量；S(t+1)表示与一个雄蜂交配后的蜂后的速度；M表示受精囊大小，α∈[0,1]为每次速度的衰减因子；

7）幼蜂的局部搜索过程：

将蜂后受精囊中每个精子结点的邻接点中结点加权网络综合特征值com-value_i最大的结点保存下来，作为发育优良的幼蜂结点；

8）蜂后的选优过程：

从发育优良的幼蜂结点中选取结点加权网络综合特征值com-value_i最大的结点作为新的蜂后结点，从而更新了聚类中心；

9）iter=iter+1，如果iter<=maxiter，并返回步骤6)，否则，转向步骤10)；

10）计算当前适应度fval，若当前适应度fval＜gfval，则令gcluster=cluster，gfval＝fval；否则gcluster和gfval不变；

11）count=count+1，如果count<=maxcount，并令所有结点的访问标记visited＝0，并返回步骤4），否则，输出全局最优聚类结果gcluster。

进一步的，步骤3）中所述的结点v_i的加权网络综合特征值com-value_i的计算过程如下：

结点v_i的加权度dw_i为与结点v_i连接的所有边的权值之和：

{dw}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{ij}

式1

结点v_i的聚集度dk_i为与结点v_i连接的近邻结点之间的连接的边的个数，定义如下：

dk_i＝|{(v_j,v_k)|(v_i,v_j)∈E，(v_i,v_k)∈E，(v_j,v_k)∈E，v_i,v_j,v_k∈E}| 式2

结点v_i的加权聚集度wk_i为：

{wk}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{jk}

式3

其中，N_j＝{v_j|(v_i,v_j)∈E,(v_i,v_k)∈E,(v_j,v_k)∈E,v_i,v_j,v_k∈E}，j∈N_j，N_k＝{v_k|(v_i,v_j)∈E,(v_i,v_k)∈E,(v_j,v_k)∈E,v_i,v_j,v_k∈E}，k∈N_k。

结点v_i的聚类系数coe_i反映了结点局部范围内的结点之间的交互作用：

{coe}_{i} = \frac{2 * {wk}_{i}}{d_{i} * (d_{i} - 1)}

式4

结点v_i的加权网络综合特征值com-value_i：

com-value_i＝α*coe_i+(1-α)*dw_i/n 式

5

其中，i,j＝1,2,......,n；(v_i,v_j)∈E表示结点v_i和vj之间有相互作用；d_i表示与结点v_i连接的结点连接成三角形的个数；n为蛋白质结点的总个数；α是[0,1]之间的一个可调节参数，一般取0.5；

进一步的，步骤3）中所述的改进的边聚集系数CC_i,j的计算过程如下：

基于度的边的聚集系数C_i,j定义为：

C_{i, j} = \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)}

式6

其中，N_i，N_j分别表示结点v_i和v_j的邻接点的集合，边的聚集系数C_i,j是一个局部变量，表示边e(i,j)的两个结点v_i和v_j的亲疏程度；C_i,j的取值为[0,1]；

加权聚集系数WC_u，v定义如下：

{WC}_{i, j} = \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式7

改进的边聚集系数CC_i,j定义为：

{CC}_{i, j} = ω * \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)} + (1 - ω) * \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式8

其中，ω是[0,1]之间的一个可调节参数，取0.5；N_i，N_j分别表示结点v_i和v_j的邻接点的集合，I_i,j表示结点v_i和v_j的邻接点的公共邻接点集合。

进一步的，步骤10）中所述的当前适应度fval的计算过程如下：

PPI网络的子图密度度量如下：

D (t) = \frac{2 c_{t}}{n_{t} * (n_{t - 1})}

式12

式中，t＝1,2,…,numclu，n_t是第t类的结点个数，c_t是PPI网络第t类的边的个数，D(t)表示第t类的密度；

两个功能模块M_s和M_t间的相似度矩阵S(M_s,M_t)通过式13计算：

S (M_{s}, M_{t}) = \frac{Σ_{x &Element; M_{s}}, y &Element; M_{t} c (x, y)}{\min (| M_{s} |, | M_{t} |)}

式13

其中，

c (x, y) = \{\begin{matrix} 1 & ifx = y \\ w_{e} (x, y) & ifx &NotEqual; yand < x, y > &Element; E \\ 0 & otherwise \end{matrix};

当前适应度fval定义如式14：

fval = \frac{1}{numclu} Σ_{t = 1}^{numclu} D (t) + \frac{Σ_{s = 1}^{numclu} Σ_{t = 1}^{numclu} S (M_{s}, M_{t})}{numclu * numclu}

式14

其中，S(M_s,M_t)为相似度矩阵，D(t)表示第t类的密度，numclu为聚类个数。

进一步的，所述步骤2）中设置蜂后婚飞时的能量E=100，蜂后婚飞时的速度S＝50，能量阈值Thred＝0.5，蜂后受精囊大小M＝50。

本发明与现有的方法相比，具有以下优点：

（1）自动确定聚类个数。本发明能够根据自身的特点和PPI网络的特性，自动确定聚类个数，避免现有方法中人为设定聚类个数的主观性，且聚类效果比较理想，经实际数据验证，利用本发明的方法确定的聚类个数比较接近标准数据库。

（2）该方法的聚类效果理想。本方法根据蜂群中三种不同功能的蜜蜂实现，模拟人工蜂群繁殖机制，其本身就具有优化的作用，因此聚类效果比较理想。

（3）识别蛋白质功能。利用本发明的聚类结果能够有效识别具有相同功能的蛋白质，从而能够使得研究人员有效识别功能未知的蛋白质功能模块，了解特定生物功能的分子机制，为探讨重大疾病的机理、疾病治疗、疾病预防和新药开发提供重要的理论基础。本发明能够从大量的蛋白质相互作用网络数据中提取出功能模块，功能模块是指在不同的时间和空间通过相互绑定的方式参与某一个特定的生化过程或者细胞功能，最终这些蛋白质复合物和功能模块再以层次嵌套的方式构成PPI网络，即蛋白质相互作用网络的聚类，该研究是生物体行为理解、蛋白质功能预测和药物设计的基础。

附图说明

图1是本发明的方法的机理模型图。

图2是本发明实施例的流程图。

以下结合附图和实施例子对本发明进一步详细说明，但本发明不限于这些实施例子。

具体实施方式

研究蛋白网络（PPI网络）的方法具有一定的普适意义，在蛋白网络中被成功应用的算法可应用到当前网络科学研究其它领域，如互联网、人际关系网和生物代谢网等，处理类似的问题，对其它领域的研究有一定的借鉴意义。

由于PPI网络的小世界特性和无尺度特性，PPI网络的聚类个数未知，因此本发明定义了基于人工蜂群繁殖机制的PPI网络聚类方法（The Clustering MethodBased on Marriage Behavior in Artificial Bees Colony，CM-MBABC），该方法根据PPI网络特性和方法本身的特点自动确定聚类个数。

（1）CM-MBABC方法的设计机理

基于繁殖机制的人工蜂群优化方法的启发，我们将蜂后、雄蜂和幼蜂三种蜜蜂分别对应到PPI网络的聚类问题中：蜂后对应聚类中心结点；雄蜂对应与聚类中心结点（蜂后结点）的改进的边的聚类系数大于零的结点；幼蜂对应候选蜂后。本发明的方法中，一个蜜蜂群体中只有一个蜂后，一个类也只有一个聚类中心。

蜂后繁殖过程包括三个过程：蜂后与雄蜂交配的过程；幼蜂的局部搜索过程；蜂后的选优过程；该三个过程对应PPI网络的聚类过程，它们的机理如下：

蜂后与雄蜂交配的过程，即婚飞过程（参见图1、图2）：首先，初始化蜂后：蜂后的初始化过程就是选取第一个聚类中心的过程，由于PPI网络的特性以及经验可知聚类中心倾向于结点信息较大的结点，但是又不能盲目地选取结点信息最大的结点，因此计算所有结点加权网络综合特征值的平均值，从大于结点加权网络综合特征值平均值的结点中随机选择一个结点作为第一个蜂后；其次，将蜂后结点的每个邻接点作为雄蜂结点依次与蜂后交配，若交配成功，将交配成功的雄蜂结点的精子加入蜂后的受精囊中（本发明中是对蛋白质结点聚类，因此本发明中雄蜂的精子结点就是雄蜂结点本身），每交配一次，蜂后的能量和速度衰减一次，直到蜂后的能量小于能量阈值或者蜂后的受精囊已满，交配停止即一次婚飞结束。蜂后的空间运动决定了蜂后婚飞的速度，蜂后婚飞能量和受精囊M的大小的决定了蜂后一次婚飞中与雄蜂交配的次数，也决定了一个类中结点个数的最大值；

幼蜂的局部搜索：将蜂后受精囊中每个精子结点的邻接点中结点加权网络综合特征值最大的结点作为发育优良的幼蜂。

蜂后的选优过程：从发育优良的幼蜂中选取结点加权网络综合特征值最大的结点作为新的蜂后。利用蜂后的选优过程来优化和更新聚类中心，得到新的蜂后。

综上，本发明的方法模仿人工蜂群繁殖机制，将三个过程即完整的聚类过程反复循环实现，最终找到最优的聚类结果。

（2）本发明方法的技术路线与实施步骤

CM-MBABC方法中蜂后与雄蜂的交配过程中，雄蜂的选取是只要某结点与蜂后结点有交互作用，该结点就可以看作是一只要与蜂后交配的雄蜂。但是结点的交互作用只反映两个结点之间的交互强度，不能反映结点局部范围内的连接强度和密度，而边的聚集系数能够很好地反映结点局部范围内的连接强度和密度以及两个结点属于同一类的可能性大小，因此利用边的聚集系数研究PPI网络聚类是非常有意义的。边的聚集系数分为度的聚集系数和加权聚集系数，结点度反映结点间的连接紧密程度，结点的加权度反映了结点间交互作用的大小。

为了综合考虑结点间的连接紧密程度和交互作用大小，本发明将边的聚集系数引入到本发明的方法的数据预处理过程中，将与蜂后结点的边的聚集系数大于零的结点按照降序排序作为一组雄蜂，依次与蜂后进行交配。根据边的聚集系数可有效地将一些与蜂后结点连接稀疏且交互作用小的结点剔除掉，避免噪声点的干扰。同时将结点的度、加权度、聚类系数和结点网络综合特征值四种不同的预处理方式用于蜂后的初始化。

如图2所示，本发明的基于人工蜂群繁殖机制的PPI网络聚类方法，具体包括如下步骤：

1）将PPI网络（即蛋白质相互作用网络）转化为一个无向加权图。

2）参数设置。

3）对PPI网络的每个结点v_i和每个边e进行预处理。

对结点v_i预处理是计算结点v_i的加权网络综合特征值com-value_i；对边e预处理是计算改进的边聚集系数CC_i，j。

结点v_i的加权度dw_i为与结点v_i连接的所有边的权值之和：

{dw}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{ij}

式1

结点v_i的加权聚集度为：

{wk}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{jk}

式3

{coe}_{i} = \frac{2 * {wk}_{i}}{d_{i} * (d_{i} - 1)}

式4

结点v_i的加权网络综合特征值com-value_i：

com-value_i＝α*coe_i+(1-α)*dw_i/n 式

5

给定一个蛋白质网络，若要从中有效识别功能模块，直觉上只要找到所有的连接不同功能模块的边即可。聚集系数用于表示网络中某结点的邻接点之间的亲疏程度。基于度的边的聚集系数C_i,j定义为：

C_{i, j} = \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)}

式6

其中，N_i，N_j分别表示结点v_i和v_j的邻接点的集合，边的聚集系数C_i,j是一个局部变量，它表示边e(i,j)的两个结点v_i和v_j的亲疏程度。C_i,j的取值为[0,1]，其取值越大表明结点v_i和v_j属于同一个功能模块的可能性越大。

为降低PPI网络中假阳性对聚类算法的影响，将边的聚集系数扩展到加权网络，Radicchi给出加权网络中的定义叫做加权聚集系数WC_u，v，定义如下：

{WC}_{i, j} = \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式7

从式6看出，基于度的边的聚集系数C_i，j是边的两个结点的公共邻接点加1与这两个结点的邻接点的最小值的比值。式7中加权聚集系数WC_i，j用两个结点各自与公共邻接点权值之和的乘积与两个结点和其邻接点权值之和的乘积的比值来刻画。另外，结点v_i和v_j之间的度反映了结点之间的连接强度，结点的加权度反映了结点v_i和v_j之间交互作用的大小，为了更准确地刻画v_i和v_j两个结点属于同一个功能模块的可能性大小，本发明将上述边的聚集系数进行了改进，定义改进的边聚集系数CC_i，j为：

{CC}_{i, j} = ω * \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)} + (1 - ω) * \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式8

5）初始化蜂后。

6）婚飞过程（聚类过程）。

给蜂后的能量和速度赋初值，令N＝100，S＝50，令蜂后与雄蜂交配成功的计数器T＝0，将与蜂后结点的改进的边的聚集系数CC_i,j大于零的结点作为要与蜂后交配的雄蜂，将所有雄蜂结点按照该雄蜂结点与蜂后结点的改进的边的聚集系数CC_i,j降序排列，排序后的雄蜂依次与蜂后交配，每交配成功一次，蜂后的速度S和能量N就以式10和式11衰减一次；通过式9计算每个雄蜂结点与蜂后结点的交配成功概率P(i)，同时生成一个[0,1]之间的随机数rand，若P(i)>rand，则交配成功，将该雄蜂的精子（雄蜂结点本身）加入蜂后的受精囊中（即归入类），并将该雄蜂结点的访问标记visited修改为1（即该雄蜂结点已经被访问），计数器T＝T+1；交配失败则继续与下一个雄蜂交配，直至蜂后的能量N小于能量阈值Thred∈[0.01,0.99]或者蜂后的受精囊中的精子数量（即计数器T的值）大于M∈[20,200]，一次婚飞过程结束；进入幼蜂的局部搜索过程；

P(i)＝exp[-Δ(f)/S(t)]

式9

其中，P(i)表示雄蜂结点v_i和蜂后交配成功的概率；-Δ(f)表示雄蜂结点v_i和蜂后结点加权网络综合特征值com-value_i的差值；S(t)是蜂后在时刻t的速度S；可以看出，蜂后开始婚飞后，速度较大时或者雄蜂和蜂后的加权网络综合特征值com-value_i相等时，交配的概率较大。随着时间的推移，蜂后的速度S和能量E以式10和式11的方式衰减：

S(t+1)＝α*S(t) 式10

N(t+1)＝N(t)-β*N(t)/M 式11

式中，β*N(t)/M是每次转移后能量的消耗量；β∈[0.5,2]为能量衰减因子；E(t)表示当前蜂后的能量，E(t+1)表示与一个雄蜂交配后的蜂后的能量；S(t+1)表示与一个雄蜂交配后的蜂后的速度；M表示受精囊大小（即类中可容纳结点的最大数量），α∈[0,1]为每次速度的衰减因子。

7）幼蜂的局部搜索过程

将蜂后受精囊中每个精子结点的邻接点中结点加权网络综合特征值com-value_i最大的结点保存下来，作为发育优良的幼蜂结点。

8）蜂后的选优过程

从发育优良的幼蜂结点中选取结点加权网络综合特征值com-value_i最大的结点作为新的蜂后结点，从而更新了聚类中心。

10）计算当前适应度fval；

PPI网络的子图密度度量如下：

D (t) = \frac{2 c_{t}}{n_{t} * (n_{t - 1})}

式12

式中t＝1,2,…,numclu，n_t是蛋白质个数（表示第t类的结点个数），c_t是PPI网络第t类的交互作用个数，D(t)表示第t类的密度；

两个功能模块M_s和M_t间的相似度矩阵S(M_s,M_t)通过式13计算：

S (M_{s}, M_{t}) = \frac{Σ_{x &Element; M_{s}}, y &Element; M_{t} c (x, y)}{\min (| M_{s} |, | M_{t} |)}

式13

其中，

c (x, y) = \{\begin{matrix} 1 & ifx = y \\ w_{e} (x, y) & ifx &NotEqual; yand < x, y > &Element; E \\ 0 & otherwise \end{matrix} .

当前适应度fval定义如式14：

fval = \frac{1}{numclu} Σ_{t = 1}^{numclu} D (t) + \frac{Σ_{s = 1}^{numclu} Σ_{t = 1}^{numclu} S (M_{s}, M_{t})}{numclu * numclu}

式14

其中，S(M_s,M_t)为相似度矩阵，D(t)表示第t类的密度，numclu为聚类个数；若当前适应度fval＜gfval，则令gcluster=cluster，gfval＝fval；否则gcluster和gfval不变；

11）count=count+1，如果count<=maxcount，并令所有结点的访问标记visited=0，并返回步骤4），否则，输出全局最优聚类结果gcluster。

实施例

本实施例以采自MIPS数据库的PPI数据集作为仿真数据集，该数据集提供了关于开放阅读框架、RNA基因和其它遗传因素的信息。实验平台为WindowsXP操作系统，Intel酷睿2双核2.0GHz处理器，2GB物理内存，用Matlab 7.7软件实现本发明的CM_MBABC方法。具体步骤如下：

(1)将PPI网络转化为一个无向加权图；

(2)进行参数设置：令count，maxcount分别表示控制外循环的当前迭代次数和外循环对应的最大迭代次数，令maxcount=100；iter，maxiter分别表示控制内循环的当前迭代次数和内循环对应的最大迭代次数，令maxiter=150；E和S分别表示蜂后婚飞时的能量和速度，令E=100，S＝50，能量阈值Thred＝0.5，蜂后受精囊大小M＝50；fval，gfval分别表示当前适应度和全局最优适应度，令gfval＝inf，cluster，gcluster分别表示当前的聚类结果和全局最优聚类结果；visited是结点的被访问标记；

(3)对PPI网络进行预处理：即计算结点加权网络综合特征值com-value_i和改进的边聚集系数CC_u,v。

(4)计算得到所有结点网络综合特征值的平均值Av-com-value，将大于Av-com-value*W的结点保存起来，作为初始化蜂后的候选结点，其中令W＝1.5，count＝1；

(5)初始化蜂后：令iter=1，蜂后代表聚类中心，初始化蜂后就是确定第一个聚类中心，从大于Av-com-value*W的结点中随机选取一个结点作为第一个蜂后，并令蜂后结点的visited＝1。

(6)婚飞过程：将与蜂后结点的改进的边的聚集系数CC_u,v大于零的结点作为要与蜂后交配的雄蜂，将雄蜂结点按照与蜂后结点的边的聚集系数降序排列，排序后的雄蜂依次与蜂后交配，每交配一次，蜂后的能量和速度就衰减一次；计算某雄蜂与蜂后的交配成功概率P(i)，若P(i)>随机数rand，则交配成功，将该雄蜂的精子（雄蜂结点本身）加入蜂后的受精囊中（即归入类），并将该雄蜂结点标记visited＝1，即标记为已访问，交配不成功则继续与下一个雄蜂交配，直至蜂后的能量E小于能量阈值Thred或者蜂后的受精囊中的雄蜂结点个数大于M，一次婚飞过程结束。

(7)将蜂后受精囊中每个精子结点的邻接点中结点加权网络综合特征值com-value_i最大的结点保存下来，作为发育优良的幼蜂。

(8)从发育优良的幼蜂中选择结点加权网络综合特征值com-value_i最大的结点作为候选蜂后，候选蜂后作为新的蜂后，从而更新了聚类中心，以访问PPI网络中没有被访问的结点。

(9)如果iter<=maxiter，令iter=iter+1，并返回步骤(6)进行下一次婚飞；否则，转向步骤（10）；

(10)计算当前适应度fval，若当前适应度fval＜gfval，则令gcluster=cluster，gfval＝fvval；否则gcluster和gfval不变；

(11）如果count<=maxcount，令count=count+1，并返回步骤4），否则输出全局最优聚类结果gcluster。

表1、表2显示了本发明的CM_MBABC方法与功能流的Flow方法两种方法的聚类效果比较。两种算法均采自MIPS数据库的PPI数据集作为仿真数据集。Flow方法参照Proceedings of 6th IEEE International Conference on Data Mining（2006年《Identification of overlapping functional modules in protein interactionnetworks:information flow-based approach》）文献中记载的Flow方法。可以看出，与功能流的Flow方法相比，在效率（运行时间）、正确率、查全率以及几何平均值等指标上，本发明的CM_MBABC方法性能均较优。

表3给出了20个类聚类正确的蛋白质和聚类错误的蛋白质，可以看出，聚类序号为11、16两类聚类完全正确即这两类中的蛋白质都各自具有相同的功能；聚类序号为1、2、3、4、7、8、9、13、17和19这10类中聚类正确的蛋白质也特别多，聚类错误的蛋白质相对比较少；其他类只有部分蛋白质聚类正确。根据表3可以有效识别具有相同功能的蛋白质，方便人们有效识别蛋白质功能模块、预测蛋白质功能，为研究人类生命生理和病理条件下的变化机制提供可能。

表1 Flow方法的聚类结果

表2 CM-MBABC方法的聚类结果

表3本发明,的方法中某一类的聚类正确的和错误的蛋白质

Claims

1.一种基于人工蜂群繁殖机制的PPI网络聚类方法，其特征在于，具体包括如下步骤：

1)将PPI网络转化为一个无向加权图：

将PPI网络转化成一个无向加权图G(V，E)，其中，V＝{v_i,i＝1,2,…,n}为结点v_i的集合，E为边e的集合，结点v_i表示蛋白质，边e表示蛋白质之间的相互作用，w_ij表示结点v_i和结点v_j之间相互作用的大小，也就是结点v_i和结点v_j之间边e_ij的权值，若v_i和v_j之间没有边则w_ij＝0，i＝1,2,…,n，j＝1,2,…,n；

2)参数设置：

令count，maxcount分别表示控制外循环的当前迭代次数和外循环对应的最大迭代次数，maxcount∈[10,1000]，并令count＝1；iter，maxiter分别表示控制内循环的当前迭代次数和内循环对应的最大迭代次数，maxiter∈[10,200]；N和S分别表示蜂后婚飞时的能量和速度，N∈[50,1000]，S∈[10,500]；fval，gfval分别表示当前适应度和全局最优适应度，令gfval＝∞；cluster，gcluster分别表示当前的聚类结果和全局最优聚类结果；visited是结点被访问标记；n表示蛋白质结点个数，T表示蜂后与雄蜂交配成功的计数器；

3)对PPI网络的每个结点v_i和每个边e进行预处理：

所述对结点v_i预处理是计算结点v_i的加权网络综合特征值com-value_i；所述对边e预处理是计算改进的边聚集系数CC_i,j；

4)计算所有结点的加权网络综合特征值com-value_i的代数平均值Av-com-value，将大于Av-com-value*W的结点保存，其中W∈[0.5,2]；

5)初始化蜂后：

令iter＝1，蜂后代表聚类中心，初始化蜂后就是确定第一个聚类中心，从大于Av-com-value*W的结点中随机选取一个结点作为第一个蜂后，并令该蜂后结点的visited＝1；

6)婚飞过程：

P(i)＝exp[-Δ(f)/S(t)] 式

9

蜂后的速度S和能量E以式10和式11的方式衰减：

S(t+1)＝α*S(t) 式10

N(t+1)＝N(t)-β*N(t)/M 式11

7)幼蜂的局部搜索过程：

8)蜂后的选优过程：

9)iter＝iter+1，如果iter<＝maxiter，并返回步骤6)，否则，转向步骤10)；

10)计算当前适应度fval；若当前适应度fval<gfval，则令gcluster＝cluster，gfval＝fval；否则gcluster和gfval不变；

11)count＝count+1，如果count<＝maxcount，并令所有结点的访问标记visited＝0，并返回步骤4)，否则，输出全局最优聚类结果gcluster。

2.如权利要求1所述的基于人工蜂群繁殖机制的PPI网络聚类方法，其特征在于，步骤3)中所述的结点v_i的加权网络综合特征值com-value_i的计算过程如下：

结点v_i(i＝1,2,…,n)的加权度dw_i为与结点v_i连接的所有边的权值之和：

{dw}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{ij}

式1

结点v_i的聚集度dk_i定义为与结点v_i连接的近邻结点之间的连接的边的个数，定义如下：

{dk}_{i} = | {(v_{j}, v_{k}) | (v_{i}, v_{j}) &Element; E, (v_{i}, v_{k}) &Element; E, (v_{j}, v_{k}) &Element; E, v_{i}, v_{j}, v_{k} &Element; E} |

式2

结点v_i的加权聚集度为：

{wk}_{i} = \underset{(v_{i}, v_{j}) &Element; E}{Σ} w_{jk}

式3

其中，N_j＝{v_j|(v_i,v_j)∈E,(v_i,v_k)∈E,(v_j,v_k)∈E,v_i,v_j,v_k∈E}，j∈N_j，N_k＝{v_k|(v_i,v_j)∈E,(v_i,v_k)∈E,(v_j,v_k)∈E,v_i,v_j,v_k∈E}，k∈N_k；

{coe}_{i} = \frac{2 * {wk}_{i}}{d_{i} * (d_{i} - 1)}

式4

结点v_i的加权网络综合特征值com-value_i：

com-value_i＝α*coe_i+(1-α)*dw_i/n式

5

其中，i＝1,2,…,n,j＝1,2,......,n；(v_i,v_j)∈E表示结点v_i和v_j之间有相互作用；d_i表示与结点v_i连接的结点连接成三角形的个数；n为蛋白质结点的总个数；α是[0,1]之间的一个可调节参数，一般取0.5。

3.如权利要求1所述的基于人工蜂群繁殖机制的PPI网络聚类方法，其特征在于，步骤3)中所述的改进的边聚集系数CC_i,j的计算过程如下：

基于度的边的聚集系数C_i,j定义为：

C_{i, j} = \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)}

式6

加权聚集系数WC_u,v定义如下：

{WC}_{i, j} = \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式7

改进的边聚集系数CC_i,j定义为：

{CC}_{i, j} = ω * \frac{| N_{i} \cap N_{j} | + 1}{\min (| N_{i} |, | N_{j} |)} + (1 - ω) * \frac{\underset{k &Element; I_{i, j}}{Σ} w_{ik} \cdot \underset{k &Element; I_{i, j}}{Σ} w_{iv}}{\underset{s &Element; N_{i}}{Σ} w_{js} \cdot \underset{t &Element; N_{j}}{Σ} w_{jt}}

式8

4.如权利要求1所述的基于人工蜂群繁殖机制的PPI网络聚类方法，其特征在于，步骤10)中所述的当前适应度fval的计算过程如下：

PPI网络的子图密度度量如下：

D (t) = \frac{2 c_{t}}{n_{t} * (n_{t - 1})}

式12

两个功能模块M_s和M_t间的相似度矩阵S(M_s,M_t)通过式13计算：

S (M_{s}, M_{t}) = \frac{Σx &Element; M_{s}, y {&Element; M}_{t} c (x, y)}{\min (| M_{s} |, | M_{t} |)}

式13

其中，

c (x, y) = \{\begin{matrix} 1 & ifx = y \\ w_{e} (x, y) & ifx &NotEqual; yand < x, y > &Element; E \\ 0 & otherwise \end{matrix};

当前适应度fval定义如式14：

fval = \frac{1}{numclu} Σ_{t = 1}^{numclu} D (t) + \frac{Σ_{s = 1}^{numclu} Σ_{t = 1}^{numclu} S (M_{s}, M_{t})}{numclu * numclu}

式14