原文地址：https://mp.weixin.qq.com/s/BkM9vtX_gbehUFmuXfshUw

WGCNA基本概念

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析 (WGCNA, Weighted correlation network analysis)是用來描述不同樣品之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法帆喇，可以用來鑒定高度協(xié)同變化的基因集, 并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定候補(bǔ)生物標(biāo)記基因或治療靶點(diǎn)。

相比于只關(guān)注差異表達(dá)的基因乌妒，WGCNA利用數(shù)千或近萬個(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侦啸，需要先理解下面幾個(gè)術(shù)語和它們?cè)赪GCNA中的定義。

共表達(dá)網(wǎng)絡(luò)：定義為加權(quán)基因網(wǎng)絡(luò)丧枪。點(diǎn)代表基因光涂，邊代表基因表達(dá)相關(guān)性。加權(quán)是指對(duì)相關(guān)性值進(jìn)行冥次運(yùn)算(冥次的值也就是軟閾值?(power, pickSoftThreshold這個(gè)函數(shù)所做的就是確定合適的power))豪诲。無向網(wǎng)絡(luò)的邊屬性計(jì)算方式為abs(cor(genex, geney)) ^ power顶捷；有向網(wǎng)絡(luò)的邊屬性計(jì)算方式為(1+cor(genex, geney)/2) ^ power; sign hybrid的邊屬性計(jì)算方式為cor(genex, geney)^power if cor>0 else 0。這種處理方式強(qiáng)化了強(qiáng)相關(guān)屎篱，弱化了弱相關(guān)或負(fù)相關(guān)服赎，使得相關(guān)性數(shù)值更符合無標(biāo)度網(wǎng)絡(luò)特征葵蒂，更具有生物意義。如果沒有合適的power重虑，一般是由于部分樣品與其它樣品因?yàn)槟撤N原因差別太大導(dǎo)致的践付，可根據(jù)具體問題移除部分樣品或查看后面的經(jīng)驗(yàn)值。

Module(模塊)：高度內(nèi)連的基因集缺厉。在無向網(wǎng)絡(luò)中永高，模塊內(nèi)是高度相關(guān)的基因。在有向網(wǎng)絡(luò)中提针，模塊內(nèi)是高度正相關(guān)的基因命爬。把基因聚類成模塊后，可以對(duì)每個(gè)模塊進(jìn)行三個(gè)層次的分析：1. 功能富集分析查看其功能特征是否與研究目的相符辐脖；2. 模塊與性狀進(jìn)行關(guān)聯(lián)分析饲宛，找出與關(guān)注性狀相關(guān)度最高的模塊；3. 模塊與樣本進(jìn)行關(guān)聯(lián)分析嗜价，找到樣品特異高表達(dá)的模塊艇抠。

基因富集相關(guān)文章?去東方，最好用的在線GO富集分析工具久锥；GO家淤、GSEA富集分析一網(wǎng)打進(jìn)；GSEA富集分析-界面操作瑟由。其它關(guān)聯(lián)后面都會(huì)提及絮重。

Connectivity (連接度)：類似于網(wǎng)絡(luò)中 “度” (degree)的概念。每個(gè)基因的連接度是與其相連的基因的邊屬性之和错妖。

Module eigengene E: 給定模型的第一主成分绿鸣，代表整個(gè)模型的基因表達(dá)譜。這個(gè)是個(gè)很巧妙的梳理暂氯，我們之前講過PCA分析的降維作用，之前主要是拿來做可視化亮蛔，現(xiàn)在用到這個(gè)地方痴施，很好的用一個(gè)向量代替了一個(gè)矩陣，方便后期計(jì)算究流。(降維除了PCA辣吃，還可以看看tSNE)

Intramodular connectivity: 給定基因與給定模型內(nèi)其他基因的關(guān)聯(lián)度，判斷基因所屬關(guān)系芬探。

Module membership: 給定基因表達(dá)譜與給定模型的eigengene的相關(guān)性神得。

Hub gene: 關(guān)鍵基因 (連接度最多或連接多個(gè)模塊的基因)。

Adjacency matrix (鄰接矩陣)：基因和基因之間的加權(quán)相關(guān)性值構(gòu)成的矩陣偷仿。

TOM (Topological overlap matrix)：把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣哩簿，以降低噪音和假相關(guān)宵蕉，獲得的新距離矩陣，這個(gè)信息可拿來構(gòu)建網(wǎng)絡(luò)或繪制TOM圖节榜。

基本分析流程

構(gòu)建基因共表達(dá)網(wǎng)絡(luò)：使用加權(quán)的表達(dá)相關(guān)性羡玛。

