CN105868582A - 采用果蝇优化方法识别蛋白质复合物 - Google Patents

Abstract

一种采用果蝇优化方法识别蛋白质复合物，由如下步骤组成：将蛋白质相互作用网络转化为无向图、对蛋白质相互作用网络边和结点预处理、构建动态蛋白质相互作用网络、设置参数、形成果蝇的位置、对应果蝇到蛋白质相互作用的网络中、确定初始化果蝇位置、确定果蝇气味浓度、更新果蝇位置、产生蛋白质复合物、过滤蛋白质复合物。该方法考虑了蛋白质网络的动态性、蛋白质复合物内部核心‐附件结构以及蛋白质相互作用网络的局部性和全局性，能准确地识别蛋白质复合物。仿真实验结果表明，正确率、查全率等指标性能较优。与其他聚类方法相比，结合蛋白质网络和蛋白质复合物特性，实现蛋白质复合物识别过程，提高了蛋白质复合物的识别准确率。

Description

采用果蝇优化方法识别蛋白质复合物

技术领域

本发明属于生物信息领域，具体涉及一种动态蛋白质相互作用网络中蛋白质复合物的识别方法。

背景技术

目前，随着高通量技术的诞生，大量的蛋白质相互作用(蛋白质相互作用)数据被检测出来，通过计算机检测蛋白质复合物以了解未知蛋白质的功能和预测疾病变得越来越重要。蛋白质之间的相互作用是随着细胞进入不同的生命周期进行着不同变化的，因此在聚类前构建一个能更真实地模拟蛋白质之间动态相互作用的网络对蛋白质复合物的识别起着至关重要的作用。目前通过计算机聚类识别蛋白质复合物的方法大多数作用于静态蛋白质相互作用网络，忽视了蛋白质之间联系的动态变化。

现有的蛋白质复合物识别方法有很多种，通过判别是否考虑到蛋白质相互作用网络的全局结构，聚类方法主要分为全局聚类方法和局部聚类方法。但由于蛋白质相互作用网络的小世界无尺度特性，方法本身对噪音数据敏感等缺陷，当前的这些蛋白质复合物识别方法准确率低，聚类效果不理想。

全局聚类方法是通过划整个蛋白质相互作用网络成为分开的独立子网来挖掘蛋白质复合物，其代表方法有G-N算法和马尔科夫(Markov Clustering algorithm,MCL)算法。G-N算法通过计算蛋白质相互作用网络中边的介数，然后不断移除网络中高介数边来分裂网络得到单个子网。该方法计算介数的代价高且未考虑复合物的局部特性，效率和蛋白质识别的准确率低。MCL算法主要是通过模拟随机游走来检测蛋白质复合物。随机游走从一个初始点开始，通过连接边的权重比例选择一个邻居节点移动。如果这个随机步进入一个高密度区域，将很难走出来。基于这个原理MCL通过(expansion和inflation)两个操作将蛋白质相互作用网络划分成不重叠的子网。MCL在蛋白质复合物的探测中表现出很好的性能，但其只能产生不重叠的子网，而蛋白质复合物之间的高度重叠，存在着很多功能复合物。全局聚类方法虽然考虑到了网络的全局性，但不能识别重叠的蛋白质复合物以及忽略掉蛋白质复合物的内部结构特性，是不可避免的缺点。

局部聚类方法检测蛋白质复合物主要是通过考虑局部邻居结点而不是全局网络进行聚类。主要有极大团(maximal clique algorithm,CMC)算法，小团体过滤法(Clique percolation Method,CPM),分子复合物检测(Molecular Complex Detection,MCODE)算法等。CMC聚类算法枚举所有最大簇的方法识别蛋白质复合物，枚举的时间复杂度高且忽略掉复合物中较小的类簇，性能低。CPM是基于团渗透的算法，将具有k-1个公共结点的全连通图进行合并以形成最终的簇。这种方法严格要求一个蛋白质复合物总是包含一个最大的团体，现实中的蛋白质复合却是大小不一的团体，会漏掉很多较小的蛋白质复合物。

