基于小波變換的代謝網絡比較方法＊

彭佳揚，羅 永，楊路明，李 敏，趙 巖，成 磊

（1.中南大學 信息科學與工程學院，湖南 長沙 410083；2.國防科學技術大學 理學院，湖南 長沙 410073）

隨著多物種的全基因組測序的完成以及基因組注釋信息的日益完善，使得可以從基因組信息較為可靠地重構物種特異的代謝網絡[1－2].代謝網絡的重構使得研究者可以從整體水平上來分析物種間的進化關系.然而物種整體水平上的代謝網絡是由很多代謝通路交聯而成，過于龐大和復雜，從而對整個網絡進行逐個元素的比對，復雜度特別高.越來越多的研究者開始關注用網絡比對方法來確定網絡間的進化距離.

現階段確定代謝網絡間進化距離的策略主要有兩種：分別是基于通路比對的通路水平上的網絡比對和物種水平上的全網絡比對.基于通路比對[3－5]的方法是從每個待研究物種的代謝網絡中選取一條共有的、且相對保守的代謝通路（如三羧酸循環或糖酵解通路），然后通過對每對代謝網絡中的該條代謝通路元素進行逐個比對，來確定通路間的差異，即代表物種間的進化距離.基于全網絡比較的方法通常將物種的整個網絡簡單地看成一個節點集合，然后通過集合論方法，通過計算集合間距離的函數，如Jaccard距離或Hamming距離等，來確定網絡間的進化距離.2004年，Agruilar等人[6]將物種的酶圖表示成酶的集合，并根據物種中酶的出現或缺失進一步將物種特異的酶集合表示成0－1向量，然后通過計算歸一化后的Hamming距離得到向量間距離，最終構建出物種間的進化距離矩陣；2004年，Ma等人[7]將重構的82個物種特異的代謝網絡表示為物種特異的酶集合，然后計算Jaccard系數、Simposon系數和Korbel系數來衡量集合間的距離，以此來確定代謝網絡間的進化距離；2006年，Forst等人[8]將代謝網絡表示為“干凈的”代謝物－反應二部圖，通過計算平凡Jaccard距離這一度量函數來獲得網絡間的進化距離；2006年，Tohsato等人[9]則將物種特異的代謝網絡表示為酶促反應集合，根據平凡Jaccard距離來計算物種特異的酶促反應集合之間的差異，從而得到代謝網絡的進化距離矩陣.2008年，周婷婷等人[10]加入了酶的特征來計算進化距離，提出了一種確定進化距離的WJD模型，并分析和評價了其合理性.

目前的這些計算代謝網絡間進化距離的方法大都是基于集合論的方法.該方法僅僅是將網絡抽象成節點的集合，并根據節點集合的差異來確定不同網絡間的距離，而沒有考慮到節點本身的特性以及網絡拓撲結構屬性方面的差異，因此這類網絡比較方法還有待改進.本文建立了代謝網絡的相關數學模型，主要運用小波時頻分析方法，結合主成分分析法對代謝網絡進行相似性分析.首先用主成分分析法提取代謝網絡節點的7個中心性拓撲屬性參數的主成分，因為主成分是無序的，接下來通過統計的手段對降維的主成分數據進行有序化，最后為了提取有序化后的主成分曲線趨勢，利用小波時頻分析方法提取低頻信號作為代謝網絡的特征曲線，用以比較不同代謝網絡的相似度，由此推算出代謝網絡的進化距離.

1 網絡節點的中心性拓撲參數

在研究復雜網絡時，通常用節點的中心性度量來衡量節點在網絡中影響能力的大小，通過網絡的拓撲屬性可了解該對象獲得、控制信息及資源的能力.在無向網絡中，度中心性測度（簡稱節點度或度）是最常用的測度，此外常用的節點測度還包括聚集系數、中介數、接近度、信息度、子圖和特征向量等，在具體分析中我們也稱這些測度為參數.

定義1 節點的度中心性（Degree Centrality，DC）[11]：某個節點的度定義為與節點i相連的節點個數，記為d（u）.在代謝網絡中表示與這個酶作用的酶的數量.從直觀上看，節點度越大意味著該節點在某種意義上越重要，在節點關系上具有某種優勢.

定義2 聚集系數中心性（Clustering Coefficient，CU）[12]：某個節點的聚集系數刻畫了這個節點周圍的節點彼此之間聯系的緊密程度，整個網絡的聚集系數CU就是所有節點聚集系數CUu的平均值.度為d（u）的節點u，其d（u）個鄰居之間實際存在的邊數Eu與它們之間最多可能邊數d（u）（d（u）－1）／2之比即為u的聚集系數，即：

