蛋白质复合体挖掘的加权组装聚类方法

摘要

本发明公开了一种蛋白质复合体挖掘的加权组装聚类方法，包括：输入一个蛋白质相互作用网络，产生一个无向图，选择m个聚类方法应用到这个网络上，得到m个聚类结果；对每个基本聚类结果，重新生成一个特征网络，得到m个特征网络，对应m个特征矩阵；对m个特征矩阵进行合成，获取合成矩阵W：其中uq是第q个特征网络的权重，uq≥0，且满足合成矩阵W对应一个新的网络，其中元素Wi,j是度量新网络中蛋白质i和蛋白质j的相似程度；采用贝叶斯非负矩阵分解算法挖掘该新网络中的聚类；把权重的学习和复合体发现整合为一个优化目标，从而通过聚类结果来优化权重，反之用权重结果来指导聚类；优化终止后获取最终的蛋白质复合体挖掘结果。

说明书

技术领域

本发明属于系统生物学领域，涉及蛋白质复合体的挖掘方法，更具体地， 涉及一种蛋白质复合体挖掘的加权组装聚类方法。

背景技术

蛋白质是分子功能的执行者和调控者，也是生命活动的主要载体。蛋白质 很少以单体的形式发挥作用，而是通过与其它蛋白质相互作用形成复合体或者 功能模块来协同地执行生物功能。蛋白质复合体的挖掘不仅有助于理解细胞的 功能组织机理，还有助于揭示复杂疾病的发生原理。尽管科研人员可以通过化 学实验测定方法较为准确的测定某一环境下的较为稳定的蛋白质复合体，但是 有些复合体内的蛋白质之间的相互作用是动态变化的，即存在不稳定蛋白质复 合体。以实验为基础的研究方法难以侦测到这类蛋白质复合体，并且需要耗费 大量的时间和昂贵的实验成本。近年来，高通量蛋白质相互作用预测技术的出 现产生了大量的蛋白质相互作用数据，从蛋白质相互作用网络出发侦测蛋白质 复合体已经成为了蛋白质复合体挖掘的一个主流方法。

目前，研究人员已经提出了一系列用于挖掘蛋白质复合体的图聚类算法。 这些算法主要基于分析蛋白质相互作用网络的拓扑结构来侦测蛋白质复合体。 依据其实现思想的不同，可以把这些算法粗略地分为三个子类：基于密度的局 部搜索算法、基于图划分的方法和基于层次聚类的方法。然而，由于每种算法 都有各自的特点，它们通常只能够捕捉到网络中某种特定的拓扑特征。例如， 基于密度的局部搜索算法通常只能侦测到蛋白质相互作用网络中紧密内连接的 子网络结构。然而研究表明，组织中的蛋白质复合体不仅限于紧密内连接的子 网络结构，蛋白质相互作用网络中也存在具有其它拓扑结构（例如星形和线形 结构）的蛋白质复合体。因此这些具有生物意义的低密度复合体往往会被传统 的基于密度的局部搜索算法所忽略。基于图划分的方法只能找到不重叠的蛋白 质复合体。而研究表明，蛋白质在不同的环境下可能执行不同的功能，因此蛋 白质复合体之间往往是相互重叠的。基于层次聚类的方法能够发现蛋白质相互 作用网络中的层次结构，但是这些方法对噪声数据相当敏感。而高通量技术得 到蛋白质相互作用数据不可避免地存在一定比例的噪声（假阳性和假阴性）。此 外，基于层次聚类的方法同样不能找到重叠的复合体。事实上，不同物种的不 同组织通过不同实验手段得到的蛋白质相互作用网络往往具有多种多样的网络 拓扑结构，因此很难找到一个算法在具有不同拓扑结构的蛋白质相互作用网络 数据上都有出色的表现。因此，如何设计一种方法来捕捉网络数据中的不同拓 扑特征，并且在不同的网络数据上都能得到较为准确的复合体侦测结果是一个 值得研究的问题。

目前有两种用于蛋白质复合体挖掘的组装聚类方法。一种是Asur等人提出 的方法，他们首先提出两种相似性度量来提高网络数据的精确度，然后使用三 种基于图划分的方法对改进的网络数据聚类，产生了六组聚类结果；最后，提 出一种基于主成分分析的一致聚类方法将上述结果融合为最终的聚类。该法需 要预先设定蛋白质复合体的个数，但是真实情况下，蛋白质相互作用网络中的 蛋白质复合体个数往往是未知的。此外，其使用的聚类方法都是基于图划分的 方法，这样，该法可能只能够捕捉到网络的某种拓扑特征，而忽略了其他重要 的特征。并且在组装不同聚类结果的时候，未对这些结果进行筛选，这样一来， 不可靠的聚类结果很可能会影响最终的聚类结果。