識(shí)別基因集：基于加權(quán)相關(guān)性，進(jìn)行層級(jí)聚類分析宗苍，并根據(jù)設(shè)定標(biāo)準(zhǔn)切分聚類結(jié)果稼稿，獲得不同的基因模塊，用聚類樹的分枝和不同顏色表示讳窟。

如果有表型信息让歼，計(jì)算基因模塊與表型的相關(guān)性，鑒定性狀相關(guān)的模塊丽啡。

研究模型之間的關(guān)系是越，從系統(tǒng)層面查看不同模型的互作網(wǎng)絡(luò)。

從關(guān)鍵模型中選擇感興趣的驅(qū)動(dòng)基因碌上，或根據(jù)模型中已知基因的功能推測(cè)未知基因的功能倚评。

導(dǎo)出TOM矩陣，繪制相關(guān)性圖馏予。

WGCNA包實(shí)戰(zhàn)

R包WGCNA是用于計(jì)算各種加權(quán)關(guān)聯(lián)分析的功能集合天梧，可用于網(wǎng)絡(luò)構(gòu)建，基因篩選霞丧，基因簇鑒定呢岗，拓?fù)涮卣饔?jì)算，數(shù)據(jù)模擬和可視化等蛹尝。

輸入數(shù)據(jù)和參數(shù)選擇

WGCNA本質(zhì)是基于相關(guān)系數(shù)的網(wǎng)絡(luò)分析方法后豫，適用于多樣品數(shù)據(jù)模式，一般要求樣本數(shù)多于15個(gè)突那。樣本數(shù)多于20時(shí)效果更好挫酿，樣本越多，結(jié)果越穩(wěn)定愕难。

基因表達(dá)矩陣: 常規(guī)表達(dá)矩陣即可早龟，即基因在行，樣品在列猫缭，進(jìn)入分析前做一個(gè)轉(zhuǎn)置葱弟。RPKM、FPKM或其它標(biāo)準(zhǔn)化方法影響不大猜丹，推薦使用Deseq2的varianceStabilizingTransformation或log2(x+1)對(duì)標(biāo)準(zhǔn)化后的數(shù)據(jù)做個(gè)轉(zhuǎn)換芝加。如果數(shù)據(jù)來自不同的批次，需要先移除批次效應(yīng) (記得上次轉(zhuǎn)錄組培訓(xùn)課講過如何操作)射窒。如果數(shù)據(jù)存在系統(tǒng)偏移藏杖，需要做下quantile normalization将塑。

性狀矩陣：用于關(guān)聯(lián)分析的性狀必須是數(shù)值型特征 (如下面示例中的Height,?Weight,Diameter)。如果是區(qū)域或分類變量制市，需要轉(zhuǎn)換為0-1矩陣的形式(1表示屬于此組或有此屬性抬旺，0表示不屬于此組或無此屬性，如樣品分組信息WT, KO, OE)祥楣。

ID ?WT ?KO ?OE Height Weight Diameter

samp1 ? 1 ? 0 ? 0 ? 1 ? 2 ? 3

samp2 ? 1 ? 0 ? 0 ? 2 ? 4 ? 6

samp3 ? 0 ? 1 ? 0 ? 10 ?20 ?50

samp4 ? 0 ? 1 ? 0 ? 15 ?30 ?80

samp5 ? 0 ? 0 ? 1 ? NA ?9 ? 8

samp6 ? 0 ? 0 ? 1 ? 4 ? 8 ? 7

推薦使用Signed network和Robust correlation (bicor)开财。(這個(gè)根據(jù)自己的需要，看看上面寫的每個(gè)網(wǎng)絡(luò)怎么計(jì)算的误褪，更知道怎么選擇)

無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)责鳍，沒有一個(gè)power值可以使無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8或平均連接度降到100以下，可能是由于部分樣品與其他樣品差別太大造成的兽间。這可能由批次效應(yīng)历葛、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)表達(dá)影響太大等造成, 可以通過繪制樣品聚類查看分組信息、關(guān)聯(lián)批次信息嘀略、處理信息和有無異常樣品 (可以使用之前講過的熱圖簡化恤溶，增加行或列屬性)。如果這確實(shí)是由有意義的生物變化引起的帜羊，也可以使用后面程序中的經(jīng)驗(yàn)power值咒程。

安裝WGCNA

WGCNA依賴的包比較多，bioconductor上的包需要自己安裝讼育，cran上依賴的包可以自動(dòng)安裝帐姻。通常在R中運(yùn)行下面4條語句就可以完成WGCNA的安裝。

建議在編譯安裝R時(shí)增加--with-blas --with-lapack提高矩陣運(yùn)算的速度奶段，具體見R和Rstudio安裝饥瓷。

#source("https://bioconductor.org/biocLite.R")

#biocLite(c("AnnotationDbi", "impute","GO.db", "preprocessCore"))