定義3 中介系數中心性（Betweenness Centrality，BC）[13]刻畫節點在節點對之間通信的影響程度，被定義為網絡所有的最短路徑中經過當前節點的數目.它反映節點在特定網絡拓撲結構中所處位置的樞紐程度，通過它可以有效地區別出模塊內、外的邊.如果ρ（i，j）表示從節點i到節點j最短路的條數，而ρ（i，j，k）是在這些最短路徑中經過節點k的數量，因此，k節點的中介數為：

定義4 接近度中心性（Closeness Centrality，CC）[14]：如果節點與網絡中所有其它節點的距離都很短（距離之和最小），那么它的接近度就高，因此該測度在一定程度上能夠反映出節點在全局或整體上的中心程度.設μuv是節點u與節點v之間的距離，那么節點u的接近度為：

定義5 信息中心性（Information Centrality，IC）[15]：Newman認為信息中心性實際上是另一種接近度，本質上是測量以節點u為端點的路徑的調和平均長度，如果u通過許多短的路徑與其他節點相連則意味著其平均路徑長度小而信息度（接近度）大.如果A是一個網絡的鄰接矩陣，D是各節點的度的對角矩陣，J矩陣是所有的元素均為1的矩陣.那么IC可定義為通過倒轉矩陣B，而B＝D－A＋J，從而得到矩陣C＝（cij）＝B－1，信息矩陣Iij＝（cii＋cjj－cij）－1，那么第i個節點的信息中心值可通過下式的調和平均值來定義：

定義6 子圖中心性（Subgraph Centrality，SC）[16]：表示網絡節點出現在不同連通子網中的次數，即：

其中μl（i）是以節點i開始和結束的步長數，即：從i節點開始的經過L步再回到起點的閉路數.（v1，v2，…，vn）是RN的一個正交基 ，RN是由A矩 陣的特征向量組成的矩陣，這些特征向量分別與相關的特征值λ1，λ2，…，λn對應.vj（i）是Vj的第i個元素.

定義7 特征向量中心性（Eigenvector Centrality，EC）[17]定義為網絡的鄰接矩陣A的主特征向量.它模擬了一種各個節點同時影響其鄰居節點的程度的機制.這種特征向量的決定公式為λe＝Ae，其中A是圖的鄰接矩陣，λ是特征值，e是特征向量.因此，節點i的EC定義為特征向量e1的第i個元素e1（i），而e1對應于A矩陣的最大的特征值λ1，如下列公式（6）所示：

如果一個節點具有較高的特征向量評分，那么這個節點被認為是關鍵的，這意味著與其毗鄰的其他節點也具有較高的得分數.

2 主成分分析

將包含各方面拓撲信息的7個節點中心性進行主成分分析，得到一個網絡中心性參數.

2.1 選取主成分

主成分分析方法（PCA）[18]是將原來較多的指標簡化為少數幾個新的綜合指標的多元統計方法.PCA能夠用較少的幾個綜合指標盡量多地反映原來較多變量指標所反映的信息，同時它們之間又是彼此獨立正交的.這些互相正交的新變量是原先變量的線性組合，叫做主成分.

記X1，X2，…，Xn為原變量指標，F1，F2，…，Fp為新變量指標，n≥p.主成分分析實際上是做一個線性變換：

滿足如下條件：

每個主成分的系數平方和為1；且主成分之間相互獨立，即無重疊的信息；Var（F1）為方差函數，每個成分的方差依次遞減，重要性依次遞減.

對代謝網絡中n個節點的p個中心性進行主成分分析 PCL（si）＝ PCL（si1，si2，…，sip），i＝ 1，2，…，n，選取其主成分si0，具體的實現步驟如圖1所示.

圖1 選取主成分實現步驟Fig.1 The process of select principal component

設原始數據矩陣為X，n為代謝網絡中節點的個數，p為中心性屬性的個數，且

1）對矩陣X中的原始數據進行標準化處理.

將原始數據矩陣X進行標準化，得到標準化后的協方差矩陣；計算平均值和標準方差的公式為：

2）建立標準化后的p個指標相關系數矩陣R＝ （rij）；i，j＝1，2，…，p，其中rij的計算公式為：

3）求相關系數矩陣R的特征值和特征向量，并對特征值和特征向量按照特征值的大小進行排序.排序后特征值為λ1≥λ2≥…≥λn≥0.相應的標準化，正交特征向量為：

