采用隐马尔科夫模型的蛋白质复合物识别研究

李 鹏，罗爱静，闵 慧，谭荪怡，郭惠敏

1.中南大学 湘雅三医院，长沙 410006

2.湖南中医药大学 信息科学与工程学院，长沙 410208

3.医学信息研究湖南省普通高等学校重点实验室（中南大学），长沙 410006

4.湖南信息职业技术学院 软件学院，长沙 410200

随着人类基因组计划以及对于多个物种全基因组测序工作的结束，目前生命科学研究的重点已经转移到了蛋白组学[1]。蛋白质是生物体生命活动的重要物质基础，没有蛋白质就没有生命。一个生物体内所有蛋白质相互作用被称为蛋白质相互作用网络（protein-protein interaction networks，PPIN），简称蛋白质网络[2]。值得注意的是，蛋白质之间的相互作用是动态的，它会随着时间环境、蛋白质的存在和降解、细胞的不同生理状态等因素的变化而变化，但由于PPIN 本身的复杂性、可利用的蛋白质相互作用数据的不完全性和存在噪声等众多问题，准确且高效地衡量蛋白质相互作用的动态性还存在很多挑战[3-4]，这也直接限制了PPIN 领域内很多问题[5-6]（复合物识别、关键蛋白识别、网络比对等）的研究进展。为此，本文将对动态蛋白质网络的构建问题进行研究，并在此基础上提出新的蛋白质复合物识别方法，这有利于从微观层面解释细胞内蛋白质之间的复杂关系，为生物学家、医学家理解生命复杂网络的内在组织和生物过程提供新的途径，在药物标靶设计、疾病诊治和预测等方面也具有重要的应用价值。

1 相关工作

目前，已有众多研究者对动态蛋白质网络中的复合物识别问题进行了研究，例如雷秀娟等人[7]首先将所有蛋白质之间的相互作用建模为一个拓扑势场，并以此来构建加权的动态蛋白质网络。在此基础上再采用马尔科夫聚类算法在多种典型的生物数据集上进行复合物识别。然而该方法存在复杂度较大的问题。赵学武等人[8]分析了在静态蛋白质网络上进行复合物识别的局限性，提出了一种基于时序基因表达数据和蚁群聚类的复合物识别算法（identify protein complexes by integrating temporal function continue feature and ant colony clustering，IPC-TFC&ACC）应用于动态蛋白质网络中。该算法首先基于蛋白质的表达活性选出复合物的种子节点，然后结合蛋白质之间的功能相似性采用蚁群聚类算法来得到最终的复合物。文献[9]从蛋白质网络的拓扑属性和生物学特征着手，提出了基于动态核依附的蛋白质复合物识别算法（dynamic core-attachment，DCA）。文献[10]首先利用蛋白质的功能注释数据和时序基因表达数据来构建动态蛋白质网络，然后基于聚类的方法得到初始的复合物，在此基础上设计出蛋白质节点的迁移和复合物间的融合策略，进而提出了一种基于动态群体决策的复合物识别算法（identifying protein complexes based on dynamic group decision，IPC-DGD）。Lei 等人[11]提出了一种改进的花授粉算法（improved flower pollination algorithm，IFPA）以识别多关系重建的动态PPI 网络中的蛋白质复合物。文中首先考虑了蛋白质在寻找本质相互作用中的重要性，然后设计了多关系重构的动态PPI 网络，并发现了网络中蛋白质复合物的潜在核心。最后，提出了一种基于花授粉机制的IFPA 算法，通过模拟花粉过程生成蛋白质复合物。该方法可在扩充过程中允许某个蛋白质重复出现，但无法识别PPI 网络中非稠密的子图结构，识别精度还有待提高。

