生物復(fù)雜網(wǎng)絡(luò)motif發(fā)現(xiàn)的并行算法

楊伏長 朱嘉富 孫佳敏 謝江

摘 要：生物復(fù)雜網(wǎng)絡(luò)motif發(fā)現(xiàn)是一種研究生物網(wǎng)絡(luò)的重要方法，它基于復(fù)雜網(wǎng)絡(luò)的理論研究，以新的視角來研究生命現(xiàn)象和生命機制，但是在處理較大的網(wǎng)絡(luò)規(guī)模或者需挖掘較大的motif時計算效率低。針對這個問題，在現(xiàn)有串行網(wǎng)絡(luò)motif發(fā)現(xiàn)算法ESU的基礎(chǔ)上，提出一種基于消息傳遞接口（MPI）的并行化ESU算法。該方法在ESU計算過程中優(yōu)化了節(jié)點值以解決節(jié)點值依賴問題，并以ESU算法的子圖發(fā)現(xiàn)策略統(tǒng)計各節(jié)點子圖數(shù)，利用動態(tài)規(guī)劃策略尋找最佳節(jié)點分配策略以解決負載不均衡問題。模擬網(wǎng)絡(luò)數(shù)據(jù)和真實生物網(wǎng)絡(luò)數(shù)據(jù)的實驗結(jié)果表明，并行化ESU算法優(yōu)化了節(jié)點值依賴問題，實現(xiàn)了基于動態(tài)規(guī)劃的負載均衡策略，其運行時間比串行算法縮短了90%，并且該并行算法對不同類型不同規(guī)模的網(wǎng)絡(luò)都具有較強的適用性，有效地提高了網(wǎng)絡(luò)motif發(fā)現(xiàn)問題的計算效率。

關(guān)鍵詞：網(wǎng)絡(luò)motif發(fā)現(xiàn);子圖枚舉;同構(gòu)比較;并行化;消息傳遞接口

中圖分類號： TP301.6

文獻標志碼：A

Abstract： Biological complex network motifs discovery， based on theoretical research of complex networks， is an important method for studying biological networks， which provides a new perspective on life phenomena and life mechanisms. However， it computes inefficiently when dealing with large networks or mining big motifs. On the basis of existing serial ESU （Enumerate SUbgraph） algorithm of network motifs discovery， a parallelized ESU algorithm based on Message Passing Interface （MPI） was proposed. The node values in ESU algorithm were optimized to solve the problem of node value dependency， the number of subgraphs was counted by using subgraph discovery strategy of ESU algorithm， and a dynamic programming method was used to determine optimal node allocation strategy to satisfy load balancing. The experiments on simulated and biological networks show that the parallelized ESU algorithm addresses node value dependency and realizes a load balancing strategy， which saves more than 90% running time compared to serial algorithm. Furthermore， the parallel algorithm is suitable for different types and different scales of networks， and effectively improves computation efficiency of network motifs discovery.

Key words： network motifs discovery; Enumerate SUbgraph （ESU）; homogeneous comparison; parallelization; Message Passing Interface （MPI）

0 引言

生物網(wǎng)絡(luò)比對對網(wǎng)絡(luò)醫(yī)學(xué)和物種演化等領(lǐng)域的發(fā)展有著重大的意義[1]，其中，一項很重要的研究就是生物網(wǎng)絡(luò)motif發(fā)現(xiàn)。生物網(wǎng)絡(luò)motif發(fā)現(xiàn)是指在生物網(wǎng)絡(luò)中尋找所有特定大小、拓撲結(jié)構(gòu)一致的子圖[2]，它包括在生物網(wǎng)絡(luò)中枚舉子圖和計算子圖之間的同構(gòu)關(guān)系等問題。生物網(wǎng)絡(luò)motif是生物網(wǎng)絡(luò)的基本單元，它所表現(xiàn)的局部性質(zhì)與網(wǎng)絡(luò)的構(gòu)成機制密切相關(guān)，對它的研究可以加深對生物網(wǎng)絡(luò)的結(jié)構(gòu)和功能的理解，有助于揭示很多重要的生物學(xué)問題[3];然而精確識別網(wǎng)絡(luò)motif是一個計算復(fù)雜度很高的問題[4]，每一個枚舉的子圖都需要執(zhí)行同構(gòu)比較，而子圖之間同構(gòu)關(guān)系的計算一直是NP（Non-deterministic Polynomial）難問題，無法通過改進算法來降低其時間復(fù)雜度。隨著生物網(wǎng)絡(luò)研究的規(guī)模和motif的大小不斷增大，參與枚舉和同構(gòu)比較的子圖數(shù)呈指數(shù)級別的增長，其計算量越來越大，給網(wǎng)絡(luò)motif的發(fā)現(xiàn)帶來了挑戰(zhàn)。