4）計算各主成分的方差貢獻率Qt和累積方差貢獻率qt.確定主成分個數，再求出各樣本的主成分值矩陣，并以每個主成分的方差貢獻率為權數計算各樣本的綜合評價值.其中，方差貢獻率qt表示第t個主成分提取原始n指標的信息量，累計方差貢獻率Qt表示前t個主成分保留的原始信息量.計算公式如下：

5）確定主成分個數有兩種方法：一種是使前t個主成分的累積方差貢獻率達到一定的要求，如可取Qt≥80%，這是實際中常用的一種方法.第二種是采用平均數原則，選用特征值大于1的所有指標.由于第一主成分在主成分中占有很大的比重，可以表示出大部分的特征信息.為易于計算比較，本章采用第一種方法，選取第一主成分為研究對象.

6）確定第一主成分中心性數值si0（i＝1，2，…，n），以累積方差貢獻率q（i）（i＝1，2，…，p）為因子作上述種中心性的線性組合得到第一主成分中心性數值：

2.2 主成分數據序列化

主成分分析得到的數據是無序的，為了分析方便，需要將主成分數據序列化.首先搜索數據的最大值M和最小值m，設梯度等級為N（N的取值大小可由網絡中節點數目多少而定，但必須是2的冪次，易于小波計算），將區間[m，M]N等分，每個區間的長度為h＝ （M－m）／N.然后判斷每個主成分數列的元素的值，將每個區間包含的主成分元素的個數記為zi，Z＝ ｛zi＝Count[（m＋i）h≤si≤m＋（i＋1）h]，i＝0，1，2，…，N－1｝，其中Count為記數函數.

圖2是對物種hsa的代謝網絡處理得到的結果.先用主成分法提取代謝網絡中所有節點中心性的第一主成分，并對其數值進行統計，將數值從最小值到最大值劃分128等分，統計每一等分內節點的個數并歸一化，即可得到此圖，從而完成數據有序化.

圖2 物種hsa代謝網絡的主成分數據有序化Fig.2 The principal components sequence of hsa

有序化后的主成分數據分布圖可以作為一種特征曲線來反映代謝網絡內在的某些特征.本文的主要工作就是將代謝網絡的這種特征曲線提取出來，通過比較不同代謝網絡的特征曲線之間的相似度，總結特點，發現規律.

3 小波分析與相似性比較

3.1 小波分析

小波分析[19]是運用傅里葉（Fourier）變換的局部化思想，進行時空序列分析的一種數學方法，是在L2（R）空間內利用小波基函數對數學表達式展開與逼近.它是一種函數的快速高效、高精度的近似方法.從信號處理的角度來看，作為一種新的時頻分析工具，小波克服了Fourier分析方法無法反映時間域上局部信息的缺陷.對于代謝網絡信息的局部性質描述和相似性比較具有重要的意義.

代謝網絡拓撲屬性的主成分序列化數據是一種典型的非平穩信號，利用小波實施時頻分析時，由于同時具有時間和頻率的局部特性以及多分辨分析特性，使得對非平穩信號的特征處理變得相對容易.

應用小波的Mallat算法對主成分序列化后的主成分信號進行小波分解.令信號T＝Hif，則其為能量有限信號f∈L2（R）在分辨率2j下的近似，Hif可以通過低通濾波器分解為f在分辨率2j－1下的近似Hi－1f，通過高通濾波器Hif分解得到分辨率2j－1與2j之間的細節Di－1f，其分解過程如圖3所示.

圖3 主成分序列化信號小波頻率分解示意圖Fig.3 The wavelets decomposition for aprincipal components sequence

圖4反映了小波強大的濾波功能，它可以將主成分序列化后的原始信號（如圖2）分解為低頻信號和高頻信號兩部分.低頻信號主要反映了原始信號中的主要部分，如圖4（a）所示，而高頻部分則是對局部細節特征的提取，如圖4（b）所示.本文主要利用小波變換提取低頻信號，研究曲線的整體趨勢.

圖4 物種hsa代謝網絡小波頻率分解實例圖Fig.4 The signal wavelets transformation of hsa

3.2 主成分特征相似性

因為主成分數據序列化后的曲線可以作為一種特征曲線來反映代謝網絡拓撲上的內在特征，所以我們提取不同物種的小波低頻信號，比較它們的相似程度，即可得到兩物種的相似度，相似度定義如下.