另一种是Greene等人提出的方法。通过设置不同的聚类个数，首先用非负 矩阵分解产生了一系列聚类结果；然后提出一个层级元聚类方法将这些聚类结 果融合为一系列不相交的“元聚类”，最后通过这些结果生成了原始网络的软层 级聚类。该方法主要是提高一个特定算法的稳定性，通过不同的初始值设置， 产生一系列聚类结果，进而综合成一个最终的聚类。由于没有对不同的聚类结 果进行筛选，而且仅使用了一种聚类方法，因此结果可能只反应出数据的部分 特征。

发明内容

本发明的主要目的是为了从不同的聚类方法得到的聚类结果中提取有效的 信息，并且产生一个更为准确和可靠的聚类结果，进而较为准确地侦测蛋白质 复合体。

为实现上述目的，本发明提出一种蛋白质复合体挖掘的加权组装聚类方法， 包括：

S1.输入一个蛋白质相互作用网络，产生一个无向图G，选择m个聚类方法 应用到这个网络上，得到m个聚类结果B_q，q=1,...,m；B=(B₁,B₂,...,B_m)；

S2.对每个基本聚类结果B_q，q=1,...,m；重新生成一个特征网络，得到m 个特征网络，m个特征网络对应m个特征矩阵；D=(D₁,D₂,...,D_m)；

S3.对上述m个特征矩阵进行加权组合，获取蛋白质相互作用网络的合成矩 阵W：其中u_q是第q个特征网络的权重，u_q≥0，q=1,...,m且 满足$\underset{q=1}{\overset{m}{\Sigma }}{u}_q=1;$

S4.合成矩阵W对应一个新的网络，其中元素W_i,j是度量新网络中蛋白质i 和蛋白质j的相似程度；采用贝叶斯非负矩阵分解算法挖掘该新网络中的聚类； 把权重的学习和复合体发现整合为一个优化目标，从而可以通过聚类结果来优 化权重，反之可以用权重结果来指导聚类；优化终止后获取最终的蛋白质复合 体挖掘结果。

其中步骤S1对输入的包含N个蛋白质的蛋白质相互作用网络建模，具体 采用无向图G：G=(N,E)来建模该网络，其中N个节点代表N个蛋白质，E 条边代表蛋白质直接交互的个数。本发明主要分为两个步骤:合成蛋白质相互作 用网络构建（提取不同聚类结果的有效信息）和蛋白质复合体挖掘。

步骤S3对m个特征矩阵进行合成，通过权重的设置有选择的组装不同聚类 方法得到的聚类结果，可以减弱不可靠聚类结果的干扰，增加可靠结果的对最 终聚类的影响。合成网络是不同基本聚类结果重构出的特征网络的加权组合， 因此可以把原始网络也作为一个特征网络加入组装聚类的过程。不仅充分利用 了原始数据的信息，还可以有效防止算法对基本聚类结果的过分依赖。合成网 络的数据特性恰好符合贝叶斯非负矩阵分解的模型假设，二者结合可以有效提 取不同基本聚类结果中的聚类信息。

本发明的模型求解过程中，通过贝叶斯推断估计模型中的参数，使用先验 分布既增强了模型的解释性，又减弱了模型对参数选择的敏感性和依赖性。最 后本发明还能发现重叠的蛋白质复合体，估算模型参数的同时自动估计出侦测 的聚类个数。

更进一步的，所述步骤S1还包括将蛋白质相互作用网络中未被第q个聚类 方法聚类的蛋白质设为单独的复合体，并添加到对应的聚类结果B_q中， q=1,...,m。此处采用将未被聚类的蛋白质设为单独的复合体，保证了每个聚 类结果都覆盖了所有的蛋白质。

更进一步的，所述特征矩阵D_q中的(D_q)_i,j代表第q个特征网络中第i个和 第j个节点间的状态，当第i个和第j个节点相连则(D_q)_i,j=1，否则， (D_q)_i,j=0，q=1,...,m。在各个特征网络中，两个节点相连当且仅当对应的 两个蛋白质至少同时出现在一个聚类中。其中的特征矩阵即为特征网络所对应 的邻接矩阵。

更进一步的，所述步骤S3是通过加权组合不同特征网络所对应的邻接矩阵 D_q（这里也称为特征矩阵），并且引入一个正则项来防止权重 过拟合某一个特征矩阵；初始化$U^{\left(0\right)}=\left(u_q^{\left(0\right)}\right),$令$u_q^{\left(0\right)}=\frac{1}{m},q=1,...,m.$

得到合成矩阵W后，该合成矩阵对于一个新的网络，其中W_i,j度量了该网 络中节点i和节点j的相似程度，而可能属于同一聚类的节点倾向具有较高的相 似度，即节点的类别信息影响了节点间的相似性，因此利用非负矩阵分解挖掘 这个网络中的聚类。利用非负矩阵分解进行聚类需要事先设定聚类个数，然而 网络中的聚类个数往往是未知的。因此在本发明中采用贝叶斯非负矩阵分解算 法挖掘该网络中的聚类，即聚类获取蛋白质复合体。