虽然上述方法各有优点，在某些生物数据集上表现出较好的计算性能和准确性，但这些算法只是单独对每个时刻的蛋白质网络快照进行蛋白质复合物识别，忽略了动态蛋白质网络中不同时刻之间的影响关系。在动态蛋白质网络中，当前时刻的网络结构信息与前一时刻或者前几个时刻的网络结构信息是相关的，前一时刻的蛋白质网络中的复合物对当前时刻蛋白质网络中的复合物存在着影响，而之前的动态算法忽略了历史信息对当前时刻网络结构的影响，使得复合物识别的准确性和稳定性偏低。针对上述识别算法存在的问题，文中借鉴隐马尔科夫模型（hidden Markov model，HMM）[12]在处理复杂网络上的优势，提出一种基于HMM 的蛋白质复合物识别算法（HMM-PC）。该方法首先引入HMM 描述蛋白质复合物结构与网络中蛋白质个体间的相互作用，将网络中的复合物结构和蛋白质节点信息分别采用状态链和观察链表示。然后通过求解HMM 中的最优状态序列实现动态蛋白质网络中的复合物识别。最后通过在多个生物数据集上的实验验证了HMM-PC 算法的有效性。

2 动态加权蛋白质网络的构建

蛋白质网络具有偏好连接特性，且随着网络规模的增加，不同节点间的连接数差异呈现着增加的趋势。因此本文通过构建加权的动态蛋白质网络来反映真实网络结构，在此基础之上识别蛋白质复合物。

2.1 初始蛋白质网络的构建

一般而言，蛋白质网络可用无向图进行建模，设G=(V,E) 表示蛋白质网络，V表示G中的蛋白质集合，|V|表示G中的蛋白质总数；E表示G中的蛋白质与蛋白质之间相互作用的集合，|E|表示G中的相互作用的总数。考虑到蛋白质只有在活性状态下才能与其他蛋白质发现相互作用，文中引入文献[13]提出的3-sigma法则来判断在一个时间序列中哪些蛋白质是活性的，然后利用蛋白质之间的基因共表达特性来构建初始蛋白质网络。文献[14]对酵母菌株的生长过程进行了耗氧量测定，验证了周期基因表达的存在，可以将蛋白质36 个时刻的基因表达分为3个周期，每12 个时刻为一个周期。设每个蛋白质的基因表达平均水平为：

其中，tv(i)表示蛋白质v在时刻i的基因表达值。计算出所有蛋白质的基因表达平均水平后，根据3-sigma 法则来判断各个蛋白质在12 个时刻的活性状态。如果蛋白质u和蛋白质v在第i(1 ≤i≤12)个时刻同时处于活性状态，则认为u、v之间在时刻i存在相互作用。即有：

为了构建第i个时刻的蛋白质网络，采用如下的思路进行构建：统计第i个时刻同时处于活性状态的所有蛋白质的数目，并采用邻接矩阵进行存储，然后根据式（2）确定邻接矩阵中对应位置的元素的值，从而可以得到一个无向图来表示蛋白质网络。

2.2 网络加权

通过基因表达数据构建得到的初始蛋白质网络只能衡量节点之间的连接关系，但不能反映出节点间的偏好连接特性。动态蛋白质网络可由多个时刻的静态蛋白质网络表示，本文将动态蛋白质网络按照时间粒度划分为多个网络快照。有如下的定义：

定义1（时序图时序图表示在时刻t内具有相互作用的所有蛋白质构成的蛋白质网络。动态蛋白质网络可以用时序图的集合来描述：

文中在2.1 节构建得到的12 个初始蛋白质网络的基础上，考虑进一步利用蛋白质的生物信息数据和拓扑结构数据来对网络进行加权。文中首先对酵母物种中部分已知的蛋白质复合物与蛋白质功能注释的关系进行了统计分析，结果如图1 所示。可以看到，随着不同蛋白质之间共享功能注释数目的增加，能够识别的蛋白质复合物数量的比例和蛋白质复合物规模的比例都在呈现上升的趋势，这表明蛋白质间的共享功能注释信息有助于识别蛋白质复合物，可以应用到复合物的识别中。为此，本文给出了如下的定义2 作为对蛋白质网络加权的因子之一。

定义2（共享功能注释）对于PPIN 中的任意两个蛋白质u和v表示u的功能注释集合；表示v的功能注释集合；表示u和v的相同功能注释集合。则有：

Fig.1 Functional annotation of protein vs protein complex图1 蛋白质复合物vs蛋白质的功能注释

考虑到大部分物种中蛋白质的已知功能注释数量不一，为了保证复合物识别的准确性，本文认为，如果某个蛋白质的功能注释数目小于2，则认为利用该功能注释来对网络加权无意义，不做考虑。

