正文

正文部分，我想盡可能的從原理出發(fā)，講解WGCNA背后的大致邏輯。我盡可能的在這部分不介紹代碼，以免喧賓奪主般的干擾我們對(duì)于WGCNA原理的思考。

先簡(jiǎn)單給個(gè)目錄，WGCNA背后邏輯得講解將從一下幾個(gè)部分進(jìn)行：

• 數(shù)據(jù)過(guò)濾
• 構(gòu)建共表達(dá)網(wǎng)絡(luò)
  • 構(gòu)建相似度矩陣
  • 構(gòu)建鄰接矩陣的參數(shù)選擇
  • 構(gòu)建拓?fù)渲丿B矩陣（TOM）
  • 利用TOM進(jìn)行聚類
  • 識(shí)別基因模塊
  • 合并相似模塊
  • 模塊功能的驗(yàn)證
• 模塊與外部性轉(zhuǎn)的關(guān)聯(lián)
• R中的WGCNA可視化
• Cytoscape的網(wǎng)絡(luò)可視化

數(shù)據(jù)過(guò)濾

對(duì)于想要分析的數(shù)據(jù)，我們要保證缺失值不能過(guò)多，如果過(guò)多我們必須將這些數(shù)據(jù)去除。好在我們不用手動(dòng)去除，只需要交給WGCNA的一個(gè)函數(shù)即可，這樣可以保證數(shù)據(jù)去除的標(biāo)準(zhǔn)更加客觀。

此外，對(duì)于多個(gè)樣本的數(shù)據(jù)，很有可能某幾個(gè)樣本與其他樣本的差異太大，這將嚴(yán)重影響我們后續(xù)的分析。因此在前期數(shù)據(jù)過(guò)濾部分，我們首先要進(jìn)行聚類，查看一下哪些樣本不符合我們的預(yù)期，提前將這些樣本去除。下面是一個(gè)簡(jiǎn)單的聚類結(jié)果

我一直覺(jué)得聚類可以單獨(dú)寫(xiě)一個(gè)教程具體講解一下，起碼了解一下常見(jiàn)的聚類方法和原理

image

其中的紅線，將樣本中的outlier與其他樣本區(qū)分開(kāi)。這個(gè)cut-off的選取可以自己設(shè)置，將自己認(rèn)為的outlier排除在外即可。

我相信即使不懂聚類的原理，也能看出來(lái)左上的那個(gè)分支所代表的樣本，和其他樣本相差甚遠(yuǎn)

構(gòu)建共表達(dá)調(diào)控網(wǎng)絡(luò)

這一部分是WGCNA的重頭戲，不得不好好說(shuō)道一下……

構(gòu)建相似度矩陣

一個(gè)網(wǎng)絡(luò)關(guān)系通常是由一個(gè)鄰接矩陣所定義的，在這個(gè)矩陣中，每個(gè)元素的取值范圍都在0-1之間，代表兩個(gè)節(jié)點(diǎn)之間的連接強(qiáng)度。比如矩陣中第i行，第j列的元素我們稱之為a ，它代表的就是gene i 與gene j之間的連接強(qiáng)度。

為了得到這個(gè)鄰接矩陣，我們需要定義一個(gè)中間變量，共表達(dá)相似性矩陣。以下圖為例：

首先，我們定義一個(gè)，3樣本10基因的表達(dá)矩陣。每一行是一個(gè)基因，每一列代表一個(gè)樣本，每個(gè)元素代表特定樣本中某基因的表達(dá)值。

有了該矩陣，我們計(jì)算Pearson相關(guān)系數(shù)，從而得到了一個(gè)10X10的相似度矩陣。這個(gè)矩陣中每個(gè)元素都代表著兩個(gè)基因之間的相似度，如第三行第四列的元素值為0.99，代表Gene 3和Gene 4的相似度是0.99。

使用皮爾遜相關(guān)系數(shù)具有一個(gè)缺點(diǎn)，對(duì)于異常值十分的敏感。為了解決該問(wèn)題，我們通常采用Jackknife相關(guān)系數(shù)，即先使用Jackknife方法進(jìn)行重采樣，再進(jìn)行相關(guān)性的分析（先抽樣，再相關(guān)性分析）

圖示中的例子采用的是硬閾值（hard thresholding）的方式構(gòu)建的鄰接矩陣，即相似度超過(guò)0.8的才認(rèn)為這倆基因之間存在連接，設(shè)置為1。而不超過(guò)0.8的則表示二者沒(méi)有連接，設(shè)置為0。

硬閾值就是我們常用的一刀切的閾值