由于W_i,j的值表示了第i个蛋白质和第j个蛋白质在基本聚类结果中被聚 类到一起的频率，即它们属于同一复合体的可能性。则步骤S4的具体实现包括：

S41.令h_i,z表示第i个蛋白质属于第z个复合体的几率，设H=(h_i,z)表示蛋 白质-复合体倾向矩阵；另共有K个复合体，则表示第i个蛋白质和第 j个蛋白质属于同一复合体的几率；即可以用近似W_i,j;

S42.通过泊松噪声模型和独立性假设，得到：

$P\left(W\right|H)=\underset{i,j=1}{\overset{N}{\Pi }}\exp (-{\left({\mathrm{HH}}^T\right)}_{i,j})·{\left({\mathrm{HH}}^T\right)}_{i,j}^{W_{i,j}}/\Gamma (W_{i,j}+1)---\left(1\right)$

其中$H\in R_{N\times K}^{+};$

S43.假设h_i,z服从参数为β_z的半正态分布：

$P\left(h_{i,z}\right|{\beta }_z)=\sqrt{\frac{2}{\pi {\beta }_z}}\exp (-\frac{1}{2{\beta }_z}{h}_{i,z}^2),i=1,...,N,z=1,...,K---\left(2\right)$

采用β_z筛选聚类；假设β_z服从参数为a和b的inverse-Gamma分布：

$P\left({\beta }_z\right|a,b)=\frac{b^a}{\Gamma \left(a\right)}{\beta }_z^{-a-1}\exp (-\frac{b}{{\beta }_z}),z=1,...,K---\left(3\right)$

其中，a和b为用户设定的模型参数；

S44.综合上述模型，得到如下的联合概率分布P(W,H,β)：

P(W,H,β)=P(W|H)P(H|β)P(β)              （4）

其中$\beta =\left({\beta }_z\right)\in R_{K\times 1}^{+};$

S45.综合上述联合概率分布并加入正则项，得到如下目标函数：

$\underset{U,H,\beta }{\min }J(U,H,\beta )=-\log P(W,H,\beta )+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q$

$=-\log P\left(W\right|H)-\log P\left(H\right|\beta )-\log P\left(\beta \right)+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q---\left(5\right)$

s.t.H≥0,以及$\underset{q=1}{\overset{m}{\Sigma }}{u}_q=1,$u_q≥0，q=1,...,m。

其中λ为用户设定的控制正则项惩罚的平衡参数；

S46.通过独立性假设，将（1），（2），（3）代入（5）并去除常数部分，得 到具体形式的目标函数：

$\underset{U,H,\beta }{\min }J(U,H,\beta )=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left({\mathrm{HH}}^T\right)}_{i,j}-{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\frac{N}{2}\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z$

$+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{2{\beta }_z}{h}_{i,z}^2+(a+1)\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z+b\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{{\beta }_z}+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q---\left(6\right)$

s.t.H≥0,以及$\underset{q=1}{\overset{m}{\Sigma }}{u}_q=1,$u_q≥0，q=1,...,m；

S47.通过迭代更新来求解上述非负限制优化问题（6）；首先固定U的取值， 通过乘法更新准则（Multiplicative Updating Rule）对H和β进行更新；令φ_i,z为 限制h_i,z≥0所对应的拉格朗日乘子，记Φ=(φ_i,z)；拉格朗日函数L为：

$L(H,\beta ,\Phi )=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left({\mathrm{HH}}^T\right)}_{i,j}-{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\frac{N}{2}\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z$

$+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{2{\beta }_z}{h}_{i,z}^2+(a+1)\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z+b\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{{\beta }_z}+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}{\phi }_{i,z}{h}_{i,z}---\left(7\right)$

拉格朗日函数L关于h_i,z和β_z的梯度分别为：

${▿}_{h_{i,z}}L(H,\beta ,\Phi )=2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}-2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}+\frac{1}{{\beta }_z}{h}_{i,z}+{\phi }_{i,z}---\left(8\right)$

${▿}_{{\beta }_z}L(H,\beta ,\Phi )=-\frac{1}{2{\beta }_z^2}\underset{i=1}{\overset{N}{\Sigma }}{h}_{i,z}^2+\frac{N}{2{\beta }_z}-\frac{b}{{\beta }_z^2}+(a+1)\frac{1}{{\beta }_z}---\left(9\right)$

h_i,z和β_z的估计满足${▿}_{h_{i,z}}L(H,\beta ,\Phi )=0$和${▿}_{{\beta }_z}L(H,\beta ,\Phi )=0,$得到：

