聽說(shuō)WGCNA是一個(gè)看上去比較厲害的轉(zhuǎn)錄組分析方法，近期又多次看到了相關(guān)的內(nèi)容，所以跟風(fēng)學(xué)習(xí)一下。對(duì)于我這種業(yè)余學(xué)生信的來(lái)說(shuō)確實(shí)復(fù)雜了些，斷斷續(xù)續(xù)看了很久，也只是勉強(qiáng)按流程跟著跑一遍。這篇只有簡(jiǎn)單的原理和介紹，沒有代碼。
一、WGCNA簡(jiǎn)介
對(duì)于兩組數(shù)據(jù)，比如實(shí)驗(yàn)組和對(duì)照組的測(cè)序數(shù)據(jù)，如果我們要找到與實(shí)驗(yàn)組條件相關(guān)的通路或者基因，可以做差異分析然后進(jìn)行GO分析、KEGG分析，或者直接使用GSEA分析。那么如果我有多個(gè)分組怎么找出每個(gè)組的特異分子呢？假設(shè)我有10個(gè)分組共60個(gè)樣品，如果兩兩比較需要比45次，這樣不夠高效。實(shí)際上臨床性狀可能更加復(fù)雜，有多種性狀混合。
加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（Weighted Gene Co-Expression Network Analysis，WGCNA），是一種適合進(jìn)行多樣本復(fù)雜數(shù)據(jù)分析的工具，通過計(jì)算基因間表達(dá)關(guān)系，鑒定表達(dá)模式相似的基因集合，解析基因集合與樣品表型之間的聯(lián)系，繪制基因集合中基因之間的調(diào)控網(wǎng)絡(luò)并鑒定關(guān)鍵調(diào)控基因。WGCNA將復(fù)雜生物過程的基因共表達(dá)網(wǎng)絡(luò)劃分為高度相關(guān)的幾個(gè)特征模塊，其代表著幾組高度協(xié)同變化的基因集，并可將模塊與特定的臨床特征建立關(guān)聯(lián)，從中尋找發(fā)揮關(guān)鍵功能的基因，幫助識(shí)別參與特定生物學(xué)過程的潛在機(jī)制以及探索候選生物標(biāo)志物。
相比于只關(guān)注差異表達(dá)的基因，WGCNA利用數(shù)千或近萬(wàn)個(gè)變化最大的基因或全部基因的信息識(shí)別感興趣的基因集，并與表型進(jìn)行顯著性關(guān)聯(lián)分析。一是充分利用了信息，二是把數(shù)千個(gè)基因與表型的關(guān)聯(lián)轉(zhuǎn)換為數(shù)個(gè)基因集與表型的關(guān)聯(lián)，免去了多重假設(shè)檢驗(yàn)校正的問題。??
二、WGCNA原理
（1）基本假設(shè)
WGCNA基于兩個(gè)假設(shè)：（1）相似表達(dá)模式的基因可能存在共調(diào)控、功能相關(guān)或處于同一通路，（2）基因網(wǎng)絡(luò)符合無(wú)尺度分布。
【一個(gè)網(wǎng)絡(luò)圖主要包含節(jié)點(diǎn)（node）和連接線（edge）兩種成分，如果兩個(gè)節(jié)點(diǎn)之間存在相互關(guān)系，就畫一條線連接起來(lái)。網(wǎng)絡(luò)可以分為有方向的和沒有方向的，對(duì)于基因網(wǎng)絡(luò)而言，一般我們預(yù)測(cè)到的都是沒有方向性的。度（degree）代表著一個(gè)節(jié)點(diǎn)和多少個(gè)其他節(jié)點(diǎn)的有連接，如果有方向就會(huì)分為in-degree和out-degree。隨機(jī)網(wǎng)絡(luò)（random network）是指在網(wǎng)絡(luò)中的所有點(diǎn)之間的關(guān)聯(lián)性都差不多，可以看做是隨機(jī)連接而成。無(wú)尺度網(wǎng)絡(luò)（scale-free network）具有嚴(yán)重的異質(zhì)性，其各節(jié)點(diǎn)之間的連接狀況（度數(shù)）具有嚴(yán)重的不均勻分布性：網(wǎng)絡(luò)中少數(shù)稱之為Hub點(diǎn)的節(jié)點(diǎn)擁有極其多的連接，而大多數(shù)節(jié)點(diǎn)只有很少量的連接。少數(shù)Hub點(diǎn)對(duì)無(wú)標(biāo)度網(wǎng)絡(luò)的運(yùn)行起著主導(dǎo)的作用。從廣義上說(shuō)，無(wú)標(biāo)度網(wǎng)絡(luò)的無(wú)標(biāo)度性是描述大量復(fù)雜系統(tǒng)整體上嚴(yán)重不均勻分布的一種內(nèi)在性質(zhì)?！?/p>
（2）計(jì)算兩個(gè)基因的相關(guān)系數(shù)
傳統(tǒng)方法上，描述兩個(gè)基因間的關(guān)聯(lián)程度，可通過計(jì)算表達(dá)值間的Pearson、Spearman等相關(guān)系數(shù)獲得。為了構(gòu)建關(guān)聯(lián)網(wǎng)絡(luò)，通常指定一個(gè)篩選閾值，如相關(guān)系數(shù)大于0.8以上，作為兩個(gè)基因間具有強(qiáng)關(guān)聯(lián)程度的依據(jù)。隨后在網(wǎng)絡(luò)中，節(jié)點(diǎn)表示基因，若兩個(gè)基因間的相關(guān)系數(shù)大于指定的閾值，則會(huì)以邊連接以表示兩者間可能存在相互作用。但是基于固定閾值法的缺點(diǎn)在于，閾值是人為定義的，將會(huì)忽略很多潛在關(guān)聯(lián)。例如，0.79就是不相關(guān)嗎？同時(shí)，這種一刀切的方法也會(huì)丟失基因的變化趨勢(shì)信息，將難以在網(wǎng)絡(luò)中描述相關(guān)性的強(qiáng)弱關(guān)系。
經(jīng)典的共表達(dá)網(wǎng)絡(luò)使用共表達(dá)相關(guān)性構(gòu)建網(wǎng)絡(luò)。相關(guān)性測(cè)量的值介于-1（完全負(fù)相關(guān)）到1（完全正相關(guān)）之間。在無(wú)方向網(wǎng)絡(luò)（unsigned）中，使用的是絕對(duì)相關(guān)系數(shù)，兩個(gè)負(fù)相關(guān)的基因也會(huì)被認(rèn)為共表達(dá)，這就導(dǎo)致負(fù)相關(guān)的基因也會(huì)被group到一起。有方向的網(wǎng)絡(luò)（signed）采取的是把相關(guān)性值scale到0和1之間，其中0-0.5代表負(fù)相關(guān)，0.5-1代表正相關(guān)。共表達(dá)網(wǎng)絡(luò)定義為無(wú)向的、加權(quán)的基因網(wǎng)絡(luò)。
【Pearson相關(guān)系數(shù)是用于表示相關(guān)性大小的最常用指標(biāo)，數(shù)值介于-1~1之間，越接近0相關(guān)性越低，越接近-1或1相關(guān)性越高。正負(fù)號(hào)表明相關(guān)方向，正號(hào)為正相關(guān)、負(fù)號(hào)為負(fù)相關(guān)。適用于兩個(gè)正態(tài)分布的連續(xù)變量。兩個(gè)變量X和Y之間的皮爾森相關(guān)系數(shù)是這樣計(jì)算的：先按基因進(jìn)行z-score標(biāo)準(zhǔn)化 (x-μ)/SD，然后用下列公式計(jì)算 兩個(gè)變量的協(xié)方差/標(biāo)準(zhǔn)差乘積?！?/p>
（3）對(duì)相關(guān)系數(shù)加權(quán)，獲得鄰接矩陣
為了解決這些問題，提出了加權(quán)的思想，WGCNA的做法是對(duì)基因表達(dá)值之間的相關(guān)系數(shù)取β次冪。直接結(jié)果是把基因間相關(guān)性的強(qiáng)弱的差異放大，使強(qiáng)弱關(guān)系更為分明（弱化弱連接，強(qiáng)化強(qiáng)連接），有利于后續(xù)聚類（模塊）識(shí)別（使得相關(guān)性數(shù)值更符合無(wú)尺度網(wǎng)絡(luò)特征，更具有生物意義）。冪次β就是軟閾值（power）。對(duì)于基因i和j，相關(guān)系數(shù)為sij，取β次冪后獲得aij = sij^β，aij表示網(wǎng)絡(luò)中基因節(jié)點(diǎn)i和j之間的連接強(qiáng)度。因此，最終網(wǎng)絡(luò)中基因間的相關(guān)強(qiáng)度，取決于β次冪的選擇，其取值的高低直接影響模塊的構(gòu)建和模塊內(nèi)基因的劃分。β值根據(jù)接近無(wú)標(biāo)度網(wǎng)絡(luò)（scale-free network）的最低值確定，pickSoftThreshold這個(gè)函數(shù)所做的就是確定合適的power，power值的選擇即為確定β值的過程。計(jì)算所有基因間表達(dá)水平的關(guān)聯(lián)強(qiáng)度（加權(quán)相關(guān)性值），獲得鄰接矩陣（Adjacency matrix）。
（4）計(jì)算拓?fù)渲丿B矩陣
為了保證信息的充分利用，進(jìn)一步計(jì)算了拓?fù)渲丿B矩陣（Topological overlap matrix，TOM），以降低噪音和假相關(guān)，獲得的新的相關(guān)性矩陣。其思想是，任何兩個(gè)基因的相關(guān)性不僅僅由它們的表達(dá)相似性直接決定，并且如果基因i和j有很多相同的鄰接基因，意味著兩基因間有相似的表達(dá)模式，TOM ij就比較高。因此用TOM矩陣用于反映基因間表達(dá)的相似性。
（5）模塊識(shí)別
通過基因之間的加權(quán)相關(guān)系數(shù)（TOM矩陣）構(gòu)建分層聚類樹（平均連接層級(jí)聚類），聚類樹的不同分支為模式相似的基因，代表不同的基因模塊。將沒有被分配到特定模塊的基因群定義為灰色模塊。這樣就可以將幾萬(wàn)個(gè)基因通過基因表達(dá)模式被分成了幾十個(gè)模塊，是一個(gè)提取歸納信息的過程。模塊（Module）是指相關(guān)性非常高的基因集。在無(wú)向網(wǎng)絡(luò)中，模塊內(nèi)是高度相關(guān)的基因。在有向網(wǎng)絡(luò)中，模塊內(nèi)是高度正相關(guān)的基因。
Module eigengene是給定模型的第一主成分，代表整個(gè)模型的基因表達(dá)譜，用一個(gè)向量代替了一個(gè)矩陣，方便后期計(jì)算。把每個(gè)模塊里面的基因進(jìn)行分析，整合得到一個(gè)值，用這個(gè)值代替整個(gè)模塊。表達(dá)矩陣原本的成千上萬(wàn)行，就變成了每個(gè)模塊的eigengene的行，矩陣就變的簡(jiǎn)單了。
Connectivity（連接度）類似于網(wǎng)絡(luò)中“度”（degree）的概念，每個(gè)基因的連接度是與其相連的基因的邊屬性之和。在共表達(dá)網(wǎng)絡(luò)中，連接度衡量一個(gè)基因與所有其他網(wǎng)絡(luò)基因的相關(guān)性。Intramodular connectivity（模塊內(nèi)連接度）衡量給定基因相對(duì)于特定模塊的基因的連接或共表達(dá)程度。
Module Membership(MM) 對(duì)于每個(gè)基因，我們通過將基因表達(dá)譜與某模塊的eigengene相關(guān)性來(lái)定義Module Membership。例如，基因i與藍(lán)色模塊eigengene的相關(guān)性MM blue(i)，如果接近0則基因i不是藍(lán)色模塊的一部分，接近1或-1則它與藍(lán)色模塊基因高度正相關(guān)或負(fù)相關(guān)。Gene significance(GS)是基因i與臨床性狀的相關(guān)性，GS越高則基因i的生物學(xué)意義越大。關(guān)鍵基因（Hub gene）是連接度最多或連接多個(gè)模塊的基因。
（6）術(shù)語(yǔ)
（7）整體流程
WGCNA旨在尋找協(xié)同表達(dá)的基因模塊，并探索基因網(wǎng)絡(luò)與關(guān)注的表型之間的關(guān)聯(lián)關(guān)系，以及網(wǎng)絡(luò)中的核心基因。
整體分析流程為：
1. 構(gòu)建基因共表達(dá)網(wǎng)絡(luò)：計(jì)算加權(quán)的表達(dá)相關(guān)性（TOM矩陣），即前面的1-4步
2. 識(shí)別基因集（模塊）：平均連接層級(jí)聚類，根據(jù)設(shè)定標(biāo)準(zhǔn)切分聚類結(jié)果，獲得不同的基因模塊，用聚類樹的分枝和不同顏色表示
得到模塊之后的分析有：
3. 根據(jù)TOM矩陣，繪制相關(guān)性圖
4. 模塊的功能富集：利用 WGCNA 算法得到的基因共表達(dá)模塊，需要進(jìn)一步的功能注釋以發(fā)現(xiàn)其生物學(xué)意義
5. 計(jì)算基因模塊與表型的相關(guān)性，鑒定性狀相關(guān)的模塊
6. 找到模塊的核心基因，選擇感興趣的關(guān)鍵基因
7. 根據(jù)模型中已知基因的功能推測(cè)未知基因的功能
三、沒有具體流程只有圖
（1）樣本聚類樹和臨床trait heatmap
輸入的基因表達(dá)矩陣，如果是芯片數(shù)據(jù)，那么常規(guī)的歸一化矩陣即可；如果是轉(zhuǎn)錄組數(shù)據(jù)，最好是FPKM值或者其它歸一化好的表達(dá)量，并進(jìn)行l(wèi)og2(x+1)轉(zhuǎn)換，也可以使用DESeq2的rlog轉(zhuǎn)換的數(shù)據(jù)?？梢赃M(jìn)行一些前處理，比如處理缺失值，按表達(dá)水平過濾低表達(dá)的基因等，避免它們的影響，同時(shí)減少后續(xù)運(yùn)算資源消耗（只要能保留感興趣的重要基因就可以）。
第一次求相關(guān)性的時(shí)候，其實(shí)跟性狀是沒有關(guān)系的；得到模塊和eigengene后第二次相關(guān)性才跟性狀有關(guān)。常規(guī)的，性狀信息應(yīng)該是數(shù)值，比如患者的年齡、生存時(shí)間。如果是分類變量，需要使用1或者0來(lái)表示。超過二分類的信息，可以使用model.matrix()制造一個(gè)設(shè)計(jì)矩陣添加進(jìn)去，這時(shí)候得到的相關(guān)模塊代表該亞組特異的模塊。除了添加常見的各種臨床信息之外，理論上講，只要是一列數(shù)值就可以添加到性狀中，比如基因A的表達(dá)量，這時(shí)候算出來(lái)的跟這個(gè)“性狀”相關(guān)的模塊就代表和基因A相關(guān)的模塊。也可以添加甲基化信息、突變信息、SNP信息等，任何能獲得的信息都可以添加。
（2）確定軟閾值power
上述圖形的橫軸均代表權(quán)重參數(shù)β（軟閾值），左圖縱軸代表對(duì)應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方，即無(wú)標(biāo)度拓?fù)鋽M合指數(shù)R^2（scale free topology fitting index R2，即統(tǒng)計(jì)結(jié)果中的SFT.R.sq列數(shù)值）。相關(guān)系數(shù)的平方越高，說(shuō)明該網(wǎng)絡(luò)越逼近無(wú)尺度分布，R^2一般要到0.85以上。之所以有負(fù)值是因?yàn)橛殖艘粤藄lope列數(shù)值的負(fù)方向，僅關(guān)注正值即可。 右圖的縱軸代表對(duì)應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值（mean.k.網(wǎng)絡(luò)的平均連接度）。 綜合考慮這些指標(biāo)選擇合適的軟閾powers值，使R^2盡可能大但也要保證連通度不能太低?？梢赃x擇R^2大于0.85或平均連接度小于100時(shí)的power值。最佳的beta值就是sft$powerEstimate。如果實(shí)際分析中的縱坐標(biāo)達(dá)不到0.85，軟閾值可選取經(jīng)驗(yàn)值。
（3）構(gòu)建網(wǎng)絡(luò)
基于基因表達(dá)值矩陣計(jì)算基因間表達(dá)的相似度，獲得鄰接度矩陣，并進(jìn)一步計(jì)算拓?fù)渲丿B矩陣TOM（Pearson相關(guān)系數(shù)→鄰接矩陣軟閾值→拓?fù)渲丿B矩陣）。TOM矩陣中的值就反映了兩兩基因間共表達(dá)的相似度，取值0-1，基因間共表達(dá)相似度越高，值越趨于1。然后推出基因間的相異性系數(shù)（1-TOM），并據(jù)此得到基因間的系統(tǒng)聚類樹。然后按照混合動(dòng)態(tài)剪切樹的標(biāo)準(zhǔn)，設(shè)置每個(gè)基因模塊最少的基因數(shù)目為30。再次分析，依次計(jì)算每個(gè)模塊的特征向量值，然后對(duì)模塊進(jìn)行聚類分析，將距離較近的模塊合并為新的模塊（基因按照相異性1-TOM聚類→混合動(dòng)態(tài)剪切樹識(shí)別模塊→得到特征向量→模塊聚類與合并）。經(jīng)過這一步，把輸入的表達(dá)矩陣的幾千個(gè)基因組歸類成了幾十個(gè)模塊。
（4）劃分模塊
（5）共表達(dá)模塊與表型的關(guān)聯(lián)分析
模塊和臨床trait的關(guān)聯(lián)，每行對(duì)應(yīng)一個(gè)模塊ME值，每列對(duì)應(yīng)一個(gè)性狀，每個(gè)單元格包含相應(yīng)的相關(guān)性和p值。這個(gè)圖就是把moduleTraitCor這個(gè)矩陣用熱圖可視化一下。
（6）模塊內(nèi)的關(guān)鍵基因
感興趣的模塊的Module membership和 Gene significance的關(guān)系圖。選取右上角的基因作為感興趣基因，或者模塊內(nèi)的連接度top10作為hub基因。
（7）TOM矩陣可視化
這一步非常耗時(shí)，可以選取其中部分基因作圖。如果3000個(gè)基因保存成pdf，打開時(shí)能把我電腦的8G內(nèi)存占用100%，然后徹底卡死。另外如果基因超過5000個(gè)，默認(rèn)會(huì)被分成幾個(gè)block，這一步有可能報(bào)錯(cuò)。其實(shí)這個(gè)圖并沒有什么用，只是看著好看而已。
（8）臨床trait和ME相關(guān)性
將自己感興趣的臨床trait數(shù)據(jù)納入ME，統(tǒng)一制作ME相關(guān)性的熱圖。最好分開畫吧，兩個(gè)圖一起字總是顯示不全，還要調(diào)半天，比例也不好看。
（9）網(wǎng)絡(luò)的可視化
這一部分留著以后學(xué)了。
最后，感謝生信技能樹的WGCNA直播課，以及一大堆生信的良心公眾號(hào)，對(duì)于初學(xué)者真的很有幫助。