我們需要注意的是，這里定義相似度矩陣時(shí)，使用的公式是 |r(Gi, Gj)| 。即，對(duì)相關(guān)系數(shù)進(jìn)行了絕對(duì)值的計(jì)算，這種方式得到的是unsigned的，即不知道是正相關(guān)還是負(fù)相關(guān)。

公式不重要，重要的是這里有一個(gè)絕對(duì)值符號(hào)；r(Gi, Gj)的取值范圍為-1到1

除此之外，我們還可以使用這個(gè)公式 |(1+ r(Gi, Gj))/2| 進(jìn)行計(jì)算。這樣的話，我們就可以將負(fù)值（負(fù)相關(guān)）轉(zhuǎn)換為正值，而原有的正值變得更大。利用這種方式，在后續(xù)分析的時(shí)候，我們更為關(guān)注的將是那些正相關(guān)的基因。

正相關(guān)反映到基因調(diào)控中，代表正向調(diào)控，gene A上調(diào)會(huì)導(dǎo)致gene B跟著上調(diào)

構(gòu)建鄰接矩陣的參數(shù)選擇

最簡(jiǎn)單的一個(gè)方法，就是如上圖中的例子一樣，選擇一個(gè)閾值（thresholding），然后一刀切。超過(guò)這個(gè)閾值的，我們認(rèn)為這倆基因存在連接，低于這個(gè)閾值的我們認(rèn)為這倆基因不存在鏈接。但是這種定義比較粗暴，比如，我定義一個(gè)cutoff為0.8，那么0.79算不算兩個(gè)基因上存在連接呢？

上面那種閾值的定義方法，我們稱之為硬閾值（hard thresholding）。為了避免硬閾值所帶來(lái)的問(wèn)題，WGCNA分析構(gòu)建共表達(dá)網(wǎng)絡(luò)時(shí)，使用軟閾值（soft thresholding）。軟閾值簡(jiǎn)單來(lái)說(shuō)，就是將相似度進(jìn)行一個(gè)冪運(yùn)算，從而使得結(jié)果符合無(wú)尺度網(wǎng)絡(luò)的標(biāo)準(zhǔn)。

該方法讓原本值大的更大，小的更小。

無(wú)尺度網(wǎng)絡(luò)模型具有一個(gè)特點(diǎn)就是對(duì)于部分節(jié)點(diǎn)的錯(cuò)誤，具有很高程度的魯棒性 （robustness，可以簡(jiǎn)單的理解為容錯(cuò)程度）。這是因?yàn)樵摼W(wǎng)絡(luò)具有很多的冗余鏈接，部分關(guān)鍵節(jié)點(diǎn)的錯(cuò)誤并不會(huì)導(dǎo)致整個(gè)網(wǎng)絡(luò)的調(diào)控出現(xiàn)問(wèn)題。這個(gè)特點(diǎn)，與生物體內(nèi)復(fù)雜的網(wǎng)絡(luò)十分相似，這也是為什么使用無(wú)尺度模型的另一個(gè)重要原因。

在無(wú)尺度網(wǎng)絡(luò)中，節(jié)點(diǎn)的度與該節(jié)點(diǎn)出現(xiàn)的概率是負(fù)相關(guān)的。比如：具有五個(gè)連接的節(jié)點(diǎn)，其出現(xiàn)的概率要小于有三個(gè)連接的節(jié)點(diǎn)出現(xiàn)的概率。如下圖：

image

其中，橫坐標(biāo)表示某個(gè)節(jié)點(diǎn)連通性的log值，縱坐標(biāo)表示該節(jié)點(diǎn)出現(xiàn)的概率的log值。當(dāng)我們選擇的冪指數(shù) (power)，讓數(shù)據(jù)十分符合無(wú)尺度網(wǎng)絡(luò)時(shí)，其R 就會(huì)十分的接近于1。在這張圖上，分布大致遵循一條直線，這被稱為近似無(wú)尺度拓?fù)洹?/p>

連通性，簡(jiǎn)單的理解為節(jié)點(diǎn)上連接的邊的個(gè)數(shù)

我們?cè)谑褂肳GCNA分析時(shí)，通常會(huì)產(chǎn)生一個(gè)如下所示的圖輔助我們進(jìn)行冪指數(shù) (power) 的選擇：

冪的英文即power，有的教程中，會(huì)使用β指代該值。

image

