基因共表達網(wǎng)絡的構建及分析方法研究綜述

汪濤等

摘 要：隨著高通量生物實驗技術的快速發(fā)展，特別是基因芯片和新一代測序技術的發(fā)展，全基因組范圍內(nèi)的基因表達數(shù)據(jù)呈爆炸式增長。利用網(wǎng)絡生物學的方法對高通量基因表達數(shù)據(jù)進行分析和挖掘已經(jīng)成為生物信息學重要的研究方向。對基因共表達網(wǎng)絡的研究與分析從系統(tǒng)層面上加深了研究人員對生物系統(tǒng)的認識。本文綜述了基因共表達網(wǎng)絡的構建和分析的常用方法，主要包括基因相似性度量方法、閾值選擇方法、拓撲分析方法、基因模塊識別及其功能注釋注釋方法，并對一些常用的分析工具進行了分析總結。

關鍵詞：基因共表達網(wǎng)絡；基因模塊；功能和拓撲分析

中圖分類號：TP391 文獻標識號：A 文章編號：2095-2163（2014）06-

Abstract： With the rapid development of high-throughput biological experiments technology， particularly the development of gene microarray and next generation sequencing technology， the genome wide gene expression data grow explosively. Network biology has become an important research direction in bioinformatics， which enhances the analysis of high throughput gene expression data. The study and analysis of gene co-expression network help the researchers to understand the biological system in system level better. This article reviewes the common methods on constructing and analyzing gene co-expression network， including the methods of measuring similarity of genes， selecting proper threshold， analyzing topological structure， detecting gene module and functional annotation on gene module. In the end， the paper summarizes some commonly used tools for analyzing gene co-expression network.

Keywords： Gene Co-expression Network； Gene Module； Functional and Topological Analysis

0 引 言

后基因組時代生命科學研究的一個主要目的是理清生物細胞內(nèi)所有分子以及分子之間的聯(lián)系，并且揭示分子之間相互作用以及如何決定細胞生命功能的內(nèi)在機理[1]。隨著系統(tǒng)生物學和復雜圖理論的發(fā)展，分子生物網(wǎng)絡的研究為探索復雜生命活動提供了有力工具。分子生物網(wǎng)絡在系統(tǒng)層面反映了生物分子的相互作用關系，因而在相當過程中有助于研究者深入理解生物細胞中各種生物分子是如何相互作用、進而行使生物功能的完整處理實現(xiàn)過程。目前人們已經(jīng)對各種類型的分子生物網(wǎng)絡進行了廣泛的研究，如基因共表達網(wǎng)絡（gene co-expression network）、基因調(diào)控網(wǎng)絡（gene regulatory network）、蛋白質(zhì)相互作用網(wǎng)絡（protein-protein interaction network）、代謝網(wǎng)絡（metabolic network）等。這些探索已經(jīng)從對單個生物分子研究層面上升到解析研究分子相互作用關系的系統(tǒng)研究層面，并且產(chǎn)生了豐碩的研究成果。

基因芯片技術以及新一代測序技術的應用，使得全基因組范圍內(nèi)基因表達數(shù)據(jù)得以快速累積。僅僅對單個基因功能水平的研究已經(jīng)限制了人們探索生物細胞行使生命功能的視野和進程。利用系統(tǒng)生物學的方法構建基因共表達網(wǎng)絡，從而由系統(tǒng)層面揭示基因之間的相互關系已經(jīng)成為一個主要的研究方向。本文中，主要對基因共表達網(wǎng)絡的構建方法以及常用的分析方法進行了關注及論述。

1基因共表達網(wǎng)絡的構建

基因共表達網(wǎng)絡大多是以基因間表達譜數(shù)據(jù)的相關性為基礎而實現(xiàn)構建的。在基因共表達網(wǎng)絡的表示中，經(jīng)常使用圖模型來描述基因之間的關系。圖中的節(jié)點代表基因，邊表示兩個基因之間的共表達相互作用關系?；蚬脖磉_網(wǎng)絡的構建主要分為兩個步驟，第一是對所有基因進行相似性度量；第二是通過閾值的選擇確定共表達網(wǎng)絡的邊。以下則對其展開具體分析。