${\phi }_{i,z}=-2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}-\frac{1}{{\beta }_z}{h}_{i,z}---\left(10\right)$

${\beta }_z=\frac{\underset{i=1}{\overset{N}{\Sigma }}{h}_{i,z}^2+2b}{N+2a+2}---\left(11\right)$

通过Karush-Kuhn-Tucker(KKT)条件，φ_i,zh_i,z=0，得到如下关于h_i,z的方程：

$h_{i,z}[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]=h_{i,z}\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]---\left(12\right)$

则得到h_i,z的如下更新准则：

$h_{i,z}\leftarrow h_{i,z}\frac{\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]}{[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]}---\left(13\right)$

β_z的更新公式可以根据（11）得到；完成H和β的一次更新之后，固定H 和β的取值，对U进行更新；令γ为限制所对应的拉格朗日乘子；拉 格朗日函数L(U,γ)为：

$L(U,\gamma )=-\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q+\gamma (\underset{q=1}{\overset{m}{\Sigma }}{u}_q-1)$

令${▿}_{u_x}L(U,\gamma )=0$得到：

$u_x=\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_x\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)\exp (-1)\exp (-\frac{\gamma }{\lambda })$

由于得到u_x的更新公式：

$u_x=\frac{\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_x\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)}{\underset{q=1}{\overset{m}{\Sigma }}\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_q\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)}---\left(15\right)$

S48.根据更新公式（11），（13）和（15），能够通过迭代更新U，H和β的 数值来求解模型参数；首先初始化H=H⁽⁰⁾，其中每个元素随机抽取自（0,1） 上的均匀分布和U=U⁽⁰⁾；第t次迭代的时候，先固定U=U^(t-1)，通过（11） 和（13）更新β^(t)和H^(t)（先根据（11）使用H^(t-1)更新得到β^(t)，再根据（13） 使用U^(t-1)，H^(t-1)和β^(t)更新得到H^(t)）；得到β^(t)和H^(t)之后，固定它们的取 值，根据（15）得到U^(t)。如此不断迭代，直到满足设定的终止条件；

S49.设置最大迭代次数T和迭代终止条件||β^(t)-β^(t-1)||<ρ，t∈N⁺，其中 T和ρ为用户给定的参数，初始设置β⁽⁰⁾=0；当二者中的一方条件满足时，停 止迭代；得到H，β和U的估计值；

S410.根据β的数值大小，筛选出合适的聚类个数，即满足 的聚类，其中ρ_K是用户设定的参数；将满足 上述条件的聚类挑选出来，即从H中取出对应的列得到H'；由于H'的每个元 素都是实数值，通过阈值τ得到蛋白质-复合体指示矩阵其中：

此处，表示第i个蛋白质属于第z个侦测的复合体；反之，表 示第i个蛋白质不属于第z个侦测的复合体。

其中步骤S47中公式（13）采用下式替换

$h_{i,z}\leftarrow \frac{1}{2}{h}_{i,z}\frac{\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]}{[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]}+\frac{h_{i,z}}{2}---\left(14\right)$

则步骤S48的可采用的替换方式为：通过更新公式（11），（14）和（15）， 能够通过迭代更新U，H和β的数值来求解模型参数；首先初始化H=H⁽⁰⁾， 其中每个元素随机抽取自（0,1）上的均匀分布和U=U⁽⁰⁾；第t次迭代的时候， 先固定U=U^(t-1)，通过（11）和（14）更新β^(t)和H^(t)（先根据（11）使用H^(t-1)更新得到β^(t)，再根据（14）使用U^(t-1)，H^(t-1)和β^(t)更新得到H^(t)）；得到β^(t)和H^(t)之后，固定它们的取值，根据（15）得到U^(t)；如此不断迭代，直到满 足设定的终止条件。

其中步骤S43中β_z取值越接近0，对应的h_i,z，i=1,...,N的取值就越接近 0，即第z个聚类为空，因此可以利用β_z筛选聚类。为了更有效地估计β_z的取 值，考虑它的共轭先验信息，所以设β_z服从参数为a和b的inverse-Gamma分 布。

本发明的目的是为了从不同聚类方法得到的聚类结果中提取有效的信息， 并且产生一个更为准确和可靠的聚类结果。利用不同图聚类方法所捕捉到的网 络拓扑特征，并且对不同方法得到的聚类结果进行筛选，该方法能够实现更为 准确和可靠的蛋白质复合体挖掘。此外，该方法能够发现重叠的蛋白质复合体， 并且能够在优化过程中自动确定预测的蛋白质复合体的个数。该方法可用于提 高单个聚类方法的精度，由于有较高的灵活性，该方法可用于解决各类基于聚 类算法的应用问题。

与现有技术相比，本发明的有益效果为：