兩張圖的橫坐標(biāo)，均代表冪指數(shù) (power) 的取值。左圖的縱坐標(biāo)是我們上面所講的R ，越高越好。而右圖的縱坐標(biāo)代表所有節(jié)點(diǎn)的平均連通數(shù)，越低越好。通常我們對(duì)于冪指數(shù) (power) 的選擇主要參考左圖，右圖用于輔助參考。

無(wú)尺度網(wǎng)絡(luò)中，大部分節(jié)點(diǎn)的連通度較低，因此平均連通數(shù)越低越好。

此外，在所有達(dá)到設(shè)定閾值的冪指數(shù)中，冪指數(shù)越小越好，不然我們直接選擇1000，或者一個(gè)非常大的值它不香么？

構(gòu)建拓?fù)渲丿B矩陣(TOM)

WGCNA認(rèn)為鄰接矩陣是不夠的，基因間的相似性應(yīng)該在表達(dá)和網(wǎng)絡(luò)拓?fù)渌缴象w現(xiàn)。于是，又利用鄰接矩陣生成了一個(gè)拓?fù)渲丿B矩陣/TOM（topological overlap matrix），此矩陣的構(gòu)建還能夠盡可能地最小化噪音和假陽(yáng)性的影響。

gene A 與 gene B之間的關(guān)系，不僅考慮AB直接的關(guān)系，還考慮可能的多個(gè)間接關(guān)系，這樣任何一個(gè)關(guān)系存在假陽(yáng)性都不至于對(duì)AB之間的關(guān)系造成巨大的影響。有點(diǎn)無(wú)尺度網(wǎng)絡(luò)那味了~

比如gene A 與 gene B，鄰接矩陣關(guān)注的是這倆基因之間直接的共表達(dá)關(guān)系，而拓?fù)渲丿B矩陣還考慮中間因素的影響，不止關(guān)注gene A對(duì)B的調(diào)控，還要關(guān)注geneA對(duì)B的次級(jí)影響。比如gene A調(diào)控 gene C，gene C又調(diào)控gene B。

既然考慮了次級(jí)調(diào)控，為什么不考慮更次一級(jí)的調(diào)控呢？這里將其稱之為三級(jí)調(diào)控。如，A不僅自己調(diào)控B，還通過(guò)調(diào)控C，C調(diào)控D，D又調(diào)控B的方式，間接的調(diào)控B。兩個(gè)原因，一個(gè)可能是計(jì)算量太大，太耗費(fèi)計(jì)算資源了。另一個(gè)可能是，三級(jí)調(diào)控對(duì)基因B的影響基本可以忽略不計(jì)了，這樣也就可以解釋為什么沒(méi)有四級(jí)調(diào)控，五級(jí)調(diào)控……

那么拓?fù)渲丿B矩陣怎么計(jì)算得到呢？具體看如下公式：

image.png

看的我腦袋疼，相信你們腦袋也疼。所以我們還是舉個(gè)例子把，這樣更方便大家的理解。具體如下圖：

image

鄰接矩陣中，A與D的關(guān)系是1。而在拓?fù)渲丿B矩陣中，A和D的關(guān)系，不僅需要考慮A和D之間的虛線，還需要A-B-D以及A-C-D。（通過(guò)B與C的兩條間接關(guān)系）

ABBC=11,=1 ACCD=11=1，因此公式中l(wèi) =1+1 = 2， a =1。l + a 就是A 到D 的直接和間接共表達(dá)的總和。而我們可以知道A的k是3 (k )，D的k也是3，因此min{k ,k }=3，因此最終w 的結(jié)果為（2+1）/（3+1-1）=1，即我們得到的拓?fù)渲丿B矩陣中，第i行j列的元素值為1。

利用TOM進(jìn)行聚類

有了拓?fù)洚悩?gòu)矩陣（TOM）之后，我們需要進(jìn)行聚類，從而對(duì)基因分塊，得到不同的模塊（module）。為了更好的展示聚類結(jié)果，我們用1減去TOM，從而得到對(duì)應(yīng)的基因不相似性的矩陣（dissTOM）。

接著利用該不相似性的矩陣，進(jìn)行聚類，聚類結(jié)果示例如下：

WGCNA直接將1-TOM得到的矩陣dissTOM，強(qiáng)制轉(zhuǎn)換為距離矩陣。因此越相關(guān)的兩個(gè)基因，其dissTOM中值越小，距離矩陣中對(duì)應(yīng)的距離越近，在圖中表示就是越容易聚到一起。

image

在這張圖上，樹(shù)上的每一個(gè)葉子節(jié)點(diǎn)都是一條小豎線，代表著一個(gè)基因。左側(cè)的縱坐標(biāo)（height）表示兩個(gè)子節(jié)點(diǎn)的距離，如gene A 與 gene B 在0.75處merge到了一起，則這倆基因的距離為0.75，不相似度為0.75。

