一種拓撲與生物功能一致的多網絡比對算法

夏金芳，陳 璟，2

1(江南大學 人工智能與計算機學院，江蘇 無錫 214122) 2(江南大學 物聯網技術應用教育部工程研究中心，江蘇 無錫 214122)

1 引 言

隨著高通量技術的發展，可用的生物數據顯著增加，人們將其建模成蛋白質相互作用(protein-protein interaction，PPI)網絡、基因調控網絡、代謝網絡和疾病網絡，并對其展開深入的研究與分析.網絡比對是一種常用的比較分析網絡的方法，通過對這些網絡的分析可以獲得更多新穎的見解，例如識別可能具有重要功能的進化保守復合物；重構物種的進化動力學；推斷系統發育樹；預測未知蛋白質的功能；驗證已知功能的蛋白質；分析網絡之間的差異等.網絡比對與子圖同構問題有關，后者力求找到一種映射關系使得一個網絡是另一個網絡的精確子圖，這是一個NP-hard問題[1].而網絡比對并不要求是精確子圖，它只需要找到一種最佳的映射關系，用以發現給定網絡之間具有相似拓撲或功能區域的節點[2].同時，由于蛋白質是細胞生命中最活躍的參與者，因此當前的網絡比對廣泛應用于PPI網絡，用以挖掘更多的生物學知識.

目前，研究者們提出了很多比對算法[3，4].根據映射關系可以分為局部比對和全局比對，也可根據比對網絡的數量分為兩個網絡的比對和多個網絡的比對.局部比對的目標是發現小的、保守的區域，但它往往犧牲了網絡的整體相似性，而全局網絡比對與之相反，它旨在找到使網絡整體相似性最大化的保守區域[5].近年來，多網絡比對受到更多的關注，雖然它與兩個網絡比對相比，其復雜度更高，但可以捕獲更多的網絡保守區域.IsoRankN[6]、SMETANA[7]、CSRW[8]、NetCoffee[9]、NetCoffee2[10]、BEAMS[11]、Fuse[12]等都是全局多網絡比對算法，其中，NetCoffee是一對一的多網絡比對，其它算法則是多對多的多網絡比對.IsoRankN將迭代譜聚類融入到IsoRank，并擴展到多個網絡比對.SMETANA基于半馬爾可夫隨機游走模型實現網絡比對，CSRW在SMETANA的基礎上使用上下文相關的隨機游走模型完成比對.NetCoffee使用類似于T-Coffee的方法，將多網絡構造成一組加權二分圖，并使用模擬退火算法來最大化相似性得分.WebNetCoffee[13]為NetCoffee提供了一個友好的用戶界面.Netcoffee2改進了NetCoffee的相似性函數，彌補了NetCoffee只能比對兩個網絡的局限性，并且在運行速度方面也有所改進.BEAMS采用backbone提取和合并方法進行比對.Fuse通過非負矩陣三因式分解(Non-negative Matrix Tri-Factorization，NMTF)方法計算蛋白質相似性得分，然后使用最大權重k分圖匹配算法進行比對.

上述算法中，SMETANA在計算有效性、比對準確性方面表現優異.它通過網絡內和網絡間的概率一致性轉移，考慮了網絡的局部和全局相似性，增強了節點比對概率，然后根據轉移后的概率得分貪心地比對節點.但是，由于蛋白質序列相似性得分高低不能直接確定兩個蛋白質功能是否相同，于是本文提出ConAlign算法，在SMETANA算法的基礎上進行三點改進，以增強節點對的相似性并增加算法使用的簡易性.

1)使用低維向量表示節點的穩態轉移概率，通過節點對之間的穩態概率轉移差異對序列相似性進行懲罰；

2)使用網絡度相似性作為網絡間的整體拓撲相似性；

3)改進網絡比對約束條件，刪除指定參數nmax.

實驗表明，ConAlign算法更高效且保證了拓撲與生物功能的一致性.

