正文

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

先簡單給個目錄帜消，WGCNA背后邏輯得講解將從一下幾個部分進行：

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

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

對于想要分析的數(shù)據(jù)，我們要保證缺失值不能過多浓体，如果過多我們必須將這些數(shù)據(jù)去除泡挺。好在我們不用手動去除，只需要交給WGCNA的一個函數(shù)即可汹碱，這樣可以保證數(shù)據(jù)去除的標準更加客觀粘衬。

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

我一直覺得聚類可以單獨寫一個教程具體講解一下屯阀，起碼了解一下常見的聚類方法和原理

image

其中的紅線缅帘，將樣本中的outlier與其他樣本區(qū)分開。這個cut-off的選取可以自己設(shè)置难衰，將自己認為的outlier排除在外即可钦无。

我相信即使不懂聚類的原理，也能看出來左上的那個分支所代表的樣本盖袭，和其他樣本相差甚遠

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

這一部分是WGCNA的重頭戲失暂，不得不好好說道一下……

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

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

為了得到這個鄰接矩陣倍踪，我們需要定義一個中間變量系宫，共表達相似性矩陣。以下圖為例：

首先建车，我們定義一個笙瑟，3樣本10基因的表達矩陣。每一行是一個基因癞志，每一列代表一個樣本，每個元素代表特定樣本中某基因的表達值框产。

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

使用皮爾遜相關(guān)系數(shù)具有一個缺點膊存，對于異常值十分的敏感。為了解決該問題忱叭，我們通常采用Jackknife相關(guān)系數(shù)隔崎，即先使用Jackknife方法進行重采樣，再進行相關(guān)性的分析（先抽樣韵丑，再相關(guān)性分析）

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

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

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

公式不重要奔坟，重要的是這里有一個絕對值符號携栋；r(Gi, Gj)的取值范圍為-1到1

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

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

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

最簡單的一個方法咽筋，就是如上圖中的例子一樣溶推，選擇一個閾值（thresholding），然后一刀切奸攻。超過這個閾值的蒜危，我們認為這倆基因存在連接，低于這個閾值的我們認為這倆基因不存在鏈接睹耐。但是這種定義比較粗暴辐赞，比如，我定義一個cutoff為0.8硝训，那么0.79算不算兩個基因上存在連接呢响委？

上面那種閾值的定義方法，我們稱之為硬閾值（hard thresholding）窖梁。為了避免硬閾值所帶來的問題赘风，WGCNA分析構(gòu)建共表達網(wǎng)絡(luò)時，使用軟閾值（soft thresholding）纵刘。軟閾值簡單來說贝次，就是將相似度進行一個冪運算，從而使得結(jié)果符合無尺度網(wǎng)絡(luò)的標準彰导。

該方法讓原本值大的更大蛔翅，小的更小敲茄。

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

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

image

其中，橫坐標表示某個節(jié)點連通性的log值掺出，縱坐標表示該節(jié)點出現(xiàn)的概率的log值徽千。當我們選擇的冪指數(shù) (power)，讓數(shù)據(jù)十分符合無尺度網(wǎng)絡(luò)時汤锨，其R 就會十分的接近于1双抽。在這張圖上，分布大致遵循一條直線闲礼，這被稱為近似無尺度拓撲牍汹。

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

我們在使用WGCNA分析時柬泽，通常會產(chǎn)生一個如下所示的圖輔助我們進行冪指數(shù) (power) 的選擇：

冪的英文即power柑贞，有的教程中，會使用β指代該值聂抢。

image

兩張圖的橫坐標，均代表冪指數(shù) (power) 的取值棠众。左圖的縱坐標是我們上面所講的R 琳疏，越高越好。而右圖的縱坐標代表所有節(jié)點的平均連通數(shù)闸拿，越低越好空盼。通常我們對于冪指數(shù) (power) 的選擇主要參考左圖，右圖用于輔助參考新荤。

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

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

構(gòu)建拓撲重疊矩陣(TOM)

WGCNA認為鄰接矩陣是不夠的，基因間的相似性應(yīng)該在表達和網(wǎng)絡(luò)拓撲水平上體現(xiàn)澄者。于是笆呆，又利用鄰接矩陣生成了一個拓撲重疊矩陣/TOM（topological overlap matrix），此矩陣的構(gòu)建還能夠盡可能地最小化噪音和假陽性的影響粱挡。

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

比如gene A 與 gene B，鄰接矩陣關(guān)注的是這倆基因之間直接的共表達關(guān)系屈留，而拓撲重疊矩陣還考慮中間因素的影響局冰，不止關(guān)注gene A對B的調(diào)控，還要關(guān)注geneA對B的次級影響灌危。比如gene A調(diào)控 gene C康二，gene C又調(diào)控gene B。

既然考慮了次級調(diào)控勇蝙，為什么不考慮更次一級的調(diào)控呢沫勿？這里將其稱之為三級調(diào)控。如味混，A不僅自己調(diào)控B产雹，還通過調(diào)控C，C調(diào)控D翁锡，D又調(diào)控B的方式蔓挖，間接的調(diào)控B。兩個原因馆衔，一個可能是計算量太大瘟判，太耗費計算資源了。另一個可能是角溃，三級調(diào)控對基因B的影響基本可以忽略不計了拷获，這樣也就可以解釋為什么沒有四級調(diào)控，五級調(diào)控……