本发明能够根据不同的网络数据，对不同的聚类方法进行评估。进而选择 较为可靠的聚类结果参与组装聚类。本发明采用加权组装的方式组合不同的聚 类结果，并且在模型优化过程中自动调整权值，使得聚类结果更为准确和可靠。

本发明充分利用基本聚类结果中的有效信息。本发明使用的算法的模型假 设是如果两个节点之间有连接，则它们很有可能属于同样的聚类。而通过特征 网络加权组装得到的合成网络中两节点之间的连接恰好反映了它们在基本聚类 结果中被分配到同一聚类的频率。二者结合可以有效提取不同基本聚类结果中 的聚类信息。通过贝叶斯推断和先验分布的假设，本发明能够在估算模型参数 的同时自动估计出侦测的聚类个数，并且对模型参数选择的敏感度较低，有较 好的稳定性。

附图说明

图1为本发明的流程图。

图2为本发明部分参数的依赖关系图。

图3-5为本发明具体实施采用不同数据的结果显示图。

图6-7为本发明具体实施采用不同数据和不同组装聚类方法的结果示意 图。

图8-11为本发明具体实施采用Collins数据库中三个已知的重叠蛋白质复 合体被四种不同算法的检测图。

具体实施方式

下面结合附图对本发明做进一步的描述，但本发明的实施方式并不限于此。

本发明主要分为两个步骤:合成蛋白质相互作用网络构建（提取不同聚类结 果的有效信息）和蛋白质复合体挖掘。总体流程图如图1所示。具体步骤如下：

1.输入一个蛋白质相互作用网络，产生一个无向图G。选择m个聚类方法 应用到这个网络上，得到m个聚类结果（这里称为基本聚类结果）， B=(B₁,B₂,...,B_m)。由于有些聚类方法并没有覆盖到所有蛋白质，因此将每个 未被聚类的蛋白质设定为单独的复合体。这样一来，每个聚类结果都覆盖了所 有的蛋白质；

2.对于每一个基本聚类结果B_q，q=1,...,m，重新生成一个特征网络。在 这个网络中，两个节点相连当且仅当对应的两个蛋白质至少同时出现在一个聚 类中。这样，可以得到m个特征网络。通过这m个特征网络，可以对应得到m 个邻接矩阵（也称为特征矩阵）D=(D₁,D₂,...,D_m)。其中，(D_q)_i,j=1当且仅 当第q个特征网络中第i个和第j个节点相连，否则，(D_q)_i,j=0；

3.令合成蛋白质相互作用网络的邻接矩阵（也称为合成矩阵）其中u_q≥0，q=1,...,m是赋予每个特征网络的权重，且满足引入 一个正则项$R=\underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q.$

4.初始化$U^{\left(0\right)}=\left(u_q^{\left(0\right)}\right),$令$u_q^{\left(0\right)}=\frac{1}{m},$q=1,...,m。

5.得到合成矩阵W之后，其对应着一个新的网络。W_i,j度量了该网络中节 点i和节点j的相似程度，可能属于同一聚类的节点倾向具有较高的相似度，即 节点的类别信息影响了节点间的相似性，因此利用非负矩阵分解挖掘这个网络 中的聚类。利用非负矩阵分解进行聚类需要事先设定聚类个数，然而网络中的 聚类个数往往是未知的。因此采用贝叶斯非负矩阵分解算法挖掘该网络中的聚 类，即可能的蛋白质复合体。

6.注意到W_i,j的值表示了第i个蛋白质和第j个蛋白质在基本聚类结果中被 聚类到一起的频率，即它们属于同一复合体的可能性。令h_i,z表示第i个蛋白质 属于第z个复合体的可能性，h_i,z的取值越大，可能性越大。设H=(h_i,z)表示 蛋白质-复合体倾向矩阵。假设一共有K个复合体，则表示第i个蛋白 质和第j个蛋白质属于同一复合体的可能性。因此，可以近似认为

$W_{i,j}\approx {\hat{W}}_{i,j}={\left({\mathrm{HH}}^T\right)}_{i,j}.$

7.通过泊松噪声模型和独立性假设，得到：

$P\left(W\right|H)=\underset{i,j=1}{\overset{N}{\Pi }}\exp (-{\left({\mathrm{HH}}^T\right)}_{i,j})·{\left({\mathrm{HH}}^T\right)}_{i,j}^{W_{i,j}}/\Gamma (W_{i,j}+1)---\left(1\right)$

其中$H\in R_{N\times K}^{+}.$

8.假设h_i,z服从参数为β_z的半正态分布：

$P\left(h_{i,z}\right|{\beta }_z)=\sqrt{\frac{2}{\pi {\beta }_z}}\exp (-\frac{1}{2{\beta }_z}{h}_{i,z}^2),i=1,...,N,z=1,...,K---\left(2\right)$