1.1數(shù)據(jù)來源及表示

在分子生物學領域，基因表達譜是指基因表達活性的有效度量。從基因表達譜的數(shù)據(jù)來源上進行分類，常用于基因共表達網(wǎng)絡構建的表達譜數(shù)據(jù)主要將劃定為兩類：一類是基因芯片（microarray）數(shù)據(jù)，另一類是RNA-seq數(shù)據(jù)?；虻谋磉_譜數(shù)據(jù)可以用一個n*m的矩陣X=[xij]來表示，數(shù)學表述如公式（1）所示。其中，第i行數(shù)據(jù)xi （i = 1，…，n）對應一個基因的表達譜，矩陣中的列則能夠反映在不同樣本或者時間點下該基因的表達水平。

因相似性度量方法

在基因共表達網(wǎng)絡的構建和分析中，經(jīng)常需要對兩個基因進行表達相似性度量?；蜷g的相似性有多種度量方式，若從計算方法上進行分類，主要可以分為基于表達譜的相似性度量和基于拓撲結構的相似性度量。具體地，基于基因表達譜的相似性度量指標主要通過計算不同基因表達譜的線性或非線性相關系數(shù)而最終得到。其中，常見的線性相關性指標主要有皮爾森相關系數(shù)（Pearson Correlation Coefficient， Pcc）、斯皮爾曼相關系數(shù)（Spearman Correlation Coefficient），以及排除其他變量影響的偏相關系數(shù)（Partial Correlation Coefficient）等。而非線性相關性指標則可列舉實例，諸如互信息（Mutual Information， MI）等。在基因模塊識別的過程中，有些研究者又利用網(wǎng)絡拓撲結構對基因相似性實現(xiàn)了進一步的度量。

1.3閾值選取方法

設定合理的閾值，將具有潛在生物功能的邊保留下來，是基因共表達網(wǎng)絡構建的重點及關鍵。以計算方法為依據(jù)，大致可以分為以下幾類：基于人工設定的閾值選取方法、基于統(tǒng)計顯著水平的閾值選取方法、基于對照試驗的閾值選取方法、基于網(wǎng)絡拓撲結構的閾值選取方法以及基于多種方法綜合的閾值選取方法等。在此，即對各類方法做以詳細的探討與分析。

首先，基于人工設定的閾值選取方法是最簡單、也是應用最廣的方法。其中有兩類模型使用最為廣泛：

（1）選取固定的閾值t，將相似性大于閾值t的基因?qū)ΡＡ粝聛?，例如Tom C Freeman[4]等人的文章中就采取了這種策略。

（2）將所有基因?qū)Π凑障嗨菩韵禂?shù)進行排序，選擇前百分之x（例如前1%）的基因?qū)M行保留，Ala等人[5]選取前0.5%作為閾值，構建保守的共表達網(wǎng)絡。但這種方法最大的缺點在于其具有的武斷性，而并未考慮到生物網(wǎng)絡本身的特性，也沒有驗證其生物統(tǒng)計顯著水平。對應地，基于統(tǒng)計顯著水平的閾值選取方法往往通過與隨機網(wǎng)絡進行比對，得到不同閾值下的顯著性水平，由此而進行閾值選擇。并且，顯著性水平往往多采用p值（p-value）、q值（q-value）等來衡量。然而這種方法卻無法保留那些相似性系數(shù)較低、但卻具有生物意義的邊。

