轉(zhuǎn)自：WGCNA算法研究筆記 - 黎嫣 - 博客園 (cnblogs.com)
研究了近半年的算法永罚，記錄下來給自己一個交代重父，也應(yīng)該是考G前地最后一篇日志了屋谭。
Weighted Gene Co-Expression Network Analysis中文名有翻譯成加權(quán)關(guān)聯(lián)網(wǎng)絡(luò)分析的滔驶，感覺不是很恰當(dāng)，英文來得比較直接亭引。本來是佟昊從老汪那拿的一個課題嚼松，因?yàn)榭雌饋肀容^有意思就把文章找來慢慢啃虾啦，到現(xiàn)在算是搗鼓出點(diǎn)名堂了。方法是UCLA的一個教授提出來的锭吨，在文章中他將其歸類到系統(tǒng)生物學(xué)的研究方法中蠢莺，不過個人認(rèn)為由于其分析水平還是只停留在DNA芯片上，并未到達(dá)系統(tǒng)的程度零如，但是方法本身還是能夠預(yù)見一些incisive idea的躏将。我打算先介紹方法的基本思想，然后把聯(lián)系網(wǎng)上公布的sample data將該算法完整實(shí)現(xiàn)一遍考蕾，并解釋一些自己在看的時候遇到的關(guān)鍵問題祸憋。
Weighted Gene Co-Expression Network Analysis(以下簡稱WGCNA)，是一種從芯片數(shù)據(jù)中挖掘模塊(module)信息的算法肖卧。在該方法中module被定義為一組具有類似表達(dá)譜的基因蚯窥，如果某些基因在一個生理過程或不同組織中總是具有相類似的表達(dá)變化，那么我們有理由認(rèn)為這些基因在功能上是相關(guān)的塞帐，可以把他們定義為一個模塊(module)拦赠。這似乎有點(diǎn)類似于進(jìn)行聚類分析所得到結(jié)果，但不同的是葵姥，WGCNA的聚類準(zhǔn)則具有生物學(xué)意義荷鼠，而非常規(guī)的聚類方法(如利用數(shù)據(jù)間的幾何距離)，因此該方法所得出的結(jié)果具有更高的可信度榔幸。當(dāng)基因module被定義出來后允乐，我們可以利用這些結(jié)果做很多進(jìn)一步的工作，如關(guān)聯(lián)性狀(隨后會以這方面的應(yīng)用為Example)削咆，代謝通路建模牍疏，建立基因互作網(wǎng)絡(luò)，甚至進(jìn)行eQTL(這個確實(shí)很方便态辛，不過前提是實(shí)驗(yàn)題有錢去雜那么多芯片)麸澜。不過我個人從中獲益最多的是能加深人們對于生物體所選擇的這種Scale-Free Topology網(wǎng)絡(luò)調(diào)控的思考(下文會提到)。
WGCNA所分析的數(shù)據(jù)是芯片數(shù)據(jù)(當(dāng)然需要雜很多芯片奏黑，比如若要研究細(xì)胞凋亡炊邦，那么使用改方法需要實(shí)驗(yàn)者提供細(xì)胞凋亡各個時期的芯片數(shù)據(jù)编矾，以了解這一生理過程中細(xì)胞內(nèi)所有基因的表達(dá)變化)。
在co-expression network中馁害，每一個基因在一個特定時間或空間的表達(dá)情況被視做一個點(diǎn)(node)窄俏，可以簡單的理解成一張芯片上的一個基因的表達(dá)狀況就是網(wǎng)絡(luò)中的一個node。如果我們做了80張芯片碘菜，每張芯片上有8000個基因凹蜈，那么我們可以用一個808000的矩陣來表示實(shí)驗(yàn)結(jié)果。為了得到基因間的關(guān)聯(lián)情況忍啸，我們需要計(jì)算任何兩個基因間的相關(guān)系數(shù)(文章中采用Person Coefficient)仰坦，在經(jīng)過該步運(yùn)算以后，我們可以得到一個80008000的實(shí)對稱陣S计雌，sij表示第i個基因和第j個基因的Person Coefficient悄晃，即兩個基因的表達(dá)譜相似性。
下一步的分析是該方法的第一個靚點(diǎn)凿滤。為了知道兩個基因的表達(dá)譜是否具有相似性妈橄，需要人為規(guī)定一個閾值，只有當(dāng)基因間的Person Coefficient達(dá)到這一閾值后(如0.8)我們才認(rèn)為這兩個基因是相似的翁脆，否則則不相似眷蚓。為此人們定義了一個Adjacency Matrix，很顯然在通過以上步驟處理S矩陣后反番，得到的Adjacency Matrix將是一個0/1矩陣(該矩陣主對角線上元素被規(guī)定為0沙热，這一點(diǎn)很重要)。但是這種分析方法存在一個很明顯的局限恬口，即我們沒有理由認(rèn)為Person Coefficient為0.8的兩個基因與Coefficient為0.79的兩個基因是有顯著差別的校读，但是以上算法卻無法避免這一處境。而WGCNA采用了一種基于軟閾值的判定方法很好地避免了這一問題祖能。軟閾值的思想是通過權(quán)函數(shù)將Adjacency Matrix中的元素連續(xù)化(所以方法才稱之為Weighted Network)歉秫，常用的權(quán)函數(shù)包括sigmoid函數(shù)和power函數(shù)。
那么aij同樣是可分解的养铸，這一點(diǎn)很好證明雁芙。擁有這一性質(zhì)的優(yōu)勢在于可以簡化計(jì)算量，當(dāng)實(shí)對稱陣是可以分解的時候钞螟，我們只需要利用其分解后的一個向量就可以表示該矩陣兔甘，因此在實(shí)際應(yīng)用中能夠簡化計(jì)算機(jī)的運(yùn)算內(nèi)存使用。 如果直接拿Adjacency Matrix中的基因相關(guān)性指標(biāo)進(jìn)行module識別則未免有些過于簡單了鳞滨，為了保證芯片信息的充分利用洞焙，文章的作者提出了計(jì)算另一矩陣——topological overlap matrix(TOM)來衡量兩個基因的相關(guān)性。建立這一矩陣的思想在于，任何兩個基因的相關(guān)性不僅僅由他們的表達(dá)相似性直接決定澡匪，它還將A基因通過B基因與C基因的作用相關(guān)性納入AC基因的TOM矩陣值中熔任，以更精確地描述基因表達(dá)譜的類似性。
利用這種方式定義TOM矩陣中各元素的值是非常巧妙的唁情，它很好地滿足了我們所希望達(dá)到的目的疑苔。分子中對于lij
的定義，表示基因i通過任何基因與j關(guān)聯(lián)的adjacency值甸鸟，并將它們相加惦费，而aij
則表示基因i與基因j間的直接關(guān)聯(lián)性。分母的定義保證了wij
始終在0抢韭，1之間薪贫，我們可以考慮極端情況。當(dāng)Adjacency Matrix中除主對角線元素外所有元均為1時篮绰，
上式直接看不容易看懂后雷，把加和表達(dá)式展開之后則容易理解。
值得注意的是吠各，WGCNA方法只考慮了一階基因關(guān)聯(lián)，更高級的關(guān)聯(lián)可以用類似的方式去表示勉抓，但是并無此必要贾漏。首先芯片數(shù)據(jù)本身存在噪音，過度地提取信息未必會得到更好地結(jié)果藕筋，而且計(jì)算高階關(guān)聯(lián)會使算法的復(fù)雜度顯著增大纵散，即使高配置的服務(wù)器也未必能滿足計(jì)算要求。
為方便后面的module identification隐圾，還需定義一個dissimilarity matrix伍掀。根據(jù)前人的研究，方式如下：
其中dij
表示dissimilarity matrix中的元暇藏，這一方程只是用1減去wij
得到dissimilarity matrix蜜笤，在dij
上加指數(shù)的原因是empirical研究的結(jié)果，當(dāng)使用指數(shù)形式的dissimilarity進(jìn)行聚類分析會得到更distinct gene module盐碱。 得到了dissimilarity matrix后把兔，我們所需要進(jìn)行的工作就是聚類了，文章中使用的是hierarchical clustering方法瓮顽，各種聚類方法的優(yōu)劣不在本文的討論范圍中县好。使用聚類分析之后，各個module的identification也就完成了暖混。
了解了整個分析流程缕贡，我們有必要再進(jìn)一步了解一下分析細(xì)節(jié)。
首先是權(quán)函數(shù)的參數(shù)選擇，由于power function有一參數(shù)
這一參數(shù)的選擇勢必影響著module identification的結(jié)果晾咪。
為了選擇一個合適的參數(shù)值黔漂，我們有必要重新審視基因互作網(wǎng)絡(luò)的構(gòu)造。 網(wǎng)絡(luò)的數(shù)學(xué)名稱是圖禀酱，在圖論中對于每一個節(jié)點(diǎn)有一個重要概念炬守，即：度。一個點(diǎn)的度是指圖中該點(diǎn)所關(guān)聯(lián)的邊數(shù)剂跟。如下圖减途，如果不加以思考，人們很容易認(rèn)為生活中常見的網(wǎng)絡(luò)會是一種random network曹洽，即每一個節(jié)點(diǎn)的度相對平均鳍置。然而第二種圖，即scale-free network才是一種更穩(wěn)定的選擇送淆。Scale-free network具有這樣的特點(diǎn)税产，即存在少數(shù)節(jié)點(diǎn)具有明顯高于一般點(diǎn)的度，這些點(diǎn)被稱為hub偷崩。由少數(shù)hub與其它節(jié)點(diǎn)關(guān)聯(lián)辟拷，最終構(gòu)成整個網(wǎng)絡(luò)。這樣的網(wǎng)絡(luò)的節(jié)點(diǎn)度數(shù)與具有該度數(shù)的節(jié)點(diǎn)個數(shù)間服從power distribution阐斜。這為我們尋找最佳參數(shù)提供了理論依據(jù)衫冻。 這里做一點(diǎn)擴(kuò)展，我認(rèn)為時非常有必要的谒出。只要我們愿意抽象隅俘，Scale-free network大量存在于的生活中。人們的社交網(wǎng)絡(luò)笤喳、生物基因蛋白質(zhì)的相互作用为居、計(jì)算機(jī)網(wǎng)絡(luò)甚至sexually transmitted diseases均有這層關(guān)系。生物體選擇scale-free network而不是random network是有它進(jìn)化上的原因的杀狡，顯然對于scale-free network蒙畴，少數(shù)關(guān)鍵基因執(zhí)行著主要功能，這種網(wǎng)絡(luò)具有非常好的魯棒性捣卤，即只要保證hub的完整性忍抽，整個生命體系的基本活動在一定刺激影響下將不會受到太大影響，而random network若受到外界刺激董朝，其受到的傷害程度將直接與刺激強(qiáng)度成正比鸠项。
Random network (a) and scale-free network (b). In the scale-free network, the larger hubs are highlighted. 圖片來源：[http://en.wikipedia.org/wiki/File:Scale-free_network_sample.png](http://en.wikipedia.org/wiki/File:Scale-free_network_sample.png)有了這一理論基礎(chǔ)，我們可以嘗試一系列權(quán)函數(shù)的參數(shù)值子姜，如
然后找出得到的網(wǎng)絡(luò)最符合scale-free network的frequency distribution的這樣一個
作為后續(xù)分析使用祟绊。但在實(shí)際尋參過程中楼入，存在一個trade-off，即在最大化topology free network的回歸系數(shù)R2
與保證節(jié)點(diǎn)的關(guān)聯(lián)數(shù)目上為此消彼長的關(guān)系牧抽，這其實(shí)可以建立一個優(yōu)化模型嘉熊，但是文章的作者沒有繼續(xù)深入研究得到客觀的尋參標(biāo)準(zhǔn)。盡管作者已經(jīng)建立了一個scale-free topology criterion扬舒，但分析過程依舊有很大的主觀成分阐肤。 Post-module analysis 當(dāng)建立完module后，為了方便處理module與其它數(shù)據(jù)的相關(guān)性讲坎，如性狀信息孕惜，有必要在每一個模塊中定義一個特征基因，這一基因能在可接受的信息損失程度下代表這一module的特征晨炕，這樣做的一個很大好處是能簡化計(jì)算衫画，即使在數(shù)據(jù)量極為龐大的時候也能快速地得到結(jié)果。 在后續(xù)的分析中瓮栗，作者還比較了hard/soft-threshold方法建立網(wǎng)絡(luò)的clustering coefficient的相關(guān)性質(zhì)削罩，以及它們對網(wǎng)絡(luò)connectivity的影響，作者這樣分析是為了說明soft-threshold方法相比hard-threshold的優(yōu)勢费奸，由于涉及到較深的圖論方面的知識弥激，并且與module的建立無關(guān)，因此不在此處進(jìn)一步討論货邓。 參考文獻(xiàn)： Bin Zhang, Steve Horvath, A General Framework for Weighted Gene Co-Expression Network Analysis, Statistical Applications in Genetics and Molecular Biology, Departments of Human Genetics and Biostatistics, University of California at Los Angeles, 2005, Volume 4, issue 1, Article 17.