#site="https://mirrors.tuna.tsinghua.edu.cn/CRAN"

#install.packages(c("WGCNA", "stringr", "reshape2"), repos=site)

WGCNA實(shí)戰(zhàn)

實(shí)戰(zhàn)采用的是官方提供的清理后的矩陣，原矩陣信息太多痹籍，容易給人誤導(dǎo)呢铆，后臺(tái)回復(fù)WGCNA?獲取數(shù)據(jù)。

數(shù)據(jù)讀入

library(WGCNA)

## Loading required package: dynamicTreeCut

## Loading required package: fastcluster

##

## Attaching package: 'fastcluster'

## The following object is masked from 'package:stats':

##

## ? ? hclust

## ==========================================================================

## *

## * ?Package WGCNA 1.63 loaded.

## *

## * ? ?Important note: It appears that your system supports multi-threading,

## * ? ?but it is not enabled within WGCNA in R.

## * ? ?To allow multi-threading within WGCNA with all available cores, use

## *

## * ? ? ? ? ?allowWGCNAThreads()

## *

## * ? ?within R. Use disableWGCNAThreads() to disable threading if necessary.

## * ? ?Alternatively, set the following environment variable on your system:

## *

## * ? ? ? ? ?ALLOW_WGCNA_THREADS=

## *

## * ? ?for example

## *

## * ? ? ? ? ?ALLOW_WGCNA_THREADS=48

## *

## * ? ?To set the environment variable in linux bash shell, type

## *

## * ? ? ? ? ? export ALLOW_WGCNA_THREADS=48

## *

## * ? ? before running R. Other operating systems or shells will

## * ? ? have a similar command to achieve the same aim.

## *

## ==========================================================================

##

## Attaching package: 'WGCNA'

## The following object is masked from 'package:stats':

##

## ? ? cor

library(reshape2)

library(stringr)

#

options(stringsAsFactors = FALSE)

# 打開多線程

enableWGCNAThreads()

## Allowing parallel execution with up to 47 working processes.

# 常規(guī)表達(dá)矩陣词裤，log2轉(zhuǎn)換后或

# Deseq2的varianceStabilizingTransformation轉(zhuǎn)換的數(shù)據(jù)

# 如果有批次效應(yīng)刺洒，需要事先移除，可使用removeBatchEffect

# 如果有系統(tǒng)偏移(可用boxplot查看基因表達(dá)分布是否一致)吼砂，

# 需要quantile normalization

exprMat <- "WGCNA/LiverFemaleClean.txt"

# 官方推薦 "signed" 或 "signed hybrid"

# 為與原文檔一致，故未修改

type = "unsigned"

# 相關(guān)性計(jì)算

# 官方推薦 biweight mid-correlation & bicor

# corType: pearson or bicor

# 為與原文檔一致鼎文，故未修改

corType = "pearson"

corFnc = ifelse(corType=="pearson", cor, bicor)

# 對(duì)二元變量渔肩，如樣本性狀信息計(jì)算相關(guān)性時(shí)，

# 或基因表達(dá)嚴(yán)重依賴于疾病狀態(tài)時(shí)拇惋，需設(shè)置下面參數(shù)

maxPOutliers = ifelse(corType=="pearson",1,0.05)

# 關(guān)聯(lián)樣品性狀的二元變量時(shí)周偎，設(shè)置

robustY = ifelse(corType=="pearson",T,F)

##導(dǎo)入數(shù)據(jù)##

dataExpr <- read.table(exprMat, sep='\t', row.names=1, header=T,

? ? ? ? ? ? ? ? ? ? quote="", comment="", check.names=F)

dim(dataExpr)

## [1] 3600 ?134

head(dataExpr)[,1:8]

## ? ? ? ? ? ? ? ? F2_2 ? ?F2_3 ? ? F2_14 ? ?F2_15 ? ?F2_19 ? ? ? F2_20

## MMT00000044 -0.01810 ?0.0642 ?6.44e-05 -0.05800 ?0.04830 -0.15197410

## MMT00000046 -0.07730 -0.0297 ?1.12e-01 -0.05890 ?0.04430 -0.09380000

## MMT00000051 -0.02260 ?0.0617 -1.29e-01 ?0.08710 -0.11500 -0.06502607

## MMT00000076 -0.00924 -0.1450 ?2.87e-02 -0.04390 ?0.00425 -0.23610000

## MMT00000080 -0.04870 ?0.0582 -4.83e-02 -0.03710 ?0.02510 ?0.08504274

## MMT00000102 ?0.17600 -0.1890 -6.50e-02 -0.00846 -0.00574 -0.01807182

## ? ? ? ? ? ? ? ?F2_23 ? ?F2_24