定義8 設經過J次小波低通濾波的主成分序列化信號為W＝ （w0，w1，…，wN－1），為了比較兩個信號W＝（w0，w1，…，wN－1）和W＊＝（w＊0，w＊1，…，w＊N－1）的近似程度，兩小波信號的相似度定義為：

其中wk為對應在第k個區間內節點的個數.由定義式可知rW，W＊∈ [－1，1]，當rW，W＊＝1時信號W＝W＊，此時信號的相似程度達到最高.

4 基于主成分分析和小波變換的代謝網絡比較方法MWD

由于主成分分析法能夠降維并進行多參數分析，因此引入了主成分分析法來處理7個中性能參數，從而達到綜合描述網絡的中心性的目的.另一方面，小波分析能研究曲線的變化趨勢，因此運用小波低頻系數來研究序列化的主成分趨勢特性，從而為研究生物進化提供理論和數學基礎.

基于主成分分析和小波分析的生物代謝網絡比較方法MWD，首先在計算得到7維（7種中心性參數）網絡節點中心性基礎上，用主成分分析法（PCA）對7維數據提取主成分處理得數據A0，并基于統計方法構造主成分有序化數列.然后運用小波得到特征信號W，最后由公式計算相似度rW，W＊.通過rW，W＊可以比較兩個物種，推算出兩物種代謝網絡的距離，并推斷出它們之間的進化關系.

方法MWD的具體步驟如下：

步驟1 由數據庫中代謝網絡的數據，利用7個中心性參數的定義，對網絡中的第i個節點用7維數組si＝ （si1，si2，…，si7），i＝1，2，…，n進行描述；

步驟2 運用主成分分析方法提取多個中心性參數 的 主 成 分，PCL（si）＝ PCL（si1，si2，…，si7）＝（si0）；

步驟3 對主成分分量｛si0，i＝1，2，…，n｝序列化，取N＝128進行序列化得到Z＝（z1，z2，…，z128）；

步驟4 對序列化數據Z＝ （z1，z2，…，z128）進行7次小波分解（這里的小波分解次數為7，只是一個經驗參數，并不是不可改變的，它是經過多次實驗比較實驗效果而確定的，經過7次小波分解的信號所得出的低頻信號平滑程度較好，較接近真實的信號）即令T＝H7f，則經過7次小波低通濾波，序列化數據Z＝ （z1，z2，…，z128）變為特征信號（記為W數列）W＝ （w1，w2，… ，w128）；

步驟5 對不同的網絡分別得到特征信號W和W＊，并由式（13）計算出不同網絡的相似度rW，W＊，對網絡比較進行量化.

5 實 驗

5.1 相似性比較結果及分析

實驗數據來自于周婷婷等[2]基于KEGG LIGAND數據庫重建的代謝網絡.本章對此數據集中的109個不同進化階段的代謝網絡進行了分析.為了探討物種的代謝網絡是否具有結構特異性，將三個進化階段物種的代謝網絡分別與三組隨機網絡進行比較；為了證明代謝網絡的進化是有一定規律的，將同一進化階段不同物種代謝網絡間進行比較、不同進化階段的不同物種代謝網絡間進行比較.

圖5為在109個物種代謝網絡的W數列三維視圖，一個曲線圖代表一個物種的W數列.此圖直觀展示了通過基于主成分分析和小波變換方法提取出的網絡特征.該圖表明，代謝網絡在整體而言存在著某種共性（曲線表現出的趨勢近似），而在局部又存在著一些差異，進而運用特征信號相似度的概念定量的闡述這種共性與差異，即根據W數列求它們的相似度.

圖5 109個物種代謝網絡的W數列三維視圖Fig.5 WSeries of 109species metabolic networks

通過對109個物種的代謝網絡與隨機網絡的相似度、109個物種的代謝網絡在生物的三個進化階段內外的相似度進行計算比較，可以得到以下結果：代謝網絡與隨機網絡的相似度的平均值為0.359 8，方差為0.037 5；任意兩個不同物種代謝網絡之間相似度的平均值為0.693 5，方差為0.034 8.不同的物種代謝網絡之間存在較高的相似性，而與隨機網絡的相似度比較低；方差值都很小，相似度都分布在平均值附近，異常值較少.說明生物代謝網絡自身結構具有特異性.

同一進化階段內的不同物種代謝網絡之間相似度的平均值為0.884 5，方差為0.029 3，存在很高的相似度，說明生物代謝網絡在短時間進化中具有一定的穩定性；不同進化階段中的不同物種代謝網絡之間相似度的平均值為0.583 3，方差為0.024 6，相似度相對較低，說明物種代謝網絡在長時間進化中發生了有規律的改變.