一些寻找局部密度子图的方法被提出，这些方法先初始化簇为一个顶点，通过不同的启发标准将周围的邻居结点聚到簇中，并过滤掉密度较小的类簇。例如MCODE算法，先通过每个结点的邻居结点的局部密度给结点附权重值，将权重高的结点选出来当种子结点，通过不断扩充种子结点来形成最终的类簇。MCODE可以挖掘重叠的簇，时间复杂度低，对噪声数据的鲁棒性比较差。基于局部密度的蛋白质识别方法还有很多，如HC-PIN,SPIC,ClusterOne等等。

上述基于局部密度的局部聚类方法的蛋白质复合物的识别结果并不是很理想，因为其忽略了蛋白质复合物的内在结构。近年的研究表明，一个蛋白质复合物由一个核心簇和多个复件蛋白构成。基于蛋白质复合物的核心附件结构提出了核心聚类法(CORE)和基于核心附件的聚类法(Core-attachment based method,COACH)。CORE通过计算每一对蛋白质的P-Value来检测核心而COACH通过检测密度子图作为核心。对于每个核心的附件蛋白，如果一个蛋白质与核心里大于一半的蛋白质相连，认为这个蛋白质是这个核心的附件蛋白。COACH对蛋白质复合物的预测性能比其他忽视了蛋白质复合物内部结构的聚类方法更好。

最近提出一些新的基于蛋白质复合物核心附件结构的聚类方法，对于核心和附件的定义没有统一标准。局部聚类方法的优点是能检测重叠的蛋白质并且容易实现。但是他们产生蛋白质复合物主要通过一些启发式规则而没能考虑到整个蛋白质相互作用的网络特性。

上述聚类方法的缺陷没有考虑到蛋白质相互作用网络的动态性，不能同时考虑到整个蛋白质相互作用网络的全局特性和局部特性以及蛋白质复合物的内部结构，蛋白质复合物识别的准确率低。

发明内容

本发明要解决的技术问题在于克服上述蛋白质复合物识别方法的缺点，提供一种采用 果蝇觅食机理，准确度高、能真实地模拟蛋白质相互作用网络的动态性的采用果蝇优化方法识别蛋白质复合物。

解决上述技术问题所采用的技术方案由下述步骤组成：

1、将蛋白质相互作用网络转化为无向图

将蛋白质相互作用网络转化成一个无向图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相连，则w_ij＝1，若v_i与v_j之间没有边，则w_ij＝0，i＝1,2,…,n，j＝1,2,…,n。

2、对蛋白质相互作用网络边和结点的预处理

对结点v_i预处理：计算结点v_i的聚集系数和度，度是与结点v_i相邻的结点的个数，按式(1)计算结点v_i的聚集系数：

式中K_v是点v_i的度，n_v是点v_i的K_v个邻居结点间的边的个数；按式(2)计算边的聚集系数：

式中，Z(v_i,v_j)表示包含边(v_i,v_j)的三角形个数，d_i,d_j分别是点v_i,v_j的度；按式(3)计算边的皮尔森相关系数：

式中，x_i,y_i表示蛋白质v_x,v_y在时间点t时的基因表达值，μ(x),μ(y)是蛋白质v_x,v_y的平均基因表达值，T为时间点的最大值。

3、构建动态蛋白质相互作用网络

将蛋白质相互作用网络分为稳定关系部分和动态关系部分：

稳定关系部分：对于蛋白质相互作用网络中的每一条边e_ij计算其皮尔森相关系数，若边的皮尔森相关系数大于阈值th1＝0.3，则边e_ij被认为是稳定边，边所连接的两个结点被认为是稳定结点V_stab；用一个N×N的矩阵S来保存蛋白质相互作用网络的稳定关系，如果蛋白质v_i和蛋白质v_j有稳定的相互作用关系，S_ij＝1，否则S_ij＝0。

动态关系部分：在时间点t时，蛋白质v_i的基因表达值GE_it若大于基因表达阈值AT(i),则被认为蛋白质v_i在时间点t具有活性；若边e_ij∈E，GE_it≥AT(i),GE_jt≥AT(j)且S_ij＝0，则边e_ij为临时边，对应的蛋白为动态蛋白质。

通过上述处理，在每个时间点t的动态网络是一个包含稳定蛋白质和在时间点t处为活性的动态蛋白质的子网。