则β_z取值越接近0，对应的h_i,z，i=1,...,N的取值就越接近0，即第z个 聚类为空。因此可以利用β_z筛选聚类。为了更有效地估计β_z的取值，考虑它的 共轭先验信息，即假设β_z服从参数为a和b的inverse-Gamma分布：

$P\left({\beta }_z\right|a,b)=\frac{b^a}{\Gamma \left(a\right)}{\beta }_z^{-a-1}\exp (-\frac{b}{{\beta }_z}),z=1,...,K---\left(3\right)$

其中，a和b为用户设定的模型参数；

9.综合上述模型，得到如下的联合概率分布P(W,H,β)：

P(W,H,β)=P(W|H)P(H|β)P(β)。       （4）

其中参数依赖关系如图2所示。

10.综合上述联合概率分布并加入正则项，得到如下目标函数：

$\underset{U,H,\beta }{\min }J(U,H,\beta )=-\log P(W,H,\beta )+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q$

$=-\log P\left(W\right|H)-\log P\left(H\right|\beta )-\log P\left(\beta \right)+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q---\left(5\right)$

s.t.H≥0,以及$\underset{q=1}{\overset{m}{\Sigma }}{u}_q=1,$uq≥0，q=1,...,m。

11.通过独立性假设，将(1)，(2)，(3)代入(5)并去除常数部分，得到具体形 式的目标函数：

$\underset{U,H,\beta }{\min }J(U,H,\beta )=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left({\mathrm{HH}}^T\right)}_{i,j}-{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\frac{N}{2}\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z$

$+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{2{\beta }_z}{h}_{i,z}^2+(a+1)\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z+b\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{{\beta }_z}+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q---\left(6\right)$

s.t.H≥0,以及$\underset{q=1}{\overset{m}{\Sigma }}{u}_q=1,$u_q≥0，q=1,...,m。

其中λ为用户设定的控制正则项惩罚的平衡参数。

12.通过迭代更新来求解上述非负限制优化问题（6）。首先固定U的取值， 通过乘法更新准则（Multiplicative Updating Rule）对H和β进行更新。令φ_i,z为 限制h_i,z≥0所对应的拉格朗日乘子且记Φ=(φ_i,z)。拉格朗日函数L为：

$L(H,\beta ,\Phi )=\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left({\mathrm{HH}}^T\right)}_{i,j}-{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\frac{N}{2}\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z$

$+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{2{\beta }_z}{h}_{i,z}^2+(a+1)\underset{z=1}{\overset{K}{\Sigma }}\log {\beta }_z+b\underset{z=1}{\overset{K}{\Sigma }}\frac{1}{{\beta }_z}+\underset{i=1}{\overset{N}{\Sigma }}\underset{z=1}{\overset{K}{\Sigma }}{\phi }_{i,z}{h}_{i,z}---\left(7\right)$

拉格朗日函数L关于h_i,z和β_z的梯度分别为：

${▿}_{h_{i,z}}L(H,\beta ,\Phi )=2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}-2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}+\frac{1}{{\beta }_z}{h}_{i,z}+{\phi }_{i,z}---\left(8\right)$

和${▿}_{{\beta }_z}L(H,\beta ,\Phi )=-\frac{1}{2{\beta }_z^2}\underset{i=1}{\overset{N}{\Sigma }}{h}_{i,z}^2+\frac{N}{2{\beta }_z}-\frac{b}{{\beta }_z^2}+(a+1)\frac{1}{{\beta }_z}---\left(9\right)$

因为h_i,z和β_z的估计应当满足${▿}_{h_{i,z}}L(H,\beta ,\Phi )=0$和${▿}_{{\beta }_z}L(H,\beta ,\Phi )=0,$得 到：${\phi }_{i,z}=-2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}-\frac{1}{{\beta }_z}{h}_{i,z}---\left(10\right)$

和${\beta }_z=\frac{\underset{i=1}{\overset{N}{\Sigma }}{h}_{i,z}^2+2b}{N+2a+2}---\left(11\right)$

通过Karush-Kuhn-Tucker(KKT)条件，φ_i,zh_i,z=0，得到如下关于h_i,z的方程：

$h_{i,z}[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]=h_{i,z}\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]---\left(12\right)$

因此很容易得到h_i,z的如下更新准则：

$h_{i,z}\leftarrow h_{i,z}\frac{\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]}{[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]}---\left(13\right)$

在实际操作中，根据Ding等的建议，以下的变换更新准则计算速度更快：

$h_{i,z}\leftarrow \frac{1}{2}{h}_{i,z}\frac{\left[2\underset{j=1}{\overset{N}{\Sigma }}{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}\frac{h_{j,z}}{{\left({\mathrm{HH}}^T\right)}_{i,j}}\right]}{[2\underset{j=1}{\overset{N}{\Sigma }}{h}_{j,z}+\frac{1}{{\beta }_z}{h}_{i,z}]}+\frac{h_{i,z}}{2}---\left(14\right)$