## MMT00000044 -0.00129 -0.23600

## MMT00000046 ?0.09340 ?0.02690

## MMT00000051 ?0.00249 -0.10200

## MMT00000076 -0.06900 ?0.01440

## MMT00000080 ?0.04450 ?0.00167

## MMT00000102 -0.12500 -0.06820

數(shù)據(jù)篩選

## 篩選中位絕對(duì)偏差前75%的基因抹剩，至少M(fèi)AD大于0.01

## 篩選后會(huì)降低運(yùn)算量，也會(huì)失去部分信息

## 也可不做篩選蓉坎，使MAD大于0即可

m.mad <- apply(dataExpr,1,mad)

dataExprVar <- dataExpr[which(m.mad >

? ? ? ? ? ? ? ? max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.01)),]

## 轉(zhuǎn)換為樣品在行澳眷，基因在列的矩陣

dataExpr <- as.data.frame(t(dataExprVar))

## 檢測(cè)缺失值

gsg = goodSamplesGenes(dataExpr, verbose = 3)

## ?Flagging genes and samples with too many missing values...

## ? ..step 1

if (!gsg$allOK){

?# Optionally, print the gene and sample names that were removed:

?if (sum(!gsg$goodGenes)>0)

? ?printFlush(paste("Removing genes:",

? ? ? ? ? ? ? ? ? ? paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));

?if (sum(!gsg$goodSamples)>0)

? ?printFlush(paste("Removing samples:",

? ? ? ? ? ? ? ? ? ? paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));

?# Remove the offending genes and samples from the data:

?dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]

}

nGenes = ncol(dataExpr)

nSamples = nrow(dataExpr)

dim(dataExpr)

## [1] ?134 2697

head(dataExpr)[,1:8]

## ? ? ? MMT00000051 MMT00000080 MMT00000102 MMT00000149 MMT00000159

## F2_2 ?-0.02260000 -0.04870000 ?0.17600000 ?0.07680000 -0.14800000

## F2_3 ? 0.06170000 ?0.05820000 -0.18900000 ?0.18600000 ?0.17700000

## F2_14 -0.12900000 -0.04830000 -0.06500000 ?0.21400000 -0.13200000

## F2_15 ?0.08710000 -0.03710000 -0.00846000 ?0.12000000 ?0.10700000

## F2_19 -0.11500000 ?0.02510000 -0.00574000 ?0.02100000 -0.11900000

## F2_20 -0.06502607 ?0.08504274 -0.01807182 ?0.06222751 -0.05497686

## ? ? ? MMT00000207 MMT00000212 MMT00000241

## F2_2 ? 0.06870000 ?0.06090000 -0.01770000

## F2_3 ? 0.10100000 ?0.05570000 -0.03690000

## F2_14 ?0.10900000 ?0.19100000 -0.15700000

## F2_15 -0.00858000 -0.12100000 ?0.06290000

## F2_19 ?0.10500000 ?0.05410000 -0.17300000

## F2_20 -0.02441415 ?0.06343181 ?0.06627665

軟閾值篩選

## 查看是否有離群樣品

sampleTree = hclust(dist(dataExpr), method = "average")

plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")

軟閾值的篩選原則是使構(gòu)建的網(wǎng)絡(luò)更符合無標(biāo)度網(wǎng)絡(luò)特征。

powers = c(c(1:10), seq(from = 12, to=30, by=2))

sft = pickSoftThreshold(dataExpr, powerVector=powers,

? ? ? ? ? ? ? ? ? ? ? ?networkType=type, verbose=5)

## pickSoftThreshold: will use block size 2697.

## ?pickSoftThreshold: calculating connectivity for given powers...

## ? ?..working on genes 1 through 2697 of 2697

## ? ?Power SFT.R.sq ?slope truncated.R.sq mean.k. median.k. max.k.

## 1 ? ? ?1 ? 0.1370 ?0.825 ? ? ? ? ?0.412 587.000 ?5.95e+02 ?922.0

## 2 ? ? ?2 ? 0.0416 -0.332 ? ? ? ? ?0.630 206.000 ?2.02e+02 ?443.0

## 3 ? ? ?3 ? 0.2280 -0.747 ? ? ? ? ?0.920 ?91.500 ?8.43e+01 ?247.0

## 4 ? ? ?4 ? 0.3910 -1.120 ? ? ? ? ?0.908 ?47.400 ?4.02e+01 ?154.0

## 5 ? ? ?5 ? 0.7320 -1.230 ? ? ? ? ?0.958 ?27.400 ?2.14e+01 ?102.0

## 6 ? ? ?6 ? 0.8810 -1.490 ? ? ? ? ?0.916 ?17.200 ?1.22e+01 ? 83.7