其次，基于對照試驗的閾值選取方法解決了計算方法帶來的隨機性。這種方法在基因芯片的實驗階段，通過在基因芯片中引入與目標物種親緣關系較遠的物種基因而形成對照組。理論上，對照組的基因由于未與任何RNA雜交而不曾產(chǎn)生任何信號，但是現(xiàn)實中隨機雜交的情況卻使得對照組將產(chǎn)生一定的噪音信號。而且，這種噪聲信號可以作為參考，用于輔助目標網(wǎng)絡的閾值選擇。Voy等人[6]證明了利用這種方法保留下的邊具有顯著的生物統(tǒng)計意義。但是，該種方法卻增加了實驗成本，并且結果將十分依賴所選擇的對照樣本。

再有，基于網(wǎng)絡拓撲結構的閾值選取方法考慮到了基因共表達網(wǎng)絡所具有的特性，如無標度（scale-free）以及小世界（small world），并通過優(yōu)化閾值的方法，使得網(wǎng)絡的某些拓撲特性得以明顯體現(xiàn)。Horvath等人[7]即提出了一種“軟”閾值選取方法，利用網(wǎng)絡達到無標度的拓撲結構而以此來確定閾值。Elo等人在文獻[8]中則提出了一種基于聚集系數(shù)（Clustering Coefficient）的閾值選擇策略。由于基因共表達網(wǎng)絡具有小世界的特性，因此往往比隨機網(wǎng)絡具有更高的聚集系數(shù)。作者利用這種特點，將閾值選擇的問題轉(zhuǎn)化為具體的關于聚集系數(shù)C*的優(yōu)化問題。Borate等人又在文章[9]中將基于最大團和基于譜圖理論的閾值選擇方法等進行了對比研究。當閾值從高到低下降時，網(wǎng)絡中的最大團的數(shù)目是指數(shù)增長的。根據(jù)這一原理，基于最大團的方法在最大團的數(shù)目隨閾值變化的趨勢線中將選擇一個關鍵拐點作為閾值。該選取拐點往往是最大團的數(shù)目突然增加2倍或者3倍時的閾值點。而與此類似，利用譜圖理論的閾值選擇方是基于網(wǎng)絡的拉普拉斯矩陣的特征值和特征向量，來發(fā)現(xiàn)網(wǎng)絡的基因模塊，并會選擇一個模塊劃分最優(yōu)時的閾值。

最后，基于多種方法綜合的閾值選取方法通常是融合了上述方法的優(yōu)點，從而用多種指標來選取閾值。例如，Langston等人[10]利用了本體距離、統(tǒng)計顯著性以及多種圖的拓撲特性來進行閾值選擇。

2基因共表達網(wǎng)絡的拓撲分析方法

生物網(wǎng)絡的快速發(fā)展表明分子網(wǎng)絡遵從著某些普適性的規(guī)則[1]。這些規(guī)則通常會在網(wǎng)絡的拓撲特性中得到相應體現(xiàn)。因此，對有關基因共表達網(wǎng)絡的拓撲展開分析即是從系統(tǒng)層面了解基因共表達關系的重要工具。根據(jù)度量拓撲結構中的基因個數(shù)的不同，可以將拓撲分析方法大致分為兩類，也就是：全局網(wǎng)絡拓撲分析和網(wǎng)絡中心性分析。

在研究進展中，為了從全局角度刻畫網(wǎng)絡的拓撲結構，一些學者即利用圖論的的方法定義了許多衡量網(wǎng)絡全局性拓撲結構的度量指標，常見的主要包括平均度（average degree）和度分布（degree distribution）、聚集系數(shù)（clustering coefficient）、平均路長（average path length）、直徑（diameter）等。這些全局性的度量指標能夠很好地反應出基因共表達網(wǎng)絡的三種全局拓撲特性，分別是：無標度分布（scale-free distribution）、小世界特性（small world property）、功能模塊網(wǎng)絡（functional modular network）。另外，基因共表達網(wǎng)絡中不同節(jié)點在網(wǎng)絡中的重要性是不同的，通常用中心性（Centrality）來度量網(wǎng)絡節(jié)點的地位。中心性往往體現(xiàn)在網(wǎng)絡的拓撲結構上，度很高的節(jié)點或者是起著關鍵連接作用的節(jié)點都可能在某些生物途徑中發(fā)揮著重要的作用，因此人們往往根據(jù)所研究問題的不同來定義節(jié)點的中心性。具體地，較為常用的中心性度量指標主要有度中心性（degree centrality）、親近中心性（closeness centrality）、居間中心性（betweenness centrality）、特征向量中心性（eigenvector centrality）等。