當前網(wǎng)絡(luò)motif發(fā)現(xiàn)研究算法仍以Wernicke[5]于2006年提出的ESU（Enumerate SUbgraph）算法為主流，其衍生算法如Ribeiro等[6]于2010提出的基于G-Tries的改進算法，通過G-Tries來存儲motif所枚舉的圖，減少了比較子圖同構(gòu)次數(shù)來降低時間、提高性能;Khakabimamaghani等[7]于2013年提出了QuateXelero算法，該算法利用四分樹結(jié)構(gòu)來減少同構(gòu)計算次數(shù);Luo等[8]于2018年提出的LCNM（Large Co-regulatory Network Motifs）算法，通過結(jié)合枚舉和隨機化ESU算法，提高了網(wǎng)絡(luò)motif發(fā)現(xiàn)的計算效率。

目前，已有的網(wǎng)絡(luò)motif發(fā)現(xiàn)算法，在計算規(guī)模上都存在瓶頸[9]。隨著生物網(wǎng)絡(luò)數(shù)據(jù)獲取的渠道越來越多，生物網(wǎng)絡(luò)規(guī)模越來越大，對計算效率的要求也會越來越高，因此，實現(xiàn)高效的并行化網(wǎng)絡(luò)motif發(fā)現(xiàn)算法是很有必要的。本文從串行的ESU算法著手，分析算法中各個部分的耗時，針對算法中耗時最長的部分，以消息傳遞接口（Message Passing Interface， MPI）為基礎(chǔ)實現(xiàn)并行化，并結(jié)合優(yōu)化節(jié)點值和動態(tài)分配節(jié)點策略改進并行算法，最后通過仿真數(shù)據(jù)和真實數(shù)據(jù)進行了分析和討論。

1 網(wǎng)絡(luò)motif發(fā)現(xiàn)

網(wǎng)絡(luò)motif發(fā)現(xiàn)問題本質(zhì)上是統(tǒng)計各種非同構(gòu)子圖在網(wǎng)絡(luò)中出現(xiàn)的頻率，并與隨機網(wǎng)絡(luò)所產(chǎn)生的子圖的頻率進行比較，出現(xiàn)頻率明顯高于隨機網(wǎng)絡(luò)的那些子圖可能會揭示網(wǎng)絡(luò)的重要功能[2]。然而，隨著估計重要性方法的出現(xiàn)，網(wǎng)絡(luò)motif發(fā)現(xiàn)問題的計算負擔就落在原始網(wǎng)絡(luò)的子圖統(tǒng)計上[10]，因此網(wǎng)絡(luò)motif發(fā)現(xiàn)過程就是各種子圖的統(tǒng)計過程，其形式化定義[11]如下：

定義1 網(wǎng)絡(luò)motif發(fā)現(xiàn)。給定一個網(wǎng)絡(luò)G和參數(shù)k，找出所有大小為k的子圖，并根據(jù)各個子圖間的同構(gòu)關(guān)系進行分類。