分別從古細菌、細菌、真核生物三個進化階段中各隨機選取4個物種的代謝網絡數據進行舉例說明.表1為12個物種代謝網絡與3組隨機網絡的相似度比較，表2為12個物種代謝網絡之間的相似度比較.

由比較結果分析可以證明，本章提出的方法可以將代謝網絡與隨機網絡區分開；還可以用于判斷某一物種大致所屬的進化階段或者某一物種與哪個階段的模式生物在進化過程中比較接近，對代謝網絡的進化研究起到一定的輔助作用，方便生物學家做進一步研究.

5.2 與集合論方法比較分析

在基于整體網絡比較策略中，集合論的方法是用Jaccard距離來確定生物網絡的進化距離[8].本章的MWD方法可以得到物種間的相似度rW，W＊.物種的相似度和物種間的距離，本身的度量單位不同.實際中的距離表示兩物體間相距的遠近，這里物種間的距離概念指兩物種間相似程度的區別，則相似程度小則距離大，相似程度大的則距離小.故可根據需要，將距離與相似度進行相互轉換.將rW，W＊轉換成JB的計算公式如公式（14）所示.

定義9 物種間的距離JB指在兩物種相似程度上的區別.根據定義8，相似度rW，W＊本質上是兩個向量的內積（余弦值），而物種間的距離JB本質上是向量之間的夾角，也就是兩物種分離開的程度即角度越大，距離越大，反之角度越小距離越小.將rW，W＊轉換成JB的計算公式如下：

表1 12個物種代謝網絡與3組隨機網絡的相似度比較Tab.1 Similarity between 12species metabolic networks and 3random networks

表2 12個物種三個進化階段的物種代謝網絡之間的相似度比較Tab.2 Similarity of 12species metabolic networks in three evolutionary stages

以公認的RNA距離[20]作為參考值，大腸桿菌E.coli為模式生物，將MWD方法計算出來的73個物種和模式生物的距離JB與用集合論的方法計算出來的Jaccard距離進行比較.

首先選取文獻[20]中涉及的73個物種（共62種細菌和11種古細菌），根據其與大腸桿菌的RNA距離，將其有序排列；其次，用基于集合論方法計算73個物種與E.coli系統層酶圖之間的進化距離即Jaccard距離，記為JC；用基于本章的MWD方法，計算73個物種與E.coli系統層酶圖之間的相似度，轉換為物種間的距離，記為JB；最后，以RNA距離為參考，將JB與JC進行比較，如表3和圖6所示.

表3 與RNA距離誤差的均值和標準差Tab.3 The mean and standard deviation of the error with RNA－distanse

表3列出了JC和JB與RNA距離之間誤差的均值與標準差，以數值的形式對比較結果進行顯示.與JC相比，JB具有更小的平均誤差，直觀地說明了本方法推算出進化距離的合理性.

圖6 以RNA距離為參考，JC與JB的比較結果圖Fig.6 The comparison of JC，JB and RNA－distance

分析JB與JC根據RNA距離的物種排序畫出的曲線比較結果，如圖6所示（x軸軸坐標分別對應73個物種按RNA距離的順序排列，y軸對應的物種間的距離）.圖中帶方塊實線表示作為參考的RNA距離，帶圈實線表示JC，帶星實線表示JB.從直觀上顯示，本方法推算出的進化距離比Jaccard距離更接近參考值.但是，曲線尾部雖然呈上升趨勢，其值比RNA小很多，反應出一些問題：

MWD方法本身進行了主成分分析、小波處理，提取出的特征信號，最后僅得出一個數值，大大壓縮了網絡信息量，數據可能失真.

小波低頻濾波保留了大部分的平穩信號，波動大的奇異點被濾掉了，而這些點往往是區別不相似的地方.所以方法本身對相似度大的測得準，相似度小的測得不準，所以距離大的波動大，曲線尾部趨勢不明顯.

定義9中JB的定義方法只是簡單直觀地從余弦求夾角，是否還可以找到更合理的定義方法還有待進一步研究.

生物數據本身不足，一個物種普遍只有幾百個節點，數據基礎不足.

當然不排除作為參考的RNA距離本身也存在著誤差.

綜上所述，用MWD算法來比較兩物種的網絡，可以很好地表現出網絡之間的相似度.但是用來測量兩物種間的距離，相似度高的可以得到很好的效果，不過相似度低的物種之間的距離不能很好地表現，還有待進一步研究.