4、设置参数

令iter，maxiter分别表示当前迭代次数和最大迭代次数，iter＝1，matxiter∈[100,800]；t,T分别表示当前的时间点和最大时间点，t＝1；V_stab表示稳定结点，V_temp表示动态结点；N_core表示核心簇的计数器,N_core＝0；count,maxcount分别表示对稳定蛋白质遍历时当前的迭代次数和最大迭代次数，count＝1，maxcount＝V_stab的数量；核心簇的密度阈值th_core＝0.7；重叠度阈值为th_os；cluster是最终聚出来的蛋白质复合物。

5、形成果蝇的位置

对V_stab中的结点先按照点的聚集系数和度进行降序排序；当前的迭代次数count<最大迭代次数maxcount时，选取第count个结点Vi，若未被访问时，标记为被访问，结点Vi与其V_stab中的邻居结点形成一个核心簇core_i,按密度计算公式计算核心簇core_i的密度；若密度小于阈值th_core，不断将簇中与Vi相连拥有最小边的聚集系数的点去掉，直至core的密度大于等于th_core的值；将core中与Vi相连的点标记为visited＝1；count＝count+1,N_core＝N_core+1；循环对V_stab中的点进行访问，直到所有的点都被访问，形成了N_core个核心簇。

6、对应果蝇到蛋白质相互作用的网络中

将时间点t时对应网络动态结点中的每一个结点作为一只果蝇。

7、确定初始化果蝇位置

随机产生核心簇的序列号，即一只果蝇随机产生一个1…N_core的序列号，作为初始化果蝇位置X(i)，初始化果蝇的当前最好位置gbest_i为当前位置X(i)。

8、确定果蝇气味浓度

按照Smell(i)＝clossness(v_i,core(X(i)))确定果蝇气味浓度，式中Smell(i)为果蝇i的气味浓度，v_i是果蝇所代表的动态蛋白质结点，core(X(i))是果蝇当前所在位置所代表的核心簇，clossness表示v_i与core(X(i))的紧密程度。

9、更新果蝇位置

更新当前位置为X(i)+randomx，randomx作为一个随机产生的1～N_core之间的整数； 计算更新位置后的果蝇气味浓度Smell(i)，若当前位置的气味浓度大于果蝇在最佳位置的气味浓度，更新果蝇的当前最佳位置gbest_i为当前位置，否则不更新最佳位置；iter＝iter+1，若iter≦maxiter，重复步骤9,否则转向下一步。

10、产生蛋白质复合物

将每只果蝇最优位置的气味浓度作为最佳气味浓度，若最佳气味浓度大于0，则将这只果蝇代表的动态蛋白质并入到其最佳位置所代表的核心簇中；若最佳气味浓度等于0，则这只果蝇所代表的动态蛋白质被认为是噪音蛋白质；处理完所有果蝇，得到时间点t处的蛋白质复合物t＝t+1，若t≦T,返回步骤6，否则转向步骤11。

11、过滤蛋白质复合物

确定所有时间点产生的蛋白质复合物之间的重叠度OS，合并重叠度大于重叠度阈值th_os的蛋白质复合物，删除只有一个蛋白质的复合物，输出所有的蛋白质复合物cluster。

在本发明的步骤3中基因表达阈值AT(i)由式(4)得到：

AT(i)＝μ(i)+3σ(i)(1-F(i)) 式(4)

式中μ(i)是蛋白质v_i平均基因表达值，σ(i)是基因表达值的标准差，F(i)＝1/(1+σ²(i))是权函数。

在本发明的步骤5中核心簇core_i密度由式(5)得到：

式中e是子图中的边数，n是子图中的顶点数。

在本发明的步骤8中动态蛋白质结点v_i与核心簇core(X(i))之间紧密程度由式(6)得到：

式中，u是簇core(X(i))中的蛋白质。

在本发明的步骤11中蛋白质复合物之间的重叠度OS由式(7)得到：

式中，C_A,C_B分别表示蛋白质复合物A和蛋白质复合物B中顶点的集合。

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

