李 鵬,羅愛靜,閔 慧,譚蓀怡,郭惠敏
1.中南大學 湘雅三醫院,長沙 410006
2.湖南中醫藥大學 信息科學與工程學院,長沙 410208
3.醫學信息研究湖南省普通高等學校重點實驗室(中南大學),長沙 410006
4.湖南信息職業技術學院 軟件學院,長沙 410200
隨著人類基因組計劃以及對于多個物種全基因組測序工作的結束,目前生命科學研究的重點已經轉移到了蛋白組學[1]。蛋白質是生物體生命活動的重要物質基礎,沒有蛋白質就沒有生命。一個生物體內所有蛋白質相互作用被稱為蛋白質相互作用網絡(protein-protein interaction networks,PPIN),簡稱蛋白質網絡[2]。值得注意的是,蛋白質之間的相互作用是動態的,它會隨著時間環境、蛋白質的存在和降解、細胞的不同生理狀態等因素的變化而變化,但由于PPIN 本身的復雜性、可利用的蛋白質相互作用數據的不完全性和存在噪聲等眾多問題,準確且高效地衡量蛋白質相互作用的動態性還存在很多挑戰[3-4],這也直接限制了PPIN 領域內很多問題[5-6](復合物識別、關鍵蛋白識別、網絡比對等)的研究進展。為此,本文將對動態蛋白質網絡的構建問題進行研究,并在此基礎上提出新的蛋白質復合物識別方法,這有利于從微觀層面解釋細胞內蛋白質之間的復雜關系,為生物學家、醫學家理解生命復雜網絡的內在組織和生物過程提供新的途徑,在藥物標靶設計、疾病診治和預測等方面也具有重要的應用價值。
目前,已有眾多研究者對動態蛋白質網絡中的復合物識別問題進行了研究,例如雷秀娟等人[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 算法的有效性。
蛋白質網絡具有偏好連接特性,且隨著網絡規模的增加,不同節點間的連接數差異呈現著增加的趨勢。因此本文通過構建加權的動態蛋白質網絡來反映真實網絡結構,在此基礎之上識別蛋白質復合物。
一般而言,蛋白質網絡可用無向圖進行建模,設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)確定鄰接矩陣中對應位置的元素的值,從而可以得到一個無向圖來表示蛋白質網絡。
通過基因表達數據構建得到的初始蛋白質網絡只能衡量節點之間的連接關系,但不能反映出節點間的偏好連接特性。動態蛋白質網絡可由多個時刻的靜態蛋白質網絡表示,本文將動態蛋白質網絡按照時間粒度劃分為多個網絡快照。有如下的定義:
定義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 個時刻的網絡都賦權完畢,則算法結束。
隱馬爾科夫模型(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的某個蛋白質節點的度。

式中,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),觀察概率為:

本文提出的蛋白質復合物識別算法充分考慮了動態蛋白質網絡中前一時刻網絡拓撲對當前時刻網絡結構的影響以及網絡中節點和復合物的固有特點,建立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列的元素取值為:

如算法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 算法完全可以高效地識別大規模蛋白質網絡中的復合物。
用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]進行了性能比較。為了保證不同算法之間識別結果對比的客觀性和有效性,這些典型算法中涉及到的相關參數沿用原文作者在論文實驗中的設置,在此不做任何修改。
本文采用了DIP 數據集、MIPS 數據集[18]和CYC2008 數據集[19]作為測試對象。其中,DIP 數據使用的是DIP20170205 版本。去掉自相互作用和重復相互作用后,該網絡中有4 995 個蛋白質和21 554 條邊,如圖3 所示;MIPS 數據集在去除了其中包含的自相互作用和冗余的相互作用后,最終的相互作用網絡包括4 546 個酵母蛋白質和12 319 對可靠的相互作用。CYC2008 數據中則包含了408 個通過生物方法預測到的蛋白質復合物,每個復合物包含兩個或兩個以上的蛋白質,可用于評估挖掘的蛋白質復合物。
為了綜合衡量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.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 算法具有不錯的運行效率,可以快速準確地識別蛋白質復合物,可以推廣應用到解決其他復雜網絡中的社區發現、子圖挖掘等問題,具有普適意義。
在動態蛋白質網絡中如何準確且高效地識別蛋白質復合物是目前生物信息學研究中的熱點之一。文中首先利用蛋白質的基因表達值、共享功能注釋信息、共享結構域信息和連接強度來構建加權的動態蛋白質網絡;然后在此基礎上將蛋白質復合物識別問題建模為基于隱馬爾科夫模型的最優序列發現問題,并采用維特比算法進行求解;最后在多個生物數據集上驗證了本文算法的有效性。在下一步工作中,將對未知蛋白質的功能進行預測研究,考慮到圖卷積理論在處理大圖數據上的優勢,擬提出一種基于圖卷積神經網絡的蛋白質功能預測算法。