2 方 法

2.1 問題定義

PPI網絡可以通過網絡圖的形式表示，本文將PPI網絡建模為無權無向圖，其中節點代表蛋白質，邊代表蛋白質之間的相互作用.假設給定兩個網絡G1(V1,E1)和G2(V2,E2)，其中V1和V2分別是兩個網絡的節點集合，E1和E2分別是兩個網絡的邊集合.兩個網絡比對是尋找一種映射關系f，如式(1)所示，將G1的節點比對到G2的節點，使得比對的得分函數最大，通常得分函數可以是兩個網絡中所有比對上的節點對的相似性之和.

f(u)={v|u∈V1,v∈V2}

(1)

多網絡比對將兩個網絡擴展到多個網絡進行比對.假設給定N個網絡G=(G1,G2,…,GN)待比對，Gi=(Vi,Ei),1≤i≤N表示G中第i個網絡，其中Vi和Ei分別是Gi網絡的節點和邊的集合.本文使用SMETANA中的最大預期準確性(Maximum Expected Accuracy，MEA)作為比對的得分函數.預期準確性(Expected Accuracy，EA)是通過比對簇中節點比對的后驗概率來估算節點正確比對的概率，如公式(2)所示.

(2)

式(2)中EA(A)表示比對A的預期準確性，A表示比對集合，P(ui～vj|G)表示網絡G中ui比對到vj的后驗概率.

2.2 網絡節點拓撲相似性

Connor等人[14]曾指出，如BLAST產生的bit-score等序列相似性值并不能準確度量蛋白質的功能相似性程度，因此本文通過網絡節點的拓撲相似性彌補序列相似性是有必要的.

SMETANA中使用半馬爾可夫隨機游走計算節點對的全局對應得分，無法體現拓撲與生物一致性高的節點對，ConAlign使用網絡節點拓撲差異性作為懲罰項，更好地突出拓撲與生物一致的節點對.ConAlign首先用低維向量π(Gi)表示Gi網絡中節點的穩態轉移概率，然后使用(ui,vj)節點對的穩態轉移概率的相似性作為其拓撲相似性T(ui,vj).節點相似性TS(ui,vj)在拓撲相似性的基礎上融合序列相似性.

(3)

(4)

式(3)中π(ui)表示節點ui簡單隨機游走的穩態分布，式(4)中S(ui,vj)是節點ui和vj的序列相似性，當π(ui)=π(vj)時，直接使用序列相似性作為網絡節點的相似性.

(5)

式中S-T(ui′,vj)表示帶歸一化的節點對得分.

2.3 網絡整體拓撲相似性

SMETANA在網絡間概率一致性轉移時，使用兩個網絡比對的后驗概率作為網絡的整體相似性，該整體相似性受前期相似性計算的影響較大，且兩個網絡比對也需要一定時間.為了防止假陽性，并加快算法效率，ConAlign通過網絡整體拓撲相似性進行網絡間概率一致性轉移.

D(Gi)表示Gi網絡中所有節點的度的集合，D(Gi)={deg(u)|u∈Vi}，(deg(u)表示u節點的度，N(u)表示u節點的鄰居節點的集合，deg(u)=|N(u)|.)由于網絡規模不同，為表現出網絡節點在網絡中的重要性，對網絡中所有節點的度通過式(6)歸一化.

(6)

對于Gi=(Vi,Ei)和Gj=(Vj,Ej)兩個網絡，假設Vi={u1,u2,…,um}，Vj={v1,v2,…,vn}，若m≠n，則通過添加度為0的節點，使得m=n.Gi與Gj之間的距離Dis(Gi,Gj)和相似性Sim(Gi,Gj)的計算分別如式(7)和式(8)所示.

Dis(Gi,Gj)=‖ND(Gi)-ND(Gj)‖F

(7)

(8)

2.4 網絡構建約束條件