根據(jù)子圖間的同構(gòu)關(guān)系進行分類是網(wǎng)絡(luò)motif發(fā)現(xiàn)過程中計算量最大的部分，目前尚無多項式級別的算法能夠完成這個分類過程，即使在較高效的算法如nauty算法[12]，其在最壞情況下，計算復(fù)雜度仍達到O（n！）。

2.1 串行ESU算法

2.1.1 算法實現(xiàn)原理

ESU算法是一種枚舉子圖的motif發(fā)現(xiàn)算法，ESU算法可以大致分為網(wǎng)絡(luò)構(gòu)建、枚舉子圖、子圖同構(gòu)比較和統(tǒng)計子圖結(jié)果四個步驟，其中ESU算法所提出的枚舉子圖的過程保證了所枚舉子圖的唯一性和子圖數(shù)的準確性，枚舉子圖的過程是后續(xù)子圖同構(gòu)比較過程的實現(xiàn)前提，因此，下文對ESU算法的枚舉子圖過程展開介紹。

ESU算法的輸入為邊集的形式，節(jié)點的序號（在后文中使用節(jié)點值來表示節(jié)點的序號）以連續(xù)的數(shù)值表示，利用基于節(jié)點值的方式構(gòu)建圖，在不斷構(gòu)建的過程中記錄節(jié)點值比節(jié)點本身大的鄰居，以防止重復(fù)搜索子圖，在搜索子圖的過程中，不斷加入節(jié)點值比自身大的鄰居，并將新加入節(jié)點的符合規(guī)則的鄰居也作為待搜索節(jié)點，以防止遺漏子圖。ESU算法偽代碼[5]描述如下：

其中，輸入用節(jié)點集V和邊集E描述的圖G，和給定一個k用來描述motif的大小，算法將會輸出所有在G中motif大小為k的子圖。算法從V中的每一個節(jié)點v開始，初始時確定VExtension集合，VExtension集合包含了所有節(jié)點值大于v的鄰居節(jié)點，偽代碼中用N來描述鄰居關(guān)系，VSubgraph集合初始化為{v}，是存儲子圖節(jié)點值的集合，然后調(diào)用ExtendSubgraph函數(shù)進行子圖搜索和同構(gòu)比較;ExtendSubgraph函數(shù)中，其關(guān)鍵步驟在于對VExtension中的每一個節(jié)點w進行如下操作：將w從VExtension集合中移出并加入VSubgraph集合，更新VExtension集合，新的VExtension集合由VSubgraph集合的所有節(jié)點值大于v的鄰居組成，偽代碼中用Nexcl來描述這種新的鄰居關(guān)系，每當VSubgraph中的節(jié)點數(shù)達到了指定的k大小，就執(zhí)行網(wǎng)絡(luò)同構(gòu)比較。

ESU算法枚舉子圖的過程如圖1所描述，在圖1中，枚舉過程被描述為樹形結(jié)構(gòu)，這種樹被稱為ESU樹[5]。ESU樹是建立ESU算法正確性的基礎(chǔ)，它能保證所枚舉的每個子圖不重復(fù)，其向下延伸的過程正好對應(yīng)了EnumerateSubgraph函數(shù)。其中，圖G是一個包含9個節(jié)點9條邊的拓撲網(wǎng)絡(luò)，其節(jié)點值依次從1到9。算法從Root開始，枚舉圖G所有存在的子圖，圖1中最后輸出16個葉子節(jié)點，對應(yīng)為所有motif大小為3的子圖。

2.1.2 算法分析

ESU算法中最重要的兩個步驟就是子圖枚舉和子圖同構(gòu)比較。

子圖枚舉是ESU算法的核心步驟，通過記錄節(jié)點值比自身大的鄰居，以防止在搜索過程中得到重復(fù)的子圖，但是，基于節(jié)點值的枚舉過程，會導(dǎo)致各個計算量集中在某些特殊的節(jié)點，比如一些鄰居很多而且鄰居節(jié)點值都比自身大的那些節(jié)點（如圖1中1，2，3節(jié)點），這種問題被稱為節(jié)點值依賴問題，這種問題在串行過程中體現(xiàn)不出來，但是在并行過程中，一旦出現(xiàn)這種問題，就會造成負載不均衡，大幅度影響并行的性能。