## 7 ? ? ?7 ? 0.8940 -1.640 ? ? ? ? ?0.869 ?11.600 ?7.29e+00 ? 75.4

## 8 ? ? ?8 ? 0.8620 -1.660 ? ? ? ? ?0.827 ? 8.250 ?4.56e+00 ? 69.2

## 9 ? ? ?9 ? 0.8200 -1.600 ? ? ? ? ?0.810 ? 6.160 ?2.97e+00 ? 64.2

## 10 ? ?10 ? 0.8390 -1.560 ? ? ? ? ?0.855 ? 4.780 ?2.01e+00 ? 60.1

## 11 ? ?12 ? 0.8020 -1.410 ? ? ? ? ?0.866 ? 3.160 ?9.61e-01 ? 53.2

## 12 ? ?14 ? 0.8470 -1.340 ? ? ? ? ?0.909 ? 2.280 ?4.84e-01 ? 47.7

## 13 ? ?16 ? 0.8850 -1.250 ? ? ? ? ?0.932 ? 1.750 ?2.64e-01 ? 43.1

## 14 ? ?18 ? 0.8830 -1.210 ? ? ? ? ?0.922 ? 1.400 ?1.46e-01 ? 39.1

## 15 ? ?20 ? 0.9110 -1.180 ? ? ? ? ?0.926 ? 1.150 ?8.35e-02 ? 35.6

## 16 ? ?22 ? 0.9160 -1.140 ? ? ? ? ?0.927 ? 0.968 ?5.02e-02 ? 32.6

## 17 ? ?24 ? 0.9520 -1.120 ? ? ? ? ?0.961 ? 0.828 ?2.89e-02 ? 29.9

## 18 ? ?26 ? 0.9520 -1.120 ? ? ? ? ?0.944 ? 0.716 ?1.77e-02 ? 27.5

## 19 ? ?28 ? 0.9380 -1.120 ? ? ? ? ?0.922 ? 0.626 ?1.08e-02 ? 25.4

## 20 ? ?30 ? 0.9620 -1.110 ? ? ? ? ?0.951 ? 0.551 ?6.49e-03 ? 23.5

par(mfrow = c(1,2))

cex1 = 0.9

# 橫軸是Soft threshold (power)蛉艾，縱軸是無標(biāo)度網(wǎng)絡(luò)的評(píng)估參數(shù)钳踊，數(shù)值越高，

# 網(wǎng)絡(luò)越符合無標(biāo)度特征 (non-scale)

plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],

? ? xlab="Soft Threshold (power)",

? ? ylab="Scale Free Topology Model Fit,signed R^2",type="n",

? ? main = paste("Scale independence"))

text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],

? ? labels=powers,cex=cex1,col="red")

# 篩選標(biāo)準(zhǔn)勿侯。R-square=0.85

abline(h=0.85,col="red")

# Soft threshold與平均連通性

plot(sft$fitIndices[,1], sft$fitIndices[,5],

? ? xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",

? ? main = paste("Mean connectivity"))

text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers,

? ? cex=cex1, col="red")

power = sft$powerEstimate

power

## [1] 6

經(jīng)驗(yàn)power?(無滿足條件的power時(shí)選用)

# 無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)拓瞪，沒有一個(gè)power值可以使

# 無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8，平均連接度較高如在100以上助琐，可能是由于

# 部分樣品與其他樣品差別太大祭埂。這可能由批次效應(yīng)、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)

# 表達(dá)影響太大等造成兵钮∏穑可以通過繪制樣品聚類查看分組信息和有無異常樣品。

# 如果這確實(shí)是由有意義的生物變化引起的掘譬，也可以使用下面的經(jīng)驗(yàn)power值泰演。

if (is.na(power)){

?power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),

? ? ? ? ?ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),

? ? ? ? ?ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),

? ? ? ? ?ifelse(type == "unsigned", 6, 12)) ? ? ?

? ? ? ? ?)

? ? ? ? ?)

}

網(wǎng)絡(luò)構(gòu)建

##一步法網(wǎng)絡(luò)構(gòu)建：One-step network construction and module detection##

# power: 上一步計(jì)算的軟閾值

# maxBlockSize: 計(jì)算機(jī)能處理的最大模塊的基因數(shù)量 (默認(rèn)5000)；

# ?4G內(nèi)存電腦可處理8000-10000個(gè)屁药，16G內(nèi)存電腦可以處理2萬個(gè)粥血，32G內(nèi)存電腦可

# ?以處理3萬個(gè)

# ?計(jì)算資源允許的情況下最好放在一個(gè)block里面。

# corType: pearson or bicor

# numericLabels: 返回?cái)?shù)字而不是顏色作為模塊的名字酿箭，后面可以再轉(zhuǎn)換為顏色