SMETANA比對構建過程中的約束條件之一：任何比對簇中來自同一個網絡的節點不能超過nmax，由于nmax參數的選擇沒有統一的標準，且選擇范圍廣，故ConAlign中對此約束進行改進.改進后的比對構建首先對通過概率一致性轉移的估計概率P降序排列，然后按P值從大到小的順序依次分3種情況進行比對.

1)如果當前節點對(ui,vj)都未被比對上，則產生一個新簇，并將這個新簇放入比對集合中；

2)如果當前節點對(ui,vj)有且僅有一個節點被比對上，假設ui已經被比對在簇Cl中，另一個節點vj沒有被比對上.判斷vj能否加入到簇Cl中的流程如圖1所示：輸入P(ui～vj)，首先判斷P(ui～vj)是否滿足公式(9)：大于等于簇Cl中來自i網絡和j網絡的節點對之間的最大值，若滿足，則將vj加入簇Cl中，否則，使用公式(10)計算未比對節點vj屬于該簇Cl的概率，式中k是簇Cl中和vj來自同一個網絡的節點數，β是一個懲罰因子，若滿足條件，則將vj加入簇Cl中，否則丟棄vj.

圖1 比對約束的流程圖Fig.1 Flow chart of alignment constraints

(9)

(10)

3)如果當前節點對(ui,vj)都已經被比對上，但未被比對在同一個簇中，則嘗試將這兩個簇合并，若合并后的簇中來自不同網絡的節點至多一個，則將這兩個簇合并，否則，不合并這兩個簇.

2.5 算法描述

定義 1.若一個節點與其它網絡中的任一節點都沒有序列相似性，則認為這個節點是非同源節點.

算法1.ConAlign算法

輸入：PPI網絡G=(G1,G2,…,GN)，網絡間所有蛋白質對的序列相似性；

輸出：比對的簇結果；

Begin

1.網絡初始化，去除非同源節點；

2.使用公式(4)網絡節點的拓撲相似性融合序列相似性計算節

點相似性；

3.網絡內部概率一致性轉移;

4.網絡間概率一致性轉移，其中網絡間的相似性使用公式(8)；

5.如2.4進行網絡構建；

End

ConAlign算法的復雜度約為O(N3z|Vmax|)，其中，N是要比對的網絡個數，z表示兩個網絡相似性的非零元素最多的個數，Vmax表示最大網絡的節點數.

3 實驗及分析

3.1 實驗數據集

本文分別在合成網絡和真實網絡上測試了ConAlign算法.合成網絡是NAPAbench[15]中的PPI網絡，這三個合成網絡分別對晶體生長(crystal growth，CG)、復制-突變-互補(duplication-mutation-complementation，DMC)、帶有隨機突變的復制(duplication with random mutation，DMR)這三種類型的8-way網絡進行了建模.CG中8個網絡均由1000個節點和3985條邊組成；DMC中8個網絡的節點均是1000，邊的數量分別是1919、1853、1923、1840、1867、1848、1818、1867；DMR中8個網絡的節點均是1000，邊的數量分別是2031、2092、1967、1977、1959、1998、2030、2056.5個真實網絡分別是蠕蟲(Caenorhabditis Elegans，CE)、蠅類(Drosophila Melanogaster，DM)、人類(Homo Sapiens，HS)、老鼠(Mus Musculus，MM)和酵母(Saccharomyces Cerevisiae，SC).真實網絡的所有數據從IsoBase[16]數據庫獲得，網絡的節點和邊信息如表1所示.生物相似性使用BLAST++工具計算得出蛋白質之間的序列相似性得分.

表1 真實網絡信息Table 1 Network information of five real networks

3.2 實驗結果與評估

目前已存在多種評估多網絡比對的質量的方法，如簇和蛋白質的覆蓋量、拓撲和生物學度量指標等.