子圖同構(gòu)比較是整個串行ESU算法耗時時間最長也是計算復(fù)雜度最高的部分[7]，在串行ESU算法中，子圖同構(gòu)使用了高效的nauty算法[12]，基于子圖枚舉的結(jié)果，對每個子圖進行同構(gòu)比較;然而nauty算法僅僅是提高了兩個子圖之間同構(gòu)比較的計算效率，它并不能解決并行過程中存在的其他問題，如負載不均衡問題，而這些問題正是提高并行效率的關(guān)鍵。

由于每一個枚舉的子圖都進行同構(gòu)比較的計算，枚舉子圖的結(jié)果就直接決定了各個節(jié)點子圖同構(gòu)比較的計算量，所以本文對子圖枚舉的過程進行并行處理，以均衡各個節(jié)點的計算量。

2.2 ESU算法并行化

由于子圖同構(gòu)比較是整個算法最耗時的部分，因此并行處理的關(guān)鍵就是均衡各個節(jié)點所需比較子圖同構(gòu)的次數(shù)。本文就均衡子圖同構(gòu)比較次數(shù)上，提出了結(jié)合優(yōu)化節(jié)點值和動態(tài)規(guī)劃分配節(jié)點的改進方法。

1）優(yōu)化節(jié)點值。由于節(jié)點值依賴問題使得在枚舉子圖的過程中，部分枚舉的子圖集中在少數(shù)的節(jié)點上，直接導(dǎo)致各個節(jié)點所需比較子圖同構(gòu)的次數(shù)差異很大，對并行性能造成很大影響，因此在2.2.1節(jié)中，優(yōu)化了節(jié)點值，改善了節(jié)點值依賴問題。

2）動態(tài)規(guī)劃分配節(jié)點。由于各個節(jié)點出發(fā)枚舉的子圖數(shù)不可能完全一致，在分配任務(wù)時就可能造成各節(jié)點計算量不一致，導(dǎo)致負載不均衡，因此在2.2.2節(jié)中，提出了基于動態(tài)規(guī)劃的節(jié)點分配策略，實現(xiàn)了負載均衡。

2.2.1 解決節(jié)點值依賴問題

由于ESU算法是根據(jù)輸入網(wǎng)絡(luò)的節(jié)點值來記錄相應(yīng)符合條件的鄰居以防止重復(fù)搜索，造成了網(wǎng)絡(luò)各個節(jié)點之間的計算存在一定的依賴關(guān)系。對于網(wǎng)絡(luò)的不同輸入形式，比如同一拓撲結(jié)構(gòu)的網(wǎng)絡(luò)，其節(jié)點賦予不同的節(jié)點值，會導(dǎo)致同一節(jié)點在不同節(jié)點值的情況下計算量差異很大，最典型的便是那些度很高的節(jié)點，如果其鄰居的節(jié)點值都大于它本身會使計算量全部集中在這些度很高的節(jié)點上，就會對并行的性能造成很大的影響，因此本文在網(wǎng)絡(luò)分發(fā)之前進行節(jié)點值優(yōu)化，步驟如下：首先統(tǒng)計各個節(jié)點值在邊集中出現(xiàn)的頻次，然后將節(jié)點值集合從大到小依次賦給出現(xiàn)頻次從大到小的節(jié)點，完成對節(jié)點值的優(yōu)化過程。具體如圖2所示。

對于圖2中兩個同一拓撲結(jié)構(gòu)的網(wǎng)絡(luò)，網(wǎng)絡(luò)A中那些度很大的節(jié)點的節(jié)點值都很小，網(wǎng)絡(luò)B那些度很大的節(jié)點的節(jié)點值都很大，兩種網(wǎng)絡(luò)分別調(diào)用ESU算法進行計算，得出各個節(jié)點出發(fā)所需計算比較同構(gòu)子圖數(shù)如表1所示。