CentiBiN[11]就是一款專門用于生物網(wǎng)絡中心性的計算和可視化軟件，其中集成了17種無向圖的中心性度量方法（針對有向圖則匯總了15種中心性度量方法），并且可以計算網(wǎng)絡直徑、平均路長等全局性網(wǎng)絡拓撲指標。另外，諸如Cytoscape[12]、Pajek[13]、Visone[14]、VisANT[15]等復雜網(wǎng)絡分析和可視化軟件工具，雖然不是專門為了計算網(wǎng)絡拓撲特性而產(chǎn)生的，但對于常用的網(wǎng)絡拓撲特性分析卻都能提供良好的支持。

3基因共表達網(wǎng)絡的模塊分析

3.1 基因模塊識別

基因功能模塊識別是基因共表達網(wǎng)絡分析中的最重要方法之一。研究證明，在基因共表達網(wǎng)絡中致密的連通子圖往往具有特定的生物學功能?；蚬脖磉_網(wǎng)絡模塊識別算法大多來源于對蛋白質(zhì)相互作用網(wǎng)絡的研究中。Junzhong Ji等人[16]對蛋白質(zhì)網(wǎng)絡中的模塊識別方法即做了詳細的綜述。在本文中，僅對目前常用于基因模塊識別的四種典型算法進行了深入分析，如基于層次的模塊識別方法（WGCNA[7]）、基于密度的模塊識別方法（MCODE[17]）、基于流模擬的模塊識別方法（MCL[18]）以及基于劃分的模塊識別方法（Qcut[19]）。下面依然對其展開綜合性論述。

首先，WGCNA[7]（Weighted Gene Co-expression Network Analysis）是以基于相關系數(shù)構建基因共表達網(wǎng)絡的代表性方法，其中集成了多種網(wǎng)絡分析方法，主要具有網(wǎng)絡構建、功能模塊探測、基因選擇、拓撲特性計算、數(shù)據(jù)模擬、可視化以及與其他軟件交互等功能。WGCNA使用了基于拓撲結構交疊的度量指標，利用無先驗知識的聚類策略來探測基因的功能模塊，實際默認使用的則是層次聚類方法——hclust[20]。層次聚類得到的系統(tǒng)樹圖的分支對應著可能的功能模塊，并且可以使用多種剪枝的方法如固定高度的剪枝方法或者兩種動態(tài)剪枝方法[21]來確定功能模塊。作為多種生物網(wǎng)絡構建的代表性方法，WGCNA現(xiàn)已成功地應用到多種研究的開展和實現(xiàn)當中。

其次，MCODE[17]（Molecular Complex Detection）算法是一種基于圖論（或網(wǎng)絡密度）的網(wǎng)絡模塊發(fā)現(xiàn)算法。算法共分為三個步驟：網(wǎng)絡節(jié)點加權、模塊預測以及模塊優(yōu)化處理。分別來說，網(wǎng)絡節(jié)點加權是指根據(jù)節(jié)點所在的、最高的k-core的密度來為網(wǎng)絡中的所有節(jié)點賦予一個權值。在第二步模塊的預測中，首先選取一個具有最高權值的節(jié)點作為種子節(jié)點，而后依次向外擴增，再將那些權值在閾值VWP之上的周邊節(jié)點依次納入模塊中，直到?jīng)]有節(jié)點再可包含進入該模塊為止。上述操作后，就將選擇剩余節(jié)點中權值最高的作為種子節(jié)點進行同樣的操作。而在第三步的模塊優(yōu)化處理中，算法會將那些不包含2-core的模塊刪除。并且用戶可以指定是否進行“fluff”和“haircut”操作來處理模塊邊緣節(jié)點。AllegroMCODE[22]是一款基于MCODE算法的Cytoscape[12]插件，可通過GPU進行加速，多將用于高效地挖掘基因模塊。