簇覆蓋量(Cluster coverage)表示簇的數量，即通過比對產生的簇的數量.蛋白質覆蓋量(Protein coverage)表示比對上的蛋白質的數量，即所有簇中蛋白質的總數.簇的k覆蓋量是指包含來自k個物種的蛋白質的簇的數量.

為了進一步度量比對結果，可以從拓撲角度度量比對結果，例如簇相互作用質量(Cluster Interaction Quality，CIQ)[11].也可以從生物學角度度量比對結果簇中蛋白質的功能相似性，例如平均歸一化熵(Mean Normalized Entropy，MNE)[17]，敏感性(Sensitivity，Sen)[11].

CIQ是對所有簇中節點對之間的一種保守相互作用度量[11].

(11)

式(11)中，EClm,Cln是簇Clm和簇Cln之間邊的數量，cs(m,n)的定義如公式(12)：

(12)

式(12)中，sm,n′是簇Clm和簇Cln中邊的網絡的數量，sm,n是簇Clm和簇Cln中存在共享節點的網絡的數量.

MNE是一種一致性度量.歸一化熵(Normalized Entropy, NE)：

(13)

式(13)中其中d是簇Cl中不同基因本體(Gene Ontology,GO)[18]注釋的數量，pi是簇Cl中有GOi注釋的蛋白質的比例.MNE是所有注釋簇的平均歸一化熵，其值越小，表明簇結果越一致.

Sen是在所有GO類別里，最接近簇中與比對節點中被GO類別注釋的蛋白質數量的比例的平均值.對于給定的GO類別，其最接近的簇是指包含最多被該GO類別注釋的蛋白質數量的簇[11].

實驗將本文算法ConAlign與IsoRankN、BEAMS和SMETANAE進行對比.ConAlign算法中α=0.9,β=0.8，與SMETANA原文保持一致，其它三種算法采用原文算法默認參數，IsoRankN：α=0.6；BEAMS在合成網絡上α=0.5,β=0.2；在真實網絡上α=0.5,β=0.4；SMETANA：α=0.9,β=0.8,nmax=10.

3.2.1 合成網絡實驗

圖2是本文算法ConAlign及另三種主流算法IsoRankN、BEAMS、SMETANA在NAPAbench合成網絡上的簇覆蓋量.雖然BEAMS在三個合成網絡上的簇覆蓋量最多，但許多簇是由2或3個物種的蛋白質組成.當產生的簇包含來自所有物種的節點時，才有可能提供所有網絡蛋白質之間的同源信息.ConAlign產生的k=8的簇覆蓋量在DMC網絡上與BEAMS相近，在CG和DMR網絡上比BEAMS多.同時，在這三個合成網絡中，ConAlign產生的k=8的簇占總的簇覆蓋量的比例均高于BEAMS.這表明ConAlign產生的簇更具有研究價值.圖3是8-way合成網絡的蛋白質覆蓋量，ConAlign產生的蛋白質覆蓋量在合成網絡上表現良好，大多數蛋白質均被比對上，即識別出了最多的同源物.

圖2 不同算法在合成網絡上的簇覆蓋量Fig.2 Cluster coverage of different algorithms on synthetic networks

圖3 不同算法在合成網絡上的蛋白質覆蓋量Fig.3 Protein coverage of different algorithms on synthetic networks

CG、DMC、DMR這三個合成網絡比對結果的拓撲和生物性能的分析分別如表2～表4所示.總體上，ConAlign在拓撲指標CIQ和生物指標MNE、Sen上的表現優異.雖然IsoRankN在CG網絡上Sen值最高，但其CIQ值和MNE值均是最差的.對其余情況而言，ConAlign均表現最佳，尤其是在DMR網絡上，ConAlign與次優的SMETANA相比，其CIQ提高了4%、MNE優化了2.6%、Sen提高了10.9%.這表明ConAlign通過拓撲相似性對序列相似性的增強是有效的.

表2 不同算法在CG網絡中的表現Table 2 Performance of different algorithms on CG networks