1、本发明基于蛋白质复合物内部的核心附件结构，根据果蝇优化算法的高优化性能 将附件聚集到对应的核心中，构建蛋白质复合物识别模型，识别出来的蛋白质复合物与标准库相近，经实际数据验证，采用本发明方法识别出蛋白质复合物的准确度要比目前采用其它蛋白质复合物识别方法要高。

2、本发明不仅考虑到单个蛋白质的基因表达还考虑到蛋白质之间的相互联系，将蛋白质分成稳定蛋白质和动态临时存在的蛋白质，针对不同蛋白质的特性进行不同处理，更真实地模拟了蛋白质相互作用网络的动态性。

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

附图说明

图1是本发明实施例1的工艺流程模型图。

图2是标准蛋白质复合物的结构示意图。

图3是采用实施例1得出的蛋白质复合物的结构示意图。

具体实施方式

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

实施例1

以12个动态蛋白质网络为例采用果蝇优化方法识别蛋白质复合物步骤如下：

本实施例以采自DIP数据库的酵母数据集(DIP 20140427版)作为仿真数据集，DIP数据包含了4995个蛋白质和21554个相互作用关系。基因表达数据集采自GEO数据库中的酵母新陈代谢表达数据集GSE3431，其中包括6777个基因，3个周期共36个时间点的基因值，覆盖了DIP中的95％的蛋白质。采用基因表达值创建12个动态蛋白质相互作用网络。实验平台为Windows 7操作系统，Intel酷睿2双核3.1GHz处理器，4GB物理内存，用Matlab R2010b软件实现本发明的FOCA方法。

1、将蛋白质相互作用网络转化为无向图

将包含4995个蛋白质和21554个相互作用关系的蛋白质相互作用网络转化成一个无 向图G＝(V，E)，其中，V＝{v_i,i＝1,2,…,4995}为结点v_i的集合，E为21554个边e的集合，结点v_i表示蛋白质，边e表示蛋白质之间的相互作用，w_ij表示结点v_i与结点v_j之间的连接情况，若结点v_i与结点v_j相连则w_ij＝1，若v_i与v_j之间没有边，则w_ij＝0，i＝1,2,…,4995，j＝1,2,…,4995。

2、对蛋白质相互作用网络边和结点的预处理

对结点v_i预处理：i＝1,2,…,4995，每给定一个确定的i，可计算出4995个聚集系数和度，度是与结点v_i相邻的结点的个数，按式(1)计算结点v_i的聚集系数：

式中K_v是点v_i的度，n_v是点v_i的K_v个邻居结点间的边的个数；按(2)计算边的聚集系数：

式中，Z(v_i,v_j)表示包含边(v_i,v_j)的三角形个数，d_i,d_j分别是点v_i,v_j的度；按式(3)计算边的皮尔森相关系数：

式中，x_i,y_i表示蛋白质v_x,v_y在时间点t时的基因表达值，μ(x),μ(y)是蛋白质v_x,v_y的平均基因表达值，T为时间点的最大值12。

3、构建动态蛋白质相互作用网络

将蛋白质相互作用网络分为稳定关系部分和动态关系部分。

稳定关系部分：对于蛋白质相互作用网络中的每一条边e_ij计算其皮尔森相关系数，若边的皮尔森相关系数大于阈值th1＝0.3，则边e_ij被认为是稳定边，边所连接的两个结点被认为是稳定结点V_stab；用一个1873×1873的矩阵S来保存蛋白质相互作用网络的稳定关系，如果蛋白质v_i和蛋白质v_j有稳定的相互作用关系，S_ij＝1，否则S_ij＝0；

动态关系部分：在时间点t时，蛋白质v_i的基因表达值GE_it若大于基因表达阈值AT(i),则被认为蛋白质v_i在时间点t具有活性；若边e_ij∈E，GE_it≥AT(i),GE_jt≥AT(j)且S_ij＝0，则边e_ij为临时边，对应的蛋白为动态蛋白质，基因表达阈值AT(i)由式(4)得到：