β_z的更新公式可以根据（11）得到。完成H和β的一次更新之后，固定它 们的取值，对U进行更新。令γ为限制所对应的拉格朗日乘子。拉格 朗日函数L(U,γ)为：

$L(U,\gamma )=-\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}[{\left(\underset{q=1}{\overset{m}{\Sigma }}{u}_q{D}_q\right)}_{i,j}·\log {\left({\mathrm{HH}}^T\right)}_{i,j}]+\lambda \underset{q=1}{\overset{m}{\Sigma }}{u}_q\log u_q+\gamma (\underset{q=1}{\overset{m}{\Sigma }}{u}_q-1)$

令${▿}_{u_x}L(U,\gamma )=0$得到：

$u_x=\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_x\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)\exp (-1)\exp (-\frac{\gamma }{\lambda })$

由于得到u_x的更新公式：

$u_x=\frac{\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_x\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)}{\underset{q=1}{\overset{m}{\Sigma }}\exp \left(\frac{1}{\lambda }\underset{i=1}{\overset{N}{\Sigma }}\underset{j=1}{\overset{N}{\Sigma }}{\left(D_q\right)}_{i,j}\log {\left({\mathrm{HH}}^T\right)}_{i,j}\right)}---\left(15\right)$

13.通过更新公式（11），（14）和（15），可以通过迭代更新U，H和β的 数值来求解模型参数。首先初始化H=H⁽⁰⁾（每个元素随机抽取自（0,1）上的 均匀分布）和U=U⁽⁰⁾。第t次迭代的时候，先固定U=U^(t-1)，通过（11）和 （14）更新β^(t)和H^(t)（先根据（11）使用H^(t-1)更新得到β^(t)，再根据（14） 使用U^(t-1)，H^(t-1)和β^(t)更新得到H^(t)）；得到β^(t)和H^(t)之后，固定它们的取 值，根据（15）得到U^(t)。如此不断迭代，直到满足设定的终止条件。

14.设置最大迭代次数T和迭代终止条件||β^(t)-β^(t-1)||<ρ，t∈N⁺（T和ρ 为用户给定的参数，初始设置β⁽⁰⁾=0）。当二者中的一方条件满足时，停止迭 代。得到H，β和U的估计值。

15.根据β的数值大小，筛选出合适的聚类个数，即满足

的聚类（这里ρ_K是用户设定的参数，取 值可以和ρ一样）。将满足上述条件的聚类挑选出来，即从H中取出对应的列 得到H'。由于H'的每个元素都是实数值，通过阈值τ得到蛋白质-复合体指示 矩阵其中：

这里，表示第i个蛋白质属于第z个侦测的复合体；反之，表 示第i个蛋白质不属于第z个侦测的复合体。

具体算法步骤如表1所示。

实施例一

选择九个经典的蛋白质复合体挖掘算法(ClusterONE,CMC,COPRA, DPClus,MCL,MCODE,MINE,RNSC,SPICi)作用在三个酵母蛋白质相互作用 网络数据库（Collins,Gavin和BioGRID）。使用两个参考蛋白质复合体数据库 （MIPS和SGD）和三个评估准则（f-measure,Jaccard和PR）来验证不同算法 结果的准确性。三个蛋白质相互作用网络数据库和两个参考数据库对应到这三 个网络的统计特征在表2和表3中显示。三个评估准则中，f-measure从蛋白 质复合体层面上度量了预测复合体和参考库中的复合体的相似程度。Jaccard 和PR从复合体-蛋白质层面上度量了预测复合体和参考库中的复合体的匹配情 况。

在描述这几个评估准则之前，我们先给出一些符号解释。另PP表示一个 算法所预测的复合体个数，PT表示参考数据库的复合体个数。C_i表示属于第i 个预测的复合体的蛋白质的集合，G_j表示属于第j个参考复合体的蛋白质的集 合。我们称预测复合体C_i和参考复合体G_j相匹配当且仅当：

$\frac{C_i\cap G_j}{C_i}>\delta$且$\frac{C_i\cap G_j}{G_j}>\delta .$

其中δ是（0，1）之间取值的阈值参数，这里固定为0.5。给定一组预测 蛋白质复合体V_P={C₁,C₂,...,C_PP}和一组参考复合体V_T={G₁,G₂,...,G_PT}， 查全率（Recall）和查准率（Precision）定义为：

为了综合考虑查全率（Recall）和查准率（Precision），定义f-measure 为查全率（Recall）和查准率（Precision）的调和平均，即：

$f-\mathrm{measure}=\frac{2\times \mathrm{Recall}\times \mathrm{Precision}}{\mathrm{Recall}+\mathrm{Precision}}.$

另外两个评估指标定义为：