那么拓撲重疊矩陣怎么計算得到呢减细？具體看如下公式：

image.png

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

image

鄰接矩陣中茧妒，A與D的關(guān)系是1。而在拓撲重疊矩陣中糠馆，A和D的關(guān)系嘶伟，不僅需要考慮A和D之間的虛線，還需要A-B-D以及A-C-D又碌。（通過B與C的兩條間接關(guān)系）

ABBC=11,=1 ACCD=11=1九昧，因此公式中l(wèi) =1+1 = 2， a =1毕匀。l + a 就是A 到D 的直接和間接共表達的總和铸鹰。而我們可以知道A的k是3 (k )，D的k也是3皂岔，因此min{k ,k }=3蹋笼，因此最終w 的結(jié)果為（2+1）/（3+1-1）=1，即我們得到的拓撲重疊矩陣中躁垛，第i行j列的元素值為1剖毯。

利用TOM進行聚類

有了拓撲異構(gòu)矩陣（TOM）之后，我們需要進行聚類教馆，從而對基因分塊逊谋，得到不同的模塊（module）。為了更好的展示聚類結(jié)果土铺，我們用1減去TOM胶滋，從而得到對應(yīng)的基因不相似性的矩陣（dissTOM）。

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

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

image

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

識別基因模塊

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

通常將模塊樹的分支切分開的方法有兩種起意，一種是指定一個高度，將結(jié)果分成若干個模塊（即：數(shù)據(jù)預(yù)處理部分的紅線）病瞳。另一種是使用 dynamic branch cut methods揽咕，該方法的優(yōu)點是可以自動化的將module區(qū)分開，不再需要手動指定高度套菜，且該方法更為靈活亲善。

使用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作為一種自動切分模塊的方法，有時可能會識別出來兩個表達譜非常相似的模塊逗柴，此時就需要我們手動將這些高度共表達的模塊合并蛹头，方便后續(xù)分析。首先需要對不同的模塊進行聚類：

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

image

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

image.png

認真看耕拷，可以觀察到，紫色和黃色兩個色塊合并到了一起肋僧，綠色和紅色兩種module合并到了一起斑胜。

模塊功能的驗證

層次聚類的一個缺點是很難確定究竟應(yīng)該將整個數(shù)據(jù)集分成幾類，盡管我們采用多種分支剪切/模塊識別的方法嫌吠，但是數(shù)據(jù)的分類問題止潘，是一個開放式的探索問題，并沒有一個固定解辫诅。即：合理即可凭戴。

那么怎么劃分模塊才是比較合理的呢？通常來說一個共表達模塊反應(yīng)著一個真實的生物學信號炕矮，比如通路么夫，比如某個生物學過程》羰樱或者反應(yīng)著噪音档痪，如組織污染，技術(shù)偏差邢滑，假陽性等腐螟。

為了判斷一個模塊的劃分是否合理，通常我們會對一個模塊進行GO富集分析，判斷是否富集在某一個或者某幾個通路乐纸。

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

在WGCNA中衬廷，我們將劃分得到的模塊，與外部性狀相關(guān)聯(lián)汽绢，來查看每個模塊分別與什么性狀高度相關(guān)吗跋，從而找到我們感興趣性狀的對應(yīng)模塊。這些性狀可以是各種指標宁昭，比如小鼠的重量跌宛，尾長，毛色等久窟；也可以是對應(yīng)植物的高矮秩冈，果實多少，顆粒飽滿度等斥扛。

但是我們知道計算顯著性入问，通常是兩個向量之間作比較。但是一個模塊是一個矩陣稀颁，矩陣和向量之間芬失，計算相關(guān)性，我們該如何分析呢匾灶？WGCNA選擇使用PCA的方法棱烂，計算模塊矩陣的主成分，并將其中的PC1定義為eigengenes阶女。

有了eigengennes 我們就可以計算模塊與外部性狀的相關(guān)性了颊糜。可視化結(jié)果如下：

image

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

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

首先我們將模塊中的每個基因與eigengenes進行相關(guān)性分析抛蚤，得到的結(jié)果我們稱之為module membership（MM）。當結(jié)果越接近于0寻狂，則我們認為該基因與其所在的模塊越不相關(guān)岁经。而結(jié)果越接近于1或者-1，則我們認為該基因與模塊基因高度相關(guān)蛇券，正負表示其是正相關(guān)還是負相關(guān)缀壤。

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

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

image

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

R中的WGCNA可視化

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

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

image

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

image

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

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

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

Cytoscape如果講起來苹熏，那將又是一堆字碟贾，還是再單獨開個教程講解吧

image

完結(jié)撒花币喧，原理/思路部分告一段落，下一部分就是WGCNA的實戰(zhàn)部分了袱耽。實戰(zhàn)部分將會跟著原理部分一步一步走下來杀餐，分別講述對應(yīng)的代碼含義和結(jié)果。如果不想等待朱巨，而且不反感英語史翘，強烈建議看看這個英文教程：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://www.reibang.com/p/fe4d3c77786e
• http://www.sci666.net/58731.html
• [http://www.biotrainee.com/thread-646-1-1.html](