定义3（共享结构域[15]）对于PPIN 中的任意两个蛋白质u和v，表示u的结构域集合；表示v的结构域集合；表示u和v的相同结构域集合。则有：

定义4（连接强度）对于PPIN 中的任意两个蛋白质u和v，dcuv是u和v之间直接相连的边的数量；nb(u)是蛋白质u拥有的邻居节点数；du是蛋白质u的度。则有：

联立定义2～4 可知，动态蛋白质网络中任意两个蛋白质u和v之间边的权值为：

综上所述，本文提出的动态加权蛋白质网络构建过程如算法1 所示。

算法1动态加权蛋白质网络构建算法

输入：蛋白质相互作用数据，基因表达数据。

输出：动态蛋白质网络G。

1.首先根据式（1）计算所有蛋白质的基因表达水平，然后采用3-sigma 法则判断各个蛋白质在第i(1 ≤i≤12)个时刻的活性状态；

2.i=1，根据式（2）找出在同一时刻存在相互作用的所有蛋白质对，构建得到第1 个时刻的初始蛋白质网络G′；

3.对G′中的所有边进行如下处理：

3.1 根据式（4）计算G′中所有蛋白质的共享功能注释信息；

3.2 根据式（5）计算G′中所有蛋白质的共享结构域信息；

3.3 根据式（6）计算G′中所有蛋白质的连接强度信息；

3.4 根据式（7）对G′中所有边进行赋权。

4.i=i+1；

5.重复执行2～3，直到所有12 个时刻的网络都赋权完毕，则算法结束。

3 基于HMM 的蛋白质复合物识别

3.1 HMM 介绍

隐马尔科夫模型（HMM）是一种用来描述随机过程统计特性的概率模型，它是一个关于时序的概率模型，描述由一个隐藏的马尔科夫链随机生成不可观测的状态随机序列，再由各个状态生成一个观测而产生观测随机序列的过程。HMM 由初始状态概率向量、状态转移概率矩阵和观测概率矩阵组成，如图2 所示。其中，x表示隐藏状态，y是观测状态。

Fig.2 Example of HMM图2 HMM 示例

假定相邻时刻的蛋白质网络结构相关，对不同时刻蛋白质网络中复合物结构的关系可用一阶Markov 链描述，由于蛋白质网络中的复合物结构不能被直接观察到，文中在考虑前一时刻蛋白质网络拓扑结构信息对当前时刻蛋白质网络拓扑结构信息影响的前提下，引入HMM 将不同时刻蛋白质网络中的复合物结构建模为HMM 中的状态链。将蛋白质网络中每个蛋白质节点的节点度分布建模为观察链。采用HMM 描述动态蛋白质网络中的复合物发现问题，通过求解HMM 中的最优状态序列得到蛋白质复合物，可以更好地理解动态蛋白质网络中的复合物结构演变情况。为了便于描述，先给出文中用到的如下几个定义：

定义5（状态序列SS）基于定义1 所示的动态蛋白质网络G，定义G中各个时刻的蛋白质复合物为HMM 中的状态序列SS。其中，SS=(C1,C2,…,Ct,…,CN)，N表示状态总数，Ct表示蛋白质节点在时刻t的状态。

定义6（可观察序列OS）基于定义1 所示的动态蛋白质网络G，定义HMM 中的观察序列为G中每个蛋白质节点的节点度大小。其中，OS=(d1,d2,…,dt,…,dM)，M表示从每一状态可能输出的不同的观察值数目，dt表示在时刻t的某个蛋白质节点的度。

3.2 基于HMM 的复合物识别建模

式中，aij表示在时刻t的状态为Ci，时刻t+1 的状态为Cj的概率；为状态j的观察概率矩阵，(bj(k))，表示状态j输出相应观察值的概率，其中：

定义7（核心节点）给定任意的蛋白质复合物Ct，设W(t,i)表示蛋白质t和蛋白质i相连的边的权值，则选取复合物Ct中具有最大权值之和的蛋白质节点为核心节点，即有：

其中，r表示复合物Ct中的蛋白质总数；l表示与Ct中某一蛋白质i存在相互作用的蛋白质数目；W(k,i)表示蛋白质k和蛋白质i之间的权值。

定义8（相似度）设时序图GS

