基因共表達網絡的構建及分析方法研究綜述

汪濤等

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

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

中圖分類號：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 引 言

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

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

1基因共表達網絡的構建

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

1.1數據來源及表示

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

因相似性度量方法

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

1.3閾值選取方法

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

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

（1）選取固定的閾值t，將相似性大于閾值t的基因對保留下來，例如Tom C Freeman[4]等人的文章中就采取了這種策略。

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

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

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

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

2基因共表達網絡的拓撲分析方法

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

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

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

3基因共表達網絡的模塊分析

3.1 基因模塊識別

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

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

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

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

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

4 結束語

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

參考文獻：

[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.