3.2 模塊的功能分析方法及工具

為了發(fā)現(xiàn)基因模塊的功能，往往需要對基因模塊進行模塊富集分析（Modular enrichment analysis， MEA）[23]。模塊富集分析多數(shù)時候是通過計算基因本體注釋數(shù)據(jù)庫[24]（Gene Ontology， GO ）或通路注釋數(shù)據(jù)庫（如KEGG[25]）中的術語（term）在每個模塊上的富集程度，再根據(jù)假設檢驗中的p值（p-value）來衡量富集到模塊上的功能的顯著性。計算p值的經(jīng)典統(tǒng)計檢驗方法主要包括卡方檢驗（Chi-square test）、Fisher精確檢驗（Fisher's exact test）、累計超幾何檢驗（cumulative hypergeometric test）等。模塊富集分析不僅繼承了單一富集分析（singular enrichment analysis， SEA）的一些特點（如在預選感興趣的基因集后，再計算各基因與注釋數(shù)據(jù)庫中的術語的富集得分），而且模塊富集分析在計算p值時還考慮到了模塊基因-基因之間的關系以及對應的術語-術語關系。這種術語-術語間關系的著重關注將可提高功能富集的敏感性和準確度。常見的能夠基于GO進行模塊功能富集工具可概略給出為：GO：：TermFinder[26]、Ontologizer[27]、 topGO[28]、 GENECODIS[29]、 ADGO[30]、GoToolBox[31]、DAVID[32]等。

4 結束語

基因共表達網(wǎng)絡的分析方法已成為高通量生物數(shù)據(jù)分析強有力的工具，并已經(jīng)廣泛應用于生物學研究中。但是目前的構建和分析方法還存在著一些不足，例如基因共表達網(wǎng)絡的構建方法繁多，時下還缺乏可靠的評價系統(tǒng)；模塊功能富集方法對連接度較低的節(jié)點的識別能力還有待進一步提高等。因此總地來說，對于基因共表達網(wǎng)絡的構建和分析也依然是目前一項頗具挑戰(zhàn)性的研究工作。

參考文獻：

[1] BARAbASI A L，OLTVAI Z N. Network biology： understanding the cell's functional organization[J]. Nature Reviews Genetics， 2004，5（2）： 101-113.

[2] RAVASZ E， et al. Hierarchical organization of modularity in metabolic networks[J]. science， 2002， 297（5586）：1551-1555.

[3] YIP A M， HORVATH S. Gene network interconnectedness and the generalized topological overlap measure[J]. BMC bioinformatics， 2007，8（1）： 22.

[4] FREEMAN T C， et al. Construction， visualisation， and clustering of transcription networks from microarray expression data[J]. PLoS computational biology， 2007， 3（10）： e206.

[5] ALA U， et al. Prediction of human disease genes by human-mouse conserved coexpression analysis[J]. PLoS computational biology， 2008， 4（3）： e1000043.

[6] VOY B H， et al. Extracting gene networks for low-dose radiation using graph theoretical algorithms[J]. PLoS computational biology， 2006， 2（7）： e89.

[7] LANGFELDER P， HORVATH S. WGCNA： an R package for weighted correlation network analysis[J]. BMC bioinformatics， 2008， 9（1）： 559.

[8] ELO L L， et al. Systematic construction of gene coexpression networks with applications to human T helper cell differentiation process[J]. Bioinformatics， 2007， 23（16）： 2096-2103.

[9] BORATE B R. Comparative Analysis of Thresholding Algorithms for Microarray-derived Gene Correlation Matrices， 2008.