識(shí)別基因模塊

得到聚類結(jié)果之后，接下來(lái)要做的就是找模塊（module），即那些連接十分緊密的基因集（共表達(dá)程度很高的一類基因），聚類后的每一個(gè)簇代表著一個(gè)module。

通常將模塊樹(shù)的分支切分開(kāi)的方法有兩種，一種是指定一個(gè)高度，將結(jié)果分成若干個(gè)模塊（即：數(shù)據(jù)預(yù)處理部分的紅線）。另一種是使用 dynamic branch cut methods，該方法的優(yōu)點(diǎn)是可以自動(dòng)化的將module區(qū)分開(kāi)，不再需要手動(dòng)指定高度，且該方法更為靈活。

使用Dynamic Branch Cut方法處理得到的結(jié)果如下：
!](https://cdn.nlark.com/yuque/0/2020/png/1111451/1607994362237-dd3d8d10-6387-4ba7-a026-9b2a3c6cef13.png#align=left&display=inline&height=1115&margin=%5Bobject%20Object%5D&name=image.png&originHeight=1115&originWidth=1499&size=240716&status=done&style=none&width=1499)

合并相似模塊

Dynamic Branch Cut作為一種自動(dòng)切分模塊的方法，有時(shí)可能會(huì)識(shí)別出來(lái)兩個(gè)表達(dá)譜非常相似的模塊，此時(shí)就需要我們手動(dòng)將這些高度共表達(dá)的模塊合并，方便后續(xù)分析。首先需要對(duì)不同的模塊進(jìn)行聚類：

Dynamic Branch Cut利用形狀參數(shù)確定分類的個(gè)數(shù)

image

這里選擇的高度為0.75，對(duì)應(yīng)的相關(guān)性為0.75。表示，我們要將相關(guān)性大于等于0.75的模塊進(jìn)行合并，最終合并結(jié)果如下：

image.png

認(rèn)真看，可以觀察到，紫色和黃色兩個(gè)色塊合并到了一起，綠色和紅色兩種module合并到了一起。

模塊功能的驗(yàn)證

層次聚類的一個(gè)缺點(diǎn)是很難確定究竟應(yīng)該將整個(gè)數(shù)據(jù)集分成幾類，盡管我們采用多種分支剪切/模塊識(shí)別的方法，但是數(shù)據(jù)的分類問(wèn)題，是一個(gè)開(kāi)放式的探索問(wèn)題，并沒(méi)有一個(gè)固定解。即：合理即可。

那么怎么劃分模塊才是比較合理的呢？通常來(lái)說(shuō)一個(gè)共表達(dá)模塊反應(yīng)著一個(gè)真實(shí)的生物學(xué)信號(hào)，比如通路，比如某個(gè)生物學(xué)過(guò)程?；蛘叻磻?yīng)著噪音，如組織污染，技術(shù)偏差，假陽(yáng)性等。

為了判斷一個(gè)模塊的劃分是否合理，通常我們會(huì)對(duì)一個(gè)模塊進(jìn)行GO富集分析，判斷是否富集在某一個(gè)或者某幾個(gè)通路。

模塊與外部性狀關(guān)聯(lián)

在WGCNA中，我們將劃分得到的模塊，與外部性狀相關(guān)聯(lián)，來(lái)查看每個(gè)模塊分別與什么性狀高度相關(guān)，從而找到我們感興趣性狀的對(duì)應(yīng)模塊。這些性狀可以是各種指標(biāo)，比如小鼠的重量，尾長(zhǎng)，毛色等；也可以是對(duì)應(yīng)植物的高矮，果實(shí)多少，顆粒飽滿度等。

但是我們知道計(jì)算顯著性，通常是兩個(gè)向量之間作比較。但是一個(gè)模塊是一個(gè)矩陣，矩陣和向量之間，計(jì)算相關(guān)性，我們?cè)撊绾畏治瞿?？WGCNA選擇使用PCA的方法，計(jì)算模塊矩陣的主成分，并將其中的PC1定義為eigengenes。

有了eigengennes 我們就可以計(jì)算模塊與外部性狀的相關(guān)性了?？梢暬Y(jié)果如下：

image

熱圖中的每個(gè)元素代表指定模塊和對(duì)應(yīng)性狀的顯著性，我們重點(diǎn)關(guān)注那些高度顯著的元素，即紅色的方塊。