i中的蛋白质节点表示为，则中任意两个蛋白质节点和之间的相似度为：

定义10（归属度）设动态蛋白质网络用邻接矩阵A表示，rij为A中的元素，则蛋白质网络中的某一蛋白质节点u归属于复合物Ct的程度可表示为：

定义11（观察概率）给定蛋白质网络中任意蛋白质节点u与复合物Ct之间的归属度M(u,Ct),t∈(1,N)，观察概率为：

3.3 蛋白质复合物识别算法

本文提出的蛋白质复合物识别算法充分考虑了动态蛋白质网络中前一时刻网络拓扑对当前时刻网络结构的影响以及网络中节点和复合物的固有特点，建立HMM 来解决动态蛋白质网络中的复合物识别问题。在当前时刻受到前一时刻网络结构影响的前提下，基于HMM 描述蛋白质复合物与网络个体间的相互关系，采用状态链和观察链分别表示蛋白质网络中的复合物和节点信息，通过求解HMM 中的最优状态序列得到动态蛋白质网络中的复合物。蛋白质复合物识别的具体过程如算法2 所示。

算法2基于HMM 的蛋白质复合物识别（HMMPC）

输入：动态蛋白质网络G=

输出：蛋白质复合物。

1.基于IPC-MCE 算法[16]对进行蛋白质复合物的识别，得到初始的蛋白质复合物C-ini；

3.在时刻t(2 ≤t≤T)，对构建矩阵，并根据式（10）计算得到每个蛋白质复合物的核心节点

5.根据式（12）计算时刻t的各个状态转移概率

6.根据式（13）计算时刻t的每个蛋白质与各个复合物之间的归属度；

7.根据式（14）计算时刻t的每个蛋白质与各个复合物之间的观察概率；

8.采用维特比算法[17]来得到最优状态序列，即中的蛋白质复合物；

9.令t=t+1，转3。直到t>12 时，算法结束。

为了便于理解，下面给出关于算法2的细节说明：

（1）在步骤1 中，IPC-MCE 算法是一种基于极大团扩展的蛋白质复合物识别算法，主要应用于静态蛋白质网络中。该算法首先找到蛋白质网络中存在的极大团（即全连通图），并将其看作蛋白质复合物的核心，然后通过定义核心与周围邻居节点的作用概率来对极大团进行扩展，进而识别出更多具有生物意义的蛋白质复合物，该算法对输入参数不敏感，且易于实现。

（2）在步骤2 中，假定任意的蛋白质节点其初始状态概率分布为随机分布，即认为该节点在初始时刻进入到某个蛋白复合网的概率是随机的。

（3）在步骤3 中，构建矩阵的目的是为了描述出蛋白质节点与复合物之间的关系。具体而言，对于来说，矩阵中第i行第j列的元素取值为：

3.4 算法复杂度分析

如算法2 所示，假设n为蛋白质子网中蛋白质节点的数目，m为蛋白质子网中边的数目，N为初始蛋白质复合物的数目。首先，根据文献[16]可知，在步骤1 中采用IPC-MCE 算法识别初始蛋白质复合物C-ini的复杂度为O(nmμ+μ2d2)，其中μ是规模大于等于3 的全连通图的个数，d是最大的稠密子图规模；步骤2～步骤3 的复杂度为O(n)；步骤4 针对每个蛋白质复合物进行操作，但最终仍是对每个蛋白质节点进行处理，其复杂度仍为O(n)；步骤5～步骤7 的复杂度为O(n×N)；步骤8 中的Viterbi 算法采用递归调用来得到最优序列，其复杂度为O(n2×lgn)；步骤9 中对步骤3～步骤8 中的操作重复进行11 次，该步骤的复杂度为一常数。综上所述，本文提出的蛋白质复合物识别算法（HMM-PC）的复杂度为O(nmμ+μ2d2)+O(n)+O(n×N)+O(n2×lgn)，它属于多项式时间，这表明HMM-PC 算法完全可以高效地识别大规模蛋白质网络中的复合物。

4 实验