結(jié)合分析表1中的數(shù)據(jù)和圖1的拓撲結(jié)構(gòu)可知，當節(jié)點度的大小與節(jié)點值的大小成正比時，各個節(jié)點出發(fā)的計算量差異相對較小，也就是說節(jié)點值優(yōu)化使得各個節(jié)點的計算量相對均衡。

2.2.2 網(wǎng)絡(luò)節(jié)點分配過程

由前面的分析可知，對ESU算法進行并行化處理的關(guān)鍵在于將搜索入口節(jié)點分發(fā)給各個進程以實現(xiàn)并行計算，而對ESU算法的并行化過程存在兩個很明顯的特點：

1）每個搜索入口頂點進行深度搜索時所需要的圖的區(qū)域是未知的，因此每個搜索入口都需要對整個圖的區(qū)域進行搜索，無法對圖進行區(qū)域分發(fā);

2）每個搜索入口頂點進行深度搜索時所得到的子圖數(shù)是未知的，那么在將所有的子圖進行同構(gòu)比較時，所需的時間也是未知，所有節(jié)點所花費的時間不一，對并行性能會造成很大的影響。

對于特點1）可知，需要將整個圖分發(fā)給所有進程，每個進程都存儲一個圖的內(nèi)存，雖然會占用一定的資源，但是一般的圖占用的內(nèi)存都很小，對于現(xiàn)代內(nèi)存大都大于10GB以上的計算機而言，這點內(nèi)存顯得微不足道。

對于特點2）可知，每個搜索入口節(jié)點所需的計算時間差異很大，對于并行算法而言，以頂點作為分發(fā)手段不是明智的選擇，需要考慮其他的負載平衡策略，當仔細分析算法計算過程時間消耗分布時，發(fā)現(xiàn)與ESU算法中深度搜索得到的子圖數(shù)這個過程相比，將所得到的每個子圖進行同構(gòu)比較的時間要遠遠大于搜索的過程，于是，本文在原有的算法上增加了一步：先從每個搜索入口節(jié)點出發(fā)，深度搜索得到每個節(jié)點所需同構(gòu)比較的子圖數(shù)，接著不進行同構(gòu)比較，而是依據(jù)各節(jié)點統(tǒng)計的子圖數(shù)對所有節(jié)點進行排序，同時使用動態(tài)規(guī)劃將頂點分配到各個進程，使得進程子圖數(shù)總和之間的方差最小。具體步驟描述如下。

1）使用ESU算法的子圖發(fā)現(xiàn)策略統(tǒng)計各節(jié)點子圖數(shù);

2）依據(jù)子圖數(shù)將節(jié)點進行降序排序;

3）確定節(jié)點與進程號之間的關(guān)系：首先將所有子圖數(shù)求和得到SubgraphSum，將其除以進程數(shù)n得到每個進程所需比較同構(gòu)子圖數(shù)的上限AverageSum，確定每個進程初始子圖數(shù)CurSum為0，不斷將已排序的節(jié)點移出并逐個賦給不同進程，并更新CurSum，對于每一個進程，如果CurSum>AverageSum，則不加入這個節(jié)點，并嘗試將節(jié)點加入下一個進程。

2.2.3 ESU算法并行框架

3 實驗與分析

3.1 實驗環(huán)境

該程序運行在上海大學(xué)高性能計算集群自強4000平臺，實驗使用2臺計算節(jié)點，每個計算節(jié)點使用2顆Intel E5-2690 CPU（2.9GHz/8-core），內(nèi)存大小為64GB，集群節(jié)點間使用標準的Clos二層Infiniband網(wǎng)絡(luò)架構(gòu)。實驗運行操作系統(tǒng)為CentOS 6.3，編程語言為C++，MPI庫版本為IntelMPI。

3.2 實驗數(shù)據(jù)