Jaccard度量：令JaccardC_i=max_jJac_i,j， JaccardG_j=max_iJac_i,j。令${\mathrm{JaccardV}}_P=\frac{\underset{i=1}{\overset{\mathrm{PP}}{\Sigma }}\left|C_i\right|·{\mathrm{JaccardC}}_i}{\underset{i=1}{\overset{\mathrm{PP}}{\Sigma }}\left|C_i\right|},$

${\mathrm{JaccardV}}_T=\frac{\underset{j=1}{\overset{\mathrm{PT}}{\Sigma }}\left|G_j\right|·{\mathrm{JaccardG}}_j}{\underset{j=1}{\overset{\mathrm{PT}}{\Sigma }}\left|G_j\right|},$$\mathrm{Jaccard}=\frac{2\times {\mathrm{JaccardV}}_T\times {\mathrm{JaccardV}}_P}{{\mathrm{JaccardV}}_T+{\mathrm{JaccardV}}_P}.$

PR度量：令${\mathrm{PR}}_{i,j}=\frac{|C_i\cap G_j|}{\left|C_i\right|}\times \frac{{|C}_i\cap G_j|}{\left|G_j\right|},$PRC_i=max_jPR_i,j， PRG_j=max_iPR_i,j。${\mathrm{PRV}}_P=\frac{\underset{i=1}{\overset{\mathrm{PP}}{\Sigma }}\left|C_i\right|·{\mathrm{PRC}}_i}{\underset{i=1}{\overset{\mathrm{PP}}{\Sigma }}\left|C_i\right|},$${\mathrm{PRV}}_T=\frac{\underset{j=1}{\overset{\mathrm{PT}}{\Sigma }}\left|G_j\right|·{\mathrm{PRG}}_j}{\underset{j=1}{\overset{\mathrm{PT}}{\Sigma }}\left|G_j\right|},$$\mathrm{PR}=\frac{2\times {\mathrm{PRV}}_T\times {\mathrm{PRV}}_P}{{\mathrm{PRV}}_T+{\mathrm{PRV}}_P}.$

表2蛋白质相互作用网络的统计特征

表3参考蛋白质复合体数据库

确定参数，对于Collins数据库，K=500，τ=0.3，a=2，b=40对于 Gavin数据库，K=500，τ=0.3，a=2，b=20；对于BioGRID数据库， K=1000，τ=0.3，a=2，b=40。最大迭代次数T=150。ρ=ρ_K=1e-6。 正则化参数λ的选择利用了先验信息只需选择λ₀即 可，这样模型对λ₀的选取就相对不敏感。对于Collins，λ₀=0.5；对于Gavin， λ₀=1；对于BioGRID,λ₀=0.5。确定了参数之后，可以对比本发明与选取的 九个蛋白质复合体挖掘算法在三个数据库上的表现。结果显示在图3-5中，可 以看出，根据不同的评估指标和不同的参考数据库，本发明在不同类型的数据 库上都有较为稳定的表现。本发明预测的蛋白质复合体与参考数据库的匹配效 果也较好。

为了进一步验证本发明的有效性，图6-7显示了本发明与其他组装聚类方 法在不同数据库上的表现。这里选取的组装聚类方法是由Greene等人提出的算 法（ENMF）。除了分解聚类个数的选择区间和最终软层级的叶节点数这两个参 数之外，其他的参数都使用软件提供的默认参数。对于Collins，聚类个数的选 择区间设为[40，100]，对于Gavin，聚类个数的选择区间设为[80，150]。由于 在BioGRID上这个算法无法在48小时内输出结果，因此没有列举在BioGRID 上的结果。对于Collins和Gavin，叶节点的个数都设为80，100和120。从图 6-7可以看出本发明在不同数据库上都有较好的表现。

本发明还能够有效挖掘蛋白质相互作用网络中的重叠蛋白质复合体。图 8-11显示了Collins数据库中三个已知的重叠蛋白质复合体被四种不同算法的检 测情况。图中圆形代表属于RNA聚合酶Ⅰ复合体的蛋白质，矩形代表属于RNA 聚合酶Ⅱ复合体的蛋白质，三角形代表属于RNA聚合酶Ⅲ复合体的蛋白质，平 行四边形代表其他功能的蛋白质，六边形表示三个复合体共有的蛋白质，菱形 表示RNA聚合酶Ⅰ和RNA聚合酶Ⅲ共有的蛋白质。图8-11中的椭圆形区域表 示不同算法检测到的聚类。8图是由DPClus检测的，9图是由ClusterONE检测 的，10图是由MCODE检测的，11图示由本发明检测到的。可以发现，本发明 的方法能够更加准确的发现重叠蛋白质复合体。

以上所述的本发明的实施方式，并不构成对本发明保护范围的限定。任何 在本发明的精神原则之内所作出的修改、等同替换和改进等，均应包含在本发 明的权利要求保护范围之内。