AT(i)＝μ(i)+3σ(i)(1-F(i)) 式(4)式中μ(i)是蛋白质v_i基因表达值，σ(i)是基因表达值的标准差，F(i)＝1/(1+σ²(i))是权函数通过 上述处理，在时间点t的动态网络G^(t)＝(V，E^(t)),t＝1,…,12,其中V表示全体蛋白质结点，E^(t)在时间点t时出现的边；边e_ij ^(t)∈E^(t)如果S_ij＝1(表示是稳定关系)或者e_ij∈E，GE_it≥AT(i),GE_jt≥AT(j)(表示临时关系)。整个网络被分成12个动态蛋白质的子网。

4、设置参数

令iter，maxiter分别表示当前迭代次数和最大迭代次数，iter＝1，matxiter∈[100,800]；t,T分别表示当前的时间点和最大时间点，t＝1，T＝12；V_stab表示稳定结点，V_temp表示动态结点；N_core表示核心簇的计数器,N_core＝0；count,maxcount分别表示对稳定蛋白质遍历时当前的迭代次数和最大迭代次数，count＝1，maxcount＝V_stab的数量1873；核心簇的密度阈值th_core＝0.7；cluster是最终聚出来的蛋白质复合物；重叠度阈值th_os＝0.8。

5、形成果蝇的位置

对V_stab中1873个结点先按照点的聚集系数和度进行降序排序；当前的迭代次数count<最大迭代次数maxcount1873时，选取第count个结点Vi，若未被访问时，标记为被访问，结点Vi与其V_stab中的邻居结点形成一个核心簇core_i,确定核心簇core的密度，核心簇core_i密度由式(5)得到：

式中e是子图中的边数，n是子图中的顶点数。若密度小于阈值0.7，不断将簇中与Vi相连拥有最小边的聚集系数的点去掉，直至core的密度大于等于0.7的值；将core中与Vi相连的点标记为visited＝1；count＝count+1,N_core＝N_core+1；循环对V_stab中的点进行访问，直到所有的点都被访问，形成了N_core个核心簇。

6、对应果蝇到蛋白质相互作用的网络中

将时间点t，t为1…12时，对应网络动态结点中的每一个结点作为一只果蝇，在不同的时间点t对应的果蝇的数量见表1。

表1 12个动态蛋白质子网分别对应的果蝇数量

7、确定初始化果蝇位置

随机产生核心簇的序列号，即一只果蝇随机产生一个1…1183的序列号，1183为核心簇的个数，作为初始化果蝇位置X(i)，初始化果蝇的当前最好位置gbest_i为当前位置。

8、确定果蝇气味浓度

按照Smell(i)＝clossness(v_i,core(X(i)))

确定果蝇气味浓度，式中Smell(i)为果蝇i的气味浓度，v_i是果蝇所代表的动态蛋白质结点，core(X(i))是果蝇当前所在位置所代表的核心簇，clossness表示v_i与core(X(i))的紧密程度，动态蛋白质结点v_i与核心簇core(X(i))之间紧密度由式6得到：

式中，u是簇core(X(i))中的蛋白质。

9、更新果蝇位置

更新当前位置为X(i)+randomx，randomx为一个随机产生的1～1183之间的整数；计算更新位置后的果蝇气味浓度Smell(i)，若当前位置的气味浓度大于果蝇在最佳位置的气味浓度，更新果蝇的当前最佳位置gbest_i为当前位置，否则不更新最佳位置；iter＝iter+1，若iter≦800，重复步骤9,否则转向下一步；

10、产生蛋白质复合物

将每只果蝇最优位置的气味浓度作为最佳气味浓度(最佳气味浓度为大于0的值)，若最佳气味浓度大于0，则将这只果蝇代表的动态蛋白质并入到其最佳位置所代表的核心簇中；若最佳气味浓度等于0，则这只果蝇所代表的动态蛋白质被认为是噪音蛋白质；处理完所有果蝇，得到时间点t处的蛋白质复合物t＝t+1，若t≦T,返回步骤6，否则转向下一步骤。

11、过滤蛋白质复合物

确定时间点产生的蛋白质复合物之间的重叠度OS，蛋白质复合物之间的重叠度OS由式7得到：

式中，C_A,C_B分别表示蛋白质复合物A和蛋白质复合物B中顶点的集合。合并重叠度大于th_os为0.8的蛋白质复合物。删除只有一个蛋白质的复合物，输出所有的蛋白质复合物cluster。