[10] LANGSTON M A， et al. Innovative computational methods for transcriptomic data analysis： A case study in the use of FPT for practical algorithm design and implementation[J]. The Computer Journal， 2008， 51（1）： 26-38.

[11] KOSCHUTZKI D. CentiBiN Version 1.4. 2， in， 2006，CentiBiN Version， 2006，1（2）： 2004-2006.

[12] SHANNON P， et al. Cytoscape： a software environment for integrated models of biomolecular interaction networks[J]. Genome research， 2003， 13（11）： 2498-2504.

[13] BATAGELJ V， MRVAR A. Pajek-program for large network analysis[J]. Connections， 1998， 21（2）： 47-57.

[14] BAUR M， et al. Visone Software for visual social network analysis. in Graph Drawing. Springer，2002.

[15] HU Z， et al. VisANT： data-integrating visual framework for biological networks and modules[J]. Nucleic acids research， 2005， 33（suppl 2）： W352-W357.

[16] JI J， et al. Survey： Functional module detection from protein-protein interaction networks. Knowledge and Data Engineering[J]. IEEE Transactions on， 2014， 26（2）： 261-277.

[17] BADER G D， HOGUE C W. An automated method for finding molecular complexes in large protein interaction networks[J]. BMC Bioinformatics， 2003， 4： 2.

[18] HWANG W， et al. A novel functional module detection algorithm for protein-protein interaction networks[J]. Algorithms for Molecular Biology， 2006， 1（24）.

[19] RUAN J， ZHANG W. Identifying network communities with a high resolution[J]. Physical Review E， 2008， 77（1）： 016104.

[20] Kaufman L， Rousseeuw P J. Finding groups in data： an introduction to cluster analysis. John Wiley & Sons，2009，34.

[21] LANGFELDER P， ZHANG B， HORVATH S. Defining clusters from a hierarchical cluster tree： the Dynamic Tree Cut package for R[J]. Bioinformatics， 2008， 24（5）： 719-720.

[22] YOON J， JUNG W H. A GPU-accelerated bioinformatics application for large-scale protein interaction networks. APBC poster presentation， 2011.

[23] HUANG D W， SHERMAN B T， LEMPICKI R A. Bioinformatics enrichment tools： paths toward the comprehensive functional analysis of large gene lists[J]. Nucleic acids research， 2009， 37（1）： 1-13.

[24] GO T. Gene Ontology： tool for the unification of biology[J]. AMERICA N， editor. Nature Genetic， 2000， 25： 25-29.

[25] KANEHISA M， GOTO S. KEGG： kyoto encyclopedia of genes and genomes. Nucleic acids research， 2000， 28（1）： 27-30.

[26] BOYLE E I， et al. GO：： TermFinder—open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes[J]. Bioinformatics， 2004， 20（18）： 3710-3715.

[27] BAUER S， et al. Ontologizer 2.0—a multifunctional tool for GO term enrichment analysis and data exploration[J]. Bioinformatics， 2008， 24（14）： 1650-1651.

[28] Alexa A， Rahnenfuhrer J. topGO： enrichment analysis for gene ontology. R package version 2.8， 2010.

[29] CARMONA-SAEZ P， et al.， GENECODIS： a web-based tool for finding significant concurrent annotations in gene lists[J]. Genome biology， 2007， 8（1）： R3.

[30] NAM D， et al. ADGO： analysis of differentially expressed gene sets using composite GO annotation[J]. Bioinformatics， 2006， 22（18）： 2249-2253.

[31] MARTIN D， et al. GOToolBox： functional analysis of gene datasets based on Gene Ontology[J]. Genome biology， 2004， 5（12）： R101.

[32] ALVORD G， et al. The DAVID Gene Functional Classification Tool： a novel biological module-centric algorithm to functionally analyze large gene lists[J]. Genome Biol， 2007，8： R183.