模塊與指定性狀的相關(guān)性我們已經(jīng)分析得到了，但是一個(gè)模塊通常有很多基因。我們?nèi)绾握业阶钕嚓P(guān)的基因呢？WGCNA通過(guò)GS和MM給出了一個(gè)合理的解決辦法。

首先我們將模塊中的每個(gè)基因與eigengenes進(jìn)行相關(guān)性分析，得到的結(jié)果我們稱之為module membership（MM）。當(dāng)結(jié)果越接近于0，則我們認(rèn)為該基因與其所在的模塊越不相關(guān)。而結(jié)果越接近于1或者-1，則我們認(rèn)為該基因與模塊基因高度相關(guān)，正負(fù)表示其是正相關(guān)還是負(fù)相關(guān)。

此外，module membership與模塊內(nèi)的連通性是高度相關(guān)的。其中，高連通性的hub genes傾向于有更高的module membership 值。有了這個(gè)因果關(guān)系，我們不需要通過(guò)cytoscape繪制網(wǎng)絡(luò)圖，也可以找到hub genes了。

模塊內(nèi)的基因，我們不僅想知道它與模塊的相關(guān)性，還想知道它與模塊對(duì)應(yīng)的性狀的相關(guān)性。于是WGCNA定義了一個(gè)新的變量GS (gene significance)，即計(jì)算模塊內(nèi)gene 與對(duì)應(yīng)性狀的顯著性。下圖展示了 module membership的相關(guān)性散點(diǎn)圖：

image

從上圖我們可以看到，GS和MM還是很相關(guān)的，這表明與性狀高度相關(guān)的基因，通常也是該性狀對(duì)應(yīng)模塊內(nèi)比較重要的基因。因此當(dāng)我們選擇感興趣的重要基因時(shí)，推薦選取散點(diǎn)圖右上角部分的基因。

R中的WGCNA可視化

我們通常使用熱圖展示我們得到的權(quán)重網(wǎng)絡(luò)，其中每一行和每一列分別代表著一個(gè)基因。使用熱圖可以給我們展示哪些基因之間更為相似，通常來(lái)說(shuō)都是一個(gè)模塊內(nèi)的基因相似度更高。如下圖：

熱圖中方塊的顏色越深（紅）表示共表達(dá)相關(guān)性越高，越淺（黃）表示相關(guān)性越弱

image

除了看基因之間的共表達(dá)程度，我們還可以看看不同模塊之間的關(guān)系。這時(shí)候，就需要我們之前計(jì)算得到的eigengenes了，利用eigengenes我們可以得到不同模塊的聚類，可視化不同模塊的相關(guān)性，如下圖：

image

熱圖中，添加了weight這一外部信息，讓我們看到該屬性最相關(guān)的module分別是哪些。

Cytoscape的網(wǎng)絡(luò)可視化

最后一步，也是大家最為喜聞樂(lè)見(jiàn)的一步，就是將我們得到的模塊結(jié)果保存，并載入到cytoscape中進(jìn)行網(wǎng)絡(luò)的可視化。每個(gè)頂點(diǎn)代表一個(gè)基因，每條邊代表兩個(gè)基因之間存在共表達(dá)關(guān)系。

Cytoscape如果講起來(lái)，那將又是一堆字，還是再單獨(dú)開(kāi)個(gè)教程講解吧

image

完結(jié)撒花，原理/思路部分告一段落，下一部分就是WGCNA的實(shí)戰(zhàn)部分了。實(shí)戰(zhàn)部分將會(huì)跟著原理部分一步一步走下來(lái)，分別講述對(duì)應(yīng)的代碼含義和結(jié)果。如果不想等待，而且不反感英語(yǔ)，強(qiáng)烈建議看看這個(gè)英文教程：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/

參考資料

• https://abego.cn/2019/05/31/what-is-wgcna-when-to-use/
• https://zh.wikipedia.org/wiki/%E5%9F%BA%E5%9B%A0%E5%85%B1%E8%A1%A8%E8%BE%BE%E7%BD%91%E7%BB%9C
• https://zh.wikipedia.org/wiki/%E5%8A%A0%E6%9D%83%E7%9B%B8%E5%85%B3%E7%BD%91%E7%BB%9C%E5%88%86%E6%9E%90
• https://zhuanlan.zhihu.com/p/36132370
• http://www.wuchangsong.com/?p=198
• http://www.sci666.net/58128.html
• https://blog.csdn.net/weixin_44377608/article/details/90489878
• http://m.itdecent.cn/p/fe4d3c77786e
• http://www.sci666.net/58731.html
• [http://www.biotrainee.com/thread-646-1-1.html](