为了验证本发明的有益效果，发明人采用本发明实施例1果蝇优化机理的蛋白质复合物识别方法对DIP数据库中的蛋白质网络进行蛋白质复合物的识别，共识别出707个蛋白质复合物，结果见表2、表3、图2、图3，表2中选取了8个蛋白质复合物中正确和错误的蛋白质以及对应的标准库中的蛋白质进行分析。

表2本发明识别的蛋白质复合物中正确的和错误的蛋白质

表2显示了本发明识别出的结果中的8个蛋白质复合物中正确和错误的蛋白质以及对应的标准库中的蛋白质复合物，蛋白质复合物3、7完全正确地识别出标准库中的蛋白质复合物；蛋白质复合物1、2、4、5、8识别的蛋白质复合物中都只有一个错误蛋白质，正确的蛋白质特别多；蛋白质复合物1、6少识别了一个蛋白质，蛋白质复合物8少识别出两个蛋白质。由表2看出，本发明能有效地识别蛋白质复合物。图2显示了一个标准发库中的蛋白质复合物和本发明识别出来的一个蛋白质复合物。图2中两带背景颜色的是本明未识别出来的蛋白质，图3中带背颜色的是识别错误的蛋白质。中间的椭圆中蛋白质为核心簇的蛋白。由此图可见，本发明能正确地识别蛋白质复合物中大部分蛋白质。

表3其他方法的对蛋白质复合物识别结果对比表

表3显示了本发明跟当前其他蛋白质复合物识别方法在性能上的对比。结果显示本发明无论是准确率还是查全率都比其它蛋白质复合物的性能要好。

Claims (5)

1.一种采用果蝇优化方法识别蛋白质复合物，其特征在于，它是由下述步骤组成：

(1)将蛋白质相互作用网络转化为无向图

将蛋白质相互作用网络转化成一个无向图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相连，则w_ij＝1，若v_i与v_j之间没有边，则w_ij＝0，i＝1,2,…,n，j＝1,2,…,n；

(2)对蛋白质相互作用网络边和结点的预处理

对结点v_i预处理：计算结点v_i的聚集系数和度，度是与结点v_i相邻的结点的个数，按式(1)计算结点v_i的聚集系数：

式中K_v是点v_i的度，n_v是点v_i的K_v个邻居结点间的边的个数；按式(2)计算边的聚集系数：

式中，Z(v_i,v_j)表示包含边(v_i,v_j)的三角形个数，d_i,d_j分别是点v_i,v_j的度；按式(3)计算边的皮尔森相关系数：

式中，x_i,y_i表示蛋白质v_x,v_y在时间点t时的基因表达值，μ(x),μ(y)是蛋白质v_x,v_y的平均基因表达值，T为时间点的最大值；

(3)构建动态蛋白质相互作用网络

将蛋白质相互作用网络分为稳定关系部分和动态关系部分：

稳定关系部分：对于蛋白质相互作用网络中的每一条边e_ij计算其皮尔森相关系数，若边的皮尔森相关系数大于阈值th1＝0.3，则边e_ij被认为是稳定边，边所连接的两个结点被认为是稳定结点V_stab；用一个N×N的矩阵S来保存蛋白质相互作用网络的稳定关系，如果蛋白质v_i和蛋白质v_j有稳定的相互作用关系，S_ij＝1，否则S_ij＝0；

动态关系部分：在时间点t时，蛋白质v_i的基因表达值GE_it若大于基因表达阈值AT(i),则被认为蛋白质v_i在时间点t具有活性；若边e_ij∈E，GE_it≥AT(i),GE_jt≥AT(j)且S_ij＝0，则边e_ij为临时边，对应的蛋白为动态蛋白质；

通过上述处理，在每个时间点t的动态网络是一个包含稳定蛋白质和在时间点t处为活性的动态蛋白质的子网；

(4)设置参数