用Python 语言实现了HMM-PC 算法，在一台8核16 线程的计算机上进行了实验。其中，处理器为Intel CoreTMi5-1035G1 CPU@1.00 GHz，内存为8 GB，操作系统为64 位的Windows10 家庭中文版。为了验证HMM-PC 的有效性，在多个数据集上将HMM-PC算法与目前较为典型的蛋白质复合物识别算法IPCTFC&ACC[8]、DCA[9]、IPC-DGD[10]和IFPA[11]进行了性能比较。为了保证不同算法之间识别结果对比的客观性和有效性，这些典型算法中涉及到的相关参数沿用原文作者在论文实验中的设置，在此不做任何修改。

4.1 实验数据集

本文采用了DIP 数据集、MIPS 数据集[18]和CYC2008 数据集[19]作为测试对象。其中，DIP 数据使用的是DIP20170205 版本。去掉自相互作用和重复相互作用后，该网络中有4 995 个蛋白质和21 554 条边，如图3 所示；MIPS 数据集在去除了其中包含的自相互作用和冗余的相互作用后，最终的相互作用网络包括4 546 个酵母蛋白质和12 319 对可靠的相互作用。CYC2008 数据中则包含了408 个通过生物方法预测到的蛋白质复合物，每个复合物包含两个或两个以上的蛋白质，可用于评估挖掘的蛋白质复合物。

4.2 评价指标

为了综合衡量HMM-PC 算法的优越性，本文采用如下的几种指标来评价多种不同的蛋白质复合物识别算法的性能：

（1）查全率（Recall）、查准率（Precision）和F-measure值的计算公式如下：

Fig.3 Example of protein network in DIP图3 DIP 中的蛋白质网络示例

其中，ER表示HMM-PC 算法识别的蛋白质复合物；RR表示通过生物实验技术测得的蛋白质复合物；MNM(ER,RR)表示ER和RR之间的最大匹配数目。综合考虑查全率和查准率两方面，可得F-measure 的计算公式为：

（2）功能富集程度。文中通过计算蛋白质复合物的P-value 值来衡量它对某个功能的富集程度，其计算公式[20]为：

其中，N表示蛋白质网络的规模；PN表示蛋白质复合物中的蛋白质数量；λ表示蛋白质复合物中具有某项功能的蛋白质数量；F表示PPIN 中具有该功能的蛋白质数量。一般而言，P-value 值越小，则表明识别的蛋白质复合物是随机具有这种功能的概率越低，即蛋白质复合物更具有生物学意义。

4.3 实验结果与分析

4.3.1 HMM-PC 算法与其他算法的比较

为了准确地分析HMM-PC 算法在文中构建的动态蛋白质网络上识别复合物的性能，将HMM-PC 算法与IPC-TFC&ACC 算法、DCA 算法、IPC-DGD 算法和IFPA 算法在DIP 数据集和MIPS 数据集上进行了性能比较。文中对每一种算法都在所构建的每个动态子网上进行复合物识别，对于最终识别得到的复合物都采用文献[19]的后处理方式进行过滤，去掉重复的、包含蛋白质数量少于2 的复合物。表1 和表2分别列出了各种算法在DIP 数据集和MIPS 数据集上的实验结果。

Table 1 Performance comparison of each algorithm on DIP data set表1 各个算法在DIP 数据集上的性能比较

Table 2 Performance comparison of each algorithm on MIPS data set表2 各个算法在MIPS 数据集上的性能比较

从表1 和表2 的结果可以看到，HMM-PC 算法在两种数据集上的查全率和查准率都要优于另外的四种算法。在DIP数据集上，HMM-PC算法的F-measure值要比IPC-TFC&ACC、DCA、IPC-DGD 和IFPA 分别高约43%、32%、23%和11%。在MIPS 数据集上，HMM-PC 算法的F-measure 值要比IPC-TFC&ACC、DCA、IPC-DGD 和IFPA 分别高约41%、35%、28%和14%。这主要是因为：（1）文中在采用基因表达数据构建12 个动态蛋白质子网的基础上，进一步使用了共享功能注释、共享结构域和连接强度等因子来对蛋白质网络进行加权，更为准确地衡量了蛋白质之间的真实相互作用关系；（2）HMM-PC 算法将蛋白质复合物的识别问题转为隐马尔科夫模型中的最优序列发现问题，通过蛋白质复合物核心节点的计算、状态转移概率、观察概率和归属度等的计算，有效地刻画了蛋白质在其整个动态变化生命周期过程中的归属情况，降低了噪声数据和网络动态变化对复合物识别带来的影响，因此取得了比其他算法更优的识别结果。