# saveTOMs：最耗費(fèi)時(shí)間的計(jì)算复亏，存儲(chǔ)起來，供后續(xù)使用

# mergeCutHeight: 合并模塊的閾值缭嫡，越大模塊越少

net = blockwiseModules(dataExpr, power = power, maxBlockSize = nGenes,

? ? ? ? ? ? ? ? ? ? ? TOMType = type, minModuleSize = 30,

? ? ? ? ? ? ? ? ? ? ? reassignThreshold = 0, mergeCutHeight = 0.25,

? ? ? ? ? ? ? ? ? ? ? numericLabels = TRUE, pamRespectsDendro = FALSE,

? ? ? ? ? ? ? ? ? ? ? saveTOMs=TRUE, corType = corType,

? ? ? ? ? ? ? ? ? ? ? maxPOutliers=maxPOutliers, loadTOM?=TRUE,

? ? ? ? ? ? ? ? ? ? ? saveTOMFileBase = paste0(exprMat, ".tom"),

? ? ? ? ? ? ? ? ? ? ? verbose = 3)

## ?Calculating module eigengenes block-wise from all genes

## ? ?Flagging genes and samples with too many missing values...

## ? ? ..step 1

## ?..Working on block 1 .

## ? ? TOM calculation: adjacency..

## ? ? ..will use 47 parallel threads.

## ? ? ?Fraction of slow calculations: 0.000000

## ? ? ..connectivity..

## ? ? ..matrix multiplication (system BLAS)..

## ? ? ..normalization..

## ? ? ..done.

## ? ?..saving TOM for block 1 into file WGCNA/LiverFemaleClean.txt.tom-block.1.RData

## ?....clustering..

## ?....detecting modules..

## ?....calculating module eigengenes..

## ?....checking kME in modules..

## ? ? ?..removing 3 genes from module 1 because their KME is too low.

## ? ? ?..removing 5 genes from module 12 because their KME is too low.

## ? ? ?..removing 1 genes from module 14 because their KME is too low.

## ?..merging modules that are too close..

## ? ? ?mergeCloseModules: Merging modules whose distance is less than 0.25

## ? ? ? ?Calculating new MEs...

# 根據(jù)模塊中基因數(shù)目的多少缔御，降序排列，依次編號(hào)為 `1-最大模塊數(shù)`妇蛀。

# **0 (grey)**表示**未**分入任何模塊的基因耕突。

table(net$colors)

##

## ? 0 ? 1 ? 2 ? 3 ? 4 ? 5 ? 6 ? 7 ? 8 ? 9 ?10 ?11 ?12 ?13

## 135 472 356 333 307 303 177 158 102 ?94 ?69 ?66 ?63 ?62

層級(jí)聚類樹展示各個(gè)模塊

## 灰色的為**未分類**到模塊的基因。

# Convert labels to colors for plotting

moduleLabels = net$colors

moduleColors = labels2colors(moduleLabels)

# Plot the dendrogram and the module colors underneath

# 如果對(duì)結(jié)果不滿意评架，還可以recutBlockwiseTrees眷茁，節(jié)省計(jì)算時(shí)間

plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],

? ? ? ? ? ? ? ? ? ?"Module colors",

? ? ? ? ? ? ? ? ? ?dendroLabels = FALSE, hang = 0.03,

? ? ? ? ? ? ? ? ? ?addGuide = TRUE, guideHang = 0.05)

繪制模塊之間相關(guān)性圖

# module eigengene, 可以繪制線圖，作為每個(gè)模塊的基因表達(dá)趨勢(shì)的展示

MEs = net$MEs

### 不需要重新計(jì)算纵诞，改下列名字就好

### 官方教程是重新計(jì)算的上祈，起始可以不用這么麻煩

MEs_col = MEs

colnames(MEs_col) = paste0("ME", labels2colors(

?as.numeric(str_replace_all(colnames(MEs),"ME",""))))

MEs_col = orderMEs(MEs_col)

# 根據(jù)基因間表達(dá)量進(jìn)行聚類所得到的各模塊間的相關(guān)性圖

# marDendro/marHeatmap 設(shè)置下、左、上登刺、右的邊距

plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap",

? ? ? ? ? ? ? ? ? ? ?marDendro = c(3,3,2,4),

? ? ? ? ? ? ? ? ? ? ?marHeatmap = c(3,4,2,2), plotDendrograms = T,

? ? ? ? ? ? ? ? ? ? ?xLabelsAngle = 90)

## 如果有表型數(shù)據(jù)籽腕，也可以跟ME數(shù)據(jù)放一起，一起出圖

#MEs_colpheno = orderMEs(cbind(MEs_col, traitData))