實驗同時使用了模擬數(shù)據(jù)和真實生物網(wǎng)絡(luò)數(shù)據(jù)進行性能分析，其中模擬數(shù)據(jù)使用NetworkX[13]的Python包模擬了三類不同的網(wǎng)絡(luò)模型，分別是無標度網(wǎng)絡(luò)模型[14]、小世界網(wǎng)絡(luò)模型[15]和規(guī)則網(wǎng)絡(luò)模型[16]，每類網(wǎng)絡(luò)模型隨機地改變邊數(shù)和節(jié)點數(shù)，產(chǎn)生具有不同平均節(jié)點度的網(wǎng)絡(luò)，以探討在不同的平均節(jié)點度下并行ESU算法的性能。實驗中使用固定的4000個節(jié)點，同時變化邊的條數(shù)，分別為20000、40000和60000條邊，來產(chǎn)生不同的節(jié)點平均度，分別是5、10和15，每種平均節(jié)點度模型均生成10個網(wǎng)絡(luò);真實生物網(wǎng)絡(luò)則選取了酵母菌代謝網(wǎng)絡(luò)、酵母菌蛋白質(zhì)網(wǎng)絡(luò)和人類蛋白質(zhì)網(wǎng)絡(luò)三個生物網(wǎng)絡(luò)數(shù)據(jù)集，其節(jié)點數(shù)和邊數(shù)如表2所示，其中：酵母菌的蛋白質(zhì)網(wǎng)絡(luò)和代謝網(wǎng)絡(luò)數(shù)據(jù)集來源于Uri Alon實驗室[2]，人類蛋白質(zhì)網(wǎng)絡(luò)數(shù)據(jù)來源于STRING（Search Tool for Recurring Instances of Neighbouring Genes）在線數(shù)據(jù)庫[17]。

3.3 實驗結(jié)果與分析

3.3.1 優(yōu)化節(jié)點值依賴的結(jié)果

為驗證有優(yōu)化節(jié)點值策略和無優(yōu)化節(jié)點值策略對網(wǎng)絡(luò)并行性能的影響，本文對前文提到的3種模擬網(wǎng)絡(luò)進行測試，每種模擬網(wǎng)絡(luò)都在不同的motif大小以及不同的平均節(jié)點度下進行測試，各種條件下的測試結(jié)果都很相似，因此，本文選取其中一種測試條件下的測試結(jié)果進行描述，選取的條件為節(jié)點數(shù)為4000、邊數(shù)為20000、motif大小為4，測試結(jié)果如圖3所示。

從圖3中可以看出，使用優(yōu)化節(jié)點值策略與否對無標度網(wǎng)絡(luò)模型的并行性能影響很大，而對小世界網(wǎng)絡(luò)和規(guī)則網(wǎng)絡(luò)的并行性能影響很小。這是由于小世界網(wǎng)絡(luò)和規(guī)則網(wǎng)絡(luò)的演化規(guī)則與節(jié)點演化的先后順序無關(guān)。小世界網(wǎng)絡(luò)中那些度很大的節(jié)點的節(jié)點值具有隨機性，使用優(yōu)化節(jié)點值策略對并行性能影響不大。規(guī)則網(wǎng)絡(luò)中，各節(jié)點的度一致，其并行性能與使用優(yōu)化節(jié)點值策略與否無關(guān)。無標度網(wǎng)絡(luò)的演化模型恰好是一種節(jié)點值依賴問題的極端情況，最先演化的那些點恰好是度最大且節(jié)點值最小的點，而計算量正是集中在這些最先演化的點上，導(dǎo)致并行性能很差。使用優(yōu)化節(jié)點值策略，對那些最度此句不通順很大且節(jié)點值很小的點的節(jié)點值重新賦值，使各個節(jié)點所枚舉子圖的數(shù)量差異相對較小，從而提高并行性能。在進程數(shù)達到32時，未使用優(yōu)化節(jié)點值策略的加速比為13.5，而使用優(yōu)化節(jié)點值依賴策略的加速比為29.7?？梢妼?jié)點值之間的依賴關(guān)系進行優(yōu)化是非常有必要的。