4.3.2 鲁棒性分析

下面以DIP 数据集为例来测试HMM-PC 算法的鲁棒性。在实验中通过随机地增加和删除蛋白质网络中边的数目来模拟DIP 数据中包含的假阳性和假阴性，实验结果如图4 所示。其中，图4 左边的Y轴表示在增加边的比例（假阳性）情况下测得的Fmeasure 值；右边的Y轴表示在删除边的比例（假阴性）情况下测得的F-measure 值。从图4 可以看到，随着数据中包含的假阳性和假阴性的增加，HMM-PC算法的F-measure 值虽然有微弱的下降，但下降幅度不超过2%，这表明HMM-PC 算法具有较好的抗噪能力，鲁棒性较强。

Fig.4 Robustness analysis of HMM-PC algorithm图4 HMM-PC 算法的鲁棒性分析

4.3.3 HMM-PC 算法的功能富集分析

为了进一步验证HMM-PC 算法识别的蛋白质复合物是否具有生物学意义，借助于TermFinder[21]工具（https://omictools.com/go-termfinder-tool）来计算文中构建的12 个时刻蛋白质子网中的蛋白质复合物的Pvalue 值。为了保证对比的公平性，仅考虑至少包含10个蛋白质的蛋白质复合物。表3分别给出了五种不同算法识别的蛋白质复合物的P-value 在[0，1E-10]，(1E-10，1E-5]，(1E-5，1E-2]和(1E-2，1E+0]区间分布的百分比情况。从表3 可以看到，IPC-TFC&ACC 算法、DCA 算法、IPC-DGD 算法和IFPA 算法识别的蛋白质复合物位于区间[0，1E-5]的比例分别为85.4%、75.1%、79.9%和87.7%，而HMM-PC 算法识别的蛋白质复合物位于区间[0，1E-5]的比例高达93%。由此可见，与其他四种算法相比，通过HMM-PC 算法识别的蛋白质复合物在生物意义方面具有明显的优势。

Table 3 Functional enrichment analysis of protein complexes identified by different algorithms表3 不同算法识别的蛋白质复合物的功能富集性分析

4.3.4 HMM-PC 算法的效率分析

最后，为了进一步衡量HMM-PC 算法的优越性，对各种蛋白质复合物识别的运行时间在DIP 数据集和MIPS 数据集上进行了实验对比，结果如图5 所示。

Fig.5 Comparison of running time of each algorithm图5 各个算法的运行时间比较

从图5 可以看到，IPC-DGD 算法的运行时间最长，DCA 算法次之。而本文提出的HMM-PC 算法在两种数据集上的运行时间大多要远远低于其他四种算法，仅有IFPA 算法在DIP 数据集上的效率略高于HMM-PC 算法。这主要是因为HMM-PC 算法采用动态规划的思路来识别复合物，它将蛋白质复合物问题建模为最优序列发现问题，通过多种概率的计算来确定蛋白质所属复合物的归属，避免了大量重复计算。此外，HMM-PC 算法在运行过程中没有涉及到额外参数的反复调整，使得算法的计算量较低，并使用整个序列的上下文来做判断，从而对包含“噪音”的序列也能进行良好的分析。总的来说，HMMPC 算法具有不错的运行效率，可以快速准确地识别蛋白质复合物，可以推广应用到解决其他复杂网络中的社区发现、子图挖掘等问题，具有普适意义。

5 结束语

在动态蛋白质网络中如何准确且高效地识别蛋白质复合物是目前生物信息学研究中的热点之一。文中首先利用蛋白质的基因表达值、共享功能注释信息、共享结构域信息和连接强度来构建加权的动态蛋白质网络；然后在此基础上将蛋白质复合物识别问题建模为基于隐马尔科夫模型的最优序列发现问题，并采用维特比算法进行求解；最后在多个生物数据集上验证了本文算法的有效性。在下一步工作中，将对未知蛋白质的功能进行预测研究，考虑到图卷积理论在处理大图数据上的优势，拟提出一种基于图卷积神经网络的蛋白质功能预测算法。