#plotEigengeneNetworks(MEs_colpheno, "Eigengene adjacency heatmap",

# ? ? ? ? ? ? ? ? ? ? ?marDendro = c(3,3,2,4),

# ? ? ? ? ? ? ? ? ? ? ?marHeatmap = c(3,4,2,2), plotDendrograms = T,

# ? ? ? ? ? ? ? ? ? ? ?xLabelsAngle = 90)

可視化基因網(wǎng)絡(luò) (TOM plot)

# 如果采用分步計(jì)算纸俭，或設(shè)置的blocksize>=總基因數(shù)皇耗，直接load計(jì)算好的TOM結(jié)果

# 否則需要再計(jì)算一遍，比較耗費(fèi)時(shí)間

# TOM = TOMsimilarityFromExpr(dataExpr, power=power, corType=corType, networkType=type)

load(net$TOMFiles[1], verbose=T)

## Loading objects:

## ? TOM

TOM <- as.matrix(TOM)

dissTOM = 1-TOM

# Transform dissTOM with a power to make moderately strong

# connections more visible in the heatmap

plotTOM = dissTOM^7

# Set diagonal to NA for a nicer plot

diag(plotTOM) = NA

# Call the plot function

# 這一部分特別耗時(shí)揍很，行列同時(shí)做層級(jí)聚類

TOMplot(plotTOM, net$dendrograms, moduleColors,

? ? ? ?main = "Network heatmap plot, all genes")

導(dǎo)出網(wǎng)絡(luò)用于Cytoscape

Cytoscape繪制網(wǎng)絡(luò)圖見我們更新版的視頻教程或https://bioinfo.ke.qq.com/郎楼。

probes = colnames(dataExpr)

dimnames(TOM) <- list(probes, probes)

# Export the network into edge and node list files Cytoscape can read

# threshold 默認(rèn)為0.5, 可以根據(jù)自己的需要調(diào)整，也可以都導(dǎo)出后在

# cytoscape中再調(diào)整

cyt = exportNetworkToCytoscape(TOM,

? ? ? ? ? ? edgeFile = paste(exprMat, ".edges.txt", sep=""),

? ? ? ? ? ? nodeFile = paste(exprMat, ".nodes.txt", sep=""),

? ? ? ? ? ? weighted = TRUE, threshold = 0,

? ? ? ? ? ? nodeNames = probes, nodeAttr = moduleColors)

關(guān)聯(lián)表型數(shù)據(jù)

trait <- "WGCNA/TraitsClean.txt"

# 讀入表型數(shù)據(jù)女轿，不是必須的

if(trait != "") {

?traitData <- read.table(file=trait, sep='\t', header=T, row.names=1,

? ? ? ? ? ? ? ? ? ? ? ? ?check.names=FALSE, comment='',quote="")

?sampleName = rownames(dataExpr)

?traitData = traitData[match(sampleName, rownames(traitData)), ]

}

### 模塊與表型數(shù)據(jù)關(guān)聯(lián)

if (corType=="pearsoon") {

?modTraitCor = cor(MEs_col, traitData, use = "p")

?modTraitP = corPvalueStudent(modTraitCor, nSamples)

} else {

?modTraitCorP = bicorAndPvalue(MEs_col, traitData, robustY=robustY)

?modTraitCor = modTraitCorP$bicor

?modTraitP ? = modTraitCorP$p

}

## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.

## Pearson correlation was used for individual columns with zero (or missing)

## MAD.

# signif表示保留幾位小數(shù)

textMatrix = paste(signif(modTraitCor, 2), "\n(", signif(modTraitP, 1), ")", sep = "")

dim(textMatrix) = dim(modTraitCor)

labeledHeatmap(Matrix = modTraitCor, xLabels = colnames(traitData),

? ? ? ? ? ? ? yLabels = colnames(MEs_col),

? ? ? ? ? ? ? cex.lab = 0.5,

? ? ? ? ? ? ? ySymbols = colnames(MEs_col), colorLabels = FALSE,

? ? ? ? ? ? ? colors = blueWhiteRed(50),

? ? ? ? ? ? ? textMatrix = textMatrix, setStdMargins = FALSE,

? ? ? ? ? ? ? cex.text = 0.5, zlim = c(-1,1),

? ? ? ? ? ? ? main = paste("Module-trait relationships"))

模塊內(nèi)基因與表型數(shù)據(jù)關(guān)聯(lián), 從上圖可以看到MEmagenta與Insulin_ug_l相關(guān)箭启，選取這兩部分進(jìn)行分析。

## 從上圖可以看到MEmagenta與Insulin_ug_l相關(guān)

## 模塊內(nèi)基因與表型數(shù)據(jù)關(guān)聯(lián)