表3 不同算法在DMC網絡中的表現Table 3 Performance of different algorithms on DMC networks

本算法在所有合成網絡上的CIQ值都最高，表明比對結果中保守相互作用多，且在大部分情況下，比對結果的MNE值較低且Sen值較高，表明ConAlign產生的比對結果更具有功能相關性.

表4 不同算法在DMR網絡中的表現Table 4 Performance of different algorithms on DMR networks

3.2.2 真實網絡實驗

為了進一步評估本算法，本文將四種算法在真實網絡上也進行了實驗.如圖4和圖5所示，不同算法在真實網絡上的簇覆蓋量和蛋白質覆蓋量與在合成網絡上呈現的趨勢近似，簇覆蓋量之間的差異主要體現在k=2和k=3，但它們不能像k=5的簇一樣提供所有物種間蛋白質之間的關系.ConAlign產生最多的蛋白質覆蓋量，也就是說比對上的蛋白質是最多的.

圖4 不同算法在真實網絡中產生的簇覆蓋量Fig.4 Cluster coverage of different algorithms on real networks

圖5 不同算法在真實網絡中產生的蛋白質覆蓋量Fig.5 Protein coverage of different algorithms on real networks

表5 不同算法在真實網絡上的表現Table 5 Performance of different algorithms on real networks

真實網絡比對結果的拓撲和生物性能分析如表5所示，雖然BEAMS產生的比對結果有最小的MNE值，但是ConAlign產生的結果有著高的CIQ值和Sen值，且其MNE值僅次于BEAMS.

圖6是真實網絡中不同算法分別對CIQ、Sen、MNE的效果從好到壞降序排列，將同種算法用直線相連.交叉的直線表示該算法在一個度量指標上表現良好，但在另一個度量指標上表現不佳.由圖6可知，BEAMS在生物方面的MNE表現優異，但其在拓撲方面表現沒有SMETANA和ConAlign好，而ConAlign是未交叉的直線中排名最好的，其產生的結果在拓撲和生物方面達到了平衡，具有拓撲與生物功能一致性.

圖6 不同算法在真實網絡中度量指標的關系Fig.6 Relationship of indicators between different algorithms on real networks

在真實網絡中，ConAlign在生物指標MNE上的表現未達到最佳，可能的原因有兩點：一是，ConAlign為了防止假陽性，采用網絡整體拓撲相似性作為網絡的相似性，并沒有融入序列相似性.雖然，在計算網絡間整體相似性時犧牲了序列相似性，但是，運行速度得到了提升，如表6所示，ConAlign在合成網絡和真實網絡上運行的時間比其它三個算法都短，而BEAMS的運行時間是最長的.二是，目前的蛋白質的功能注釋可能不完善，因此，很難完全可靠評估預測比對的生物性能.

表6 不同算法在不同網絡上的運行時間Table 6 Runtimes of different algorithms on different networks

ConAlign算法在合成網絡和真實網絡上的實驗結果表明，其比對結果在拓撲指標和生物指標這兩方面達到了整體最優，體現了拓撲與生物功能一致性.

4 結 論

本文提出了一種全局多對多的多網絡比對算法ConAlign，該算法在SMETANA算法的基礎上進行了改進，首先，使用節點對拓撲相似性補充序列相似性，并結合網絡整體拓撲相似性進行網絡間概率一致性轉移.其次，通過修改比對約束條件，刪除了難以確定的參數.本算法分別在合成網絡和真實網絡上進行實驗，并與IsoRankN、BEAMS、SMETANA算法進行比較.實驗表明，ConAlign是一種運行速度快，參數少，蛋白質覆蓋量高且具有拓撲與生物功能一致性的比對算法.

為了保證生物網絡比對具有高的生物相似性，本算法是在序列相似性的基礎上完成.由于序列相似性數據的不完整性，在接下來的工作中，將嘗試僅用拓撲相似性來完成網絡比對.