令iter，maxiter分别表示当前迭代次数和最大迭代次数，iter＝1，matxiter∈[100,800]；t,T分别表示当前的时间点和最大时间点，t＝1；V_stab表示稳定结点，V_temp表示动态结点；N_core表示核心簇的计数器,N_core＝0；count,maxcount分别表示对稳定蛋白质遍历时当前的迭代次数和最大迭代次数，count＝1，maxcount＝V_stab的数量；核心簇的密度阈值th_core＝0.7；重叠度阈值为th_os；cluster是最终聚出来的蛋白质复合物；

(5)形成果蝇的位置

对V_stab中的结点先按照点的聚集系数和度进行降序排序；当前的迭代次数count<最大迭代次数maxcount时，选取第count个结点Vi，若未被访问时，标记为被访问，结点Vi与其V_stab中的邻居结点形成一个核心簇core_i,按密度计算公式计算核心簇core_i的密度；若密度小于阈值th_core，不断将簇中与Vi相连拥有最小边的聚集系数的点去掉，直至core的密度大于等于th_core的值；将core中与Vi相连的点标记为visited＝1；count＝count+1,N_core＝N_core+1；循环对V_stab中的点进行访问，直到所有的点都被访问，形成了N_core个核心簇；

(6)对应果蝇到蛋白质相互作用的网络中

将时间点t时对应网络动态结点中的每一个结点作为一只果蝇；

(7)确定初始化果蝇位置

随机产生核心簇的序列号，即一只果蝇随机产生一个1…N_core的序列号，作为初始化果蝇位置X(i)，初始化果蝇的当前最好位置gbest_i为当前位置X(i)；

(8)确定果蝇气味浓度

按照Smell(i)＝clossness(v_i,core(X(i)))确定果蝇气味浓度，式中Smell(i)为果蝇i的气味浓度，v_i是果蝇所代表的动态蛋白质结点，core(X(i))是果蝇当前所在位置所代表的核心簇，clossness表示v_i与core(X(i))的紧密程度；

(9)更新果蝇位置

更新当前位置为X(i)+randomx，randomx作为一个随机产生的1～N_core之间的整数；计算更新位置后的果蝇气味浓度Smell(i)，若当前位置的气味浓度大于果蝇在最佳位置的气味浓度，更新果蝇的当前最佳位置gbest_i为当前位置，否则不更新最佳位置；iter＝iter+1，若iter≦maxiter，重复步骤(9),否则转向下一步；

(10)产生蛋白质复合物

将每只果蝇最优位置的气味浓度作为最佳气味浓度，若最佳气味浓度大于0，则将这只果蝇代表的动态蛋白质并入到其最佳位置所代表的核心簇中；若最佳气味浓度等于0，则这只果蝇所代表的动态蛋白质被认为是噪音蛋白质；处理完所有果蝇，得到时间点t处的蛋白质复合物t＝t+1，若t≦T,返回步骤(6)，否则转向步骤(11)；

(11)过滤蛋白质复合物

确定所有时间点产生的蛋白质复合物之间的重叠度OS，合并重叠度大于重叠度阈值th_os的蛋白质复合物，删除只有一个蛋白质的复合物，输出所有的蛋白质复合物cluster。

2.根据权利要求1所述的采用果蝇优化方法识别蛋白质复合物，其特征在于所述的步骤(3)中基因表达阈值AT(i)由式(4)得到：

AT(i)＝μ(i)+3σ(i)(1-F(i)) 式(4)

式中μ(i)是蛋白质v_i平均基因表达值，σ(i)是基因表达值的标准差，F(i)＝1/(1+σ²(i))是权函数。

3.根据权利要求1所述的采用果蝇优化方法识别蛋白质复合物，其特征在于所述的步骤(5)中核心簇core_i密度由式(5)得到：

式中e是子图中的边数，n是子图中的顶点数。

4.根据权利要求1所述的采用果蝇优化方法识别蛋白质复合物，其特征在于所述的步骤(8)中动态蛋白质结点v_i与核心簇core(X(i))之间紧密程度由式(6)得到：

式中，u是簇core(X(i))中的蛋白质。

5.根据权利要求1所述的采用果蝇优化方法识别蛋白质复合物，其特征在于所述的步骤(11)中蛋白质复合物之间的重叠度OS由式(7)得到：

式中，C_A,C_B分别表示蛋白质复合物A和蛋白质复合物B中顶点的集合。