# 性狀跟模塊雖然求出了相關(guān)性蛉迹，可以挑選最相關(guān)的那些模塊來分析傅寡，

# 但是模塊本身仍然包含非常多的基因，還需進(jìn)一步的尋找最重要的基因北救。

# 所有的模塊都可以跟基因算出相關(guān)系數(shù)荐操，所有的連續(xù)型性狀也可以跟基因的表達(dá)

# 值算出相關(guān)系數(shù)。

# 如果跟性狀顯著相關(guān)基因也跟某個(gè)模塊顯著相關(guān)珍策，那么這些基因可能就非常重要

# 托启。

### 計(jì)算模塊與基因的相關(guān)性矩陣

if (corType=="pearsoon") {

?geneModuleMembership = as.data.frame(cor(dataExpr, MEs_col, use = "p"))

?MMPvalue = as.data.frame(corPvalueStudent(

? ? ? ? ? ? as.matrix(geneModuleMembership), nSamples))

} else {

?geneModuleMembershipA = bicorAndPvalue(dataExpr, MEs_col, robustY=robustY)

?geneModuleMembership = geneModuleMembershipA$bicor

?MMPvalue ? = geneModuleMembershipA$p

}

# 計(jì)算性狀與基因的相關(guān)性矩陣

## 只有連續(xù)型性狀才能進(jìn)行計(jì)算写隶，如果是離散變量竿奏，在構(gòu)建樣品表時(shí)就轉(zhuǎn)為0-1矩陣了牛。

if (corType=="pearsoon") {

?geneTraitCor = as.data.frame(cor(dataExpr, traitData, use = "p"))

?geneTraitP = as.data.frame(corPvalueStudent(

? ? ? ? ? ? as.matrix(geneTraitCor), nSamples))

} else {

?geneTraitCorA = bicorAndPvalue(dataExpr, traitData, robustY=robustY)

?geneTraitCor = as.data.frame(geneTraitCorA$bicor)

?geneTraitP ? = as.data.frame(geneTraitCorA$p)

}

## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.

## Pearson correlation was used for individual columns with zero (or missing)

## MAD.

# 最后把兩個(gè)相關(guān)性矩陣聯(lián)合起來,指定感興趣模塊進(jìn)行分析

module = "magenta"

pheno = "Insulin_ug_l"

modNames = substring(colnames(MEs_col), 3)

# 獲取關(guān)注的列

module_column = match(module, modNames)

pheno_column = match(pheno,colnames(traitData))

# 獲取模塊內(nèi)的基因

moduleGenes = moduleColors == module

sizeGrWindow(7, 7)

par(mfrow = c(1,1))

# 與性狀高度相關(guān)的基因屎飘，也是與性狀相關(guān)的模型的關(guān)鍵基因

verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),

? ? ? ? ? ? ? ? ? abs(geneTraitCor[moduleGenes, pheno_column]),

? ? ? ? ? ? ? ? ? xlab = paste("Module Membership in", module, "module"),

? ? ? ? ? ? ? ? ? ylab = paste("Gene significance for", pheno),

? ? ? ? ? ? ? ? ? main = paste("Module membership vs. gene significance\n"),

? ? ? ? ? ? ? ? ? cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

分步法展示每一步都做了什么

### 計(jì)算鄰接矩陣

adjacency = adjacency(dataExpr, power = power)

### 把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣，以降低噪音和假相關(guān)厌处，獲得距離矩陣烙懦。

TOM = TOMsimilarity(adjacency)

dissTOM = 1-TOM

### 層級(jí)聚類計(jì)算基因之間的距離樹

geneTree = hclust(as.dist(dissTOM), method = "average")

### 模塊合并

# We like large modules, so we set the minimum module size relatively high:

minModuleSize = 30

# Module identification using dynamic tree cut:

dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,

? ? ? ? ? ? ? ? ? ? ? ? ? ?deepSplit = 2, pamRespectsDendro = FALSE,

? ? ? ? ? ? ? ? ? ? ? ? ? ?minClusterSize = minModuleSize)

# Convert numeric lables into colors

dynamicColors = labels2colors(dynamicMods)

### 通過計(jì)算模塊的代表性模式和模塊之間的定量相似性評(píng)估器躏，合并表達(dá)圖譜相似

的模塊

MEList = moduleEigengenes(datExpr, colors = dynamicColors)

MEs = MEList$eigengenes

# Calculate dissimilarity of module eigengenes

MEDiss = 1-cor(MEs)

# Cluster module eigengenes

METree = hclust(as.dist(MEDiss), method = "average")

MEDissThres = 0.25

# Call an automatic merging function

merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)

# The merged module colors

mergedColors = merge$colors;

# Eigengenes of the new merged

## 分步法完結(jié)