6 結 論

本章提出的基于主成分分析和小波變換的代謝網絡比較方法，首次結合基于主成分分析方法分析代謝網絡的拓撲結構屬性，將小波分析運用于代謝網絡的進化距離分析中.該方法能有效地分析不同物種代謝網絡結構的相似程度，以揭示代謝網絡的物種特異性，為代謝網絡的進化研究提供數學基礎；并且能夠通過比較其相似性，較有效地區分不同進化階段物種的代謝網絡，可以對不同進化階段的代謝網絡的區分提供參考作用，探索進化規律；運用該方法將某物種的代謝網絡與模式物種代謝網絡的距離進行比較，發現該方法適用于推算相似度較高的網絡距離，而對相似度較低的網絡距離效果不明顯.

[1] MA H，ZENG A P.Reconstruction of metabolic networks from genome data and analysis of their global structure for various organisms[J].Bioinformatics，2003，19（2）：270－277.

[2] ZHOU Ting－ting，YUNG Kin－fung.MetaGen：apromising tool for modeling metabolic networks from KEGG [J].Progress in Biochemistry and Biophysics，2010，37（1）：63－68.

[3] FORST C V，SCHULTEN K.Phylogenetic analysis of metabolic pathways[J].Journal of Molecular Evolution，2001，52（6）：471－489.

[4] PINTER R Y，ROKHLENKO O，YEGER－LOTEM E，etal.A－lignment of metabolic pathways[J].Bioinformatics，2005，21（16）：3401－3408.

[5] Oh S J，JOUNG J G，CHANG J H，etal.Construction of phylogenetic trees by kernel－based comparative analysis of metabolic networks[J].BMC Bioinformatics，2006，7：284.

[6] AGUILAR D，AVILES F X，QUEROL E，etal.Analysis of phenetic trees based on metabolic capabilites across the three domains of life[J].Journal of Molecular Biology，2004，340（3）：491－512.

[7] MA H W，ZENG A P.Phylogenetic comparison of metabolic capacities of organisms at genome level[J].Molecular Phylogenetics and E－volution，2004，31（1）：204－213.

[8] FORST C V，FLAMM C，HOFACKER I L，etal.Algebraic comparison of metabolic networks，phylogenetic inference，and metabolic innovation[J].BMC Bioinformatics，2006，7（1）：67.

[9] TOHSATO Y.A method for species comparison of metabolic networks using reaction profile[J].IPSJ Digital Courier，2006，2：685－690.

[10] ZHOU Ting－ting，KEITH C C Chan，PAN Yi，etal.An approach for determining evolutionary distance in network－based phylogenetic analysis[J].Bioinformatics Research and Applications Lecture Notes in Computer Science，2008，4983：38－49.

[11] JEONG H，MASON S，BARABáSI A L，etal.Lethality and centrality in protein networks[J].Nature，2001，411（6833）：41－42.

[12] WATTS D J，STROGATZ S H.Collective dynamics of“smallworld”networks[J].Nature，1998，393：440－442.

[13] JOY M，BROCK A，INGBER D E，etal.High－betweenness proteins in the yeast protein interaction network[J].Journal of Biomedicine and Biotechnology，2005，2005（2）：96－103.

[14] WUCHTY S，STADLER P F.Centers of complex networks[J].Journal of Theoretical Biology，2003，223（1）：45－53.

[15] STEVENSON K，ZELEN M.Rethinking centrality：methods and examples[J].Social Networks，1989，11（1）：1－37.

[16] ESTRADA E，RODRíGUEZ－VELáZQUEZ J A.Subgraph centrality in complex networks[J].Phys Rev E，2005，71（5），056103.

[17] BONACICH P F.Power and centrality：a family of measures[J].A－merican Journal of Sociology，1987，92（5）：1170－1182.

[18] 陳峰.主成分回歸分析[J].中國衛生統計，1991，8（1）：20－22.CHEN Feng.Principal component regression analysis[J].Chinese Journal of Health Statistics，1991，8（1）：20－22.（In Chinese）

[19] DAUBCHIES I.Orthonormal bases of compactly supported wavelets[J].Communications on Pure and Applied Mathematics，1988，41（7）：909－996.

[20] ZHANG Y，ZHANG Z，LING L，etal.Conservation analysis of small RNA genes in escherichia coli[J].Bioinformatics，2004，20（5）：599－603.