WGCNA分析，簡單全面的最新教程（可以在線做了）

生信學習的正確姿勢（第三版）

NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫(yī)學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

高顏值免費在線繪圖工具升級版來了~~~ （可以在線做 WGCNA 分析

WGCNA基本概念

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

相比于只關注差異表達的基因，WGCNA利用數千或近萬個變化最大的基因或全部基因的信息識別感興趣的基因集，并與表型進行顯著性關聯(lián)分析。一是充分利用了信息，二是把數千個基因與表型的關聯(lián)轉換為數個基因集與表型的關聯(lián)，免去了多重假設檢驗校正的問題。

理解WGCNA，需要先理解下面幾個術語和它們在WGCNA中的定義。

• 共表達網絡：定義為加權基因網絡。點代表基因，邊代表基因表達相關性。加權是指對相關性值進行冥次運算(冥次的值也就是軟閾值?(power, pickSoftThreshold這個函數所做的就是確定合適的power))。無向網絡的邊屬性計算方式為abs(cor(genex, geney)) ^ power；有向網絡的邊屬性計算方式為(1+cor(genex, geney)/2) ^ power; sign hybrid的邊屬性計算方式為cor(genex, geney)^power if cor>0 else 0。這種處理方式強化了強相關，弱化了弱相關或負相關，使得相關性數值更符合無標度網絡特征，更具有生物意義。如果沒有合適的power，一般是由于部分樣品與其它樣品因為某種原因差別太大導致的，可根據具體問題移除部分樣品或查看后面的經驗值。

• Module(模塊)：高度內連的基因集。在無向網絡中，模塊內是高度相關的基因。在有向網絡中，模塊內是高度正相關的基因。把基因聚類成模塊后，可以對每個模塊進行三個層次的分析：1. 功能富集分析查看其功能特征是否與研究目的相符；2. 模塊與性狀進行關聯(lián)分析，找出與關注性狀相關度最高的模塊；3. 模塊與樣本進行關聯(lián)分析，找到樣品特異高表達的模塊。

  基因富集相關文章?去東方，最好用的在線GO富集分析工具；GO、GSEA富集分析一網打進；GSEA富集分析-界面操作。其它關聯(lián)后面都會提及。

• Connectivity (連接度)：類似于網絡中 “度” (degree)的概念。每個基因的連接度是與其相連的基因的邊屬性之和。

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

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

• Module membership: 給定基因表達譜與給定模型的eigengene的相關性。

• Hub gene: 關鍵基因 (連接度最多或連接多個模塊的基因)。

• Adjacency matrix (鄰接矩陣)：基因和基因之間的加權相關性值構成的矩陣。

• TOM (Topological overlap matrix)：把鄰接矩陣轉換為拓撲重疊矩陣，以降低噪音和假相關，獲得的新距離矩陣，這個信息可拿來構建網絡或繪制TOM圖。

基本分析流程

1. 構建基因共表達網絡：使用加權的表達相關性。

2. 識別基因集：基于加權相關性，進行層級聚類分析，并根據設定標準切分聚類結果，獲得不同的基因模塊，用聚類樹的分枝和不同顏色表示。

3. 如果有表型信息，計算基因模塊與表型的相關性，鑒定性狀相關的模塊。

4. 研究模型之間的關系，從系統(tǒng)層面查看不同模型的互作網絡。

5. 從關鍵模型中選擇感興趣的驅動基因，或根據模型中已知基因的功能推測未知基因的功能。

6. 導出TOM矩陣，繪制相關性圖。

WGCNA包實戰(zhàn)

R包WGCNA是用于計算各種加權關聯(lián)分析的功能集合，可用于網絡構建，基因篩選，基因簇鑒定，拓撲特征計算，數據模擬和可視化等。

輸入數據和參數選擇

1. WGCNA本質是基于相關系數的網絡分析方法，適用于多樣品數據模式，一般要求樣本數多于15個。樣本數多于20時效果更好，樣本越多，結果越穩(wěn)定。

2. 基因表達矩陣: 常規(guī)表達矩陣即可，即基因在行，樣品在列，進入分析前做一個轉置。RPKM、FPKM或其它標準化方法影響不大，推薦使用Deseq2的varianceStabilizingTransformation或log2(x+1)對標準化后的數據做個轉換。如果數據來自不同的批次，需要先移除批次效應 (記得上次轉錄組培訓課講過如何操作)。如果數據存在系統(tǒng)偏移，需要做下quantile normalization。

3. 性狀矩陣：用于關聯(lián)分析的性狀必須是數值型特征 (如下面示例中的Height,?Weight,Diameter)。如果是區(qū)域或分類變量，需要轉換為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

4. 推薦使用Signed network和Robust correlation (bicor)。(這個根據自己的需要，看看上面寫的每個網絡怎么計算的，更知道怎么選擇)

5. 無向網絡在power小于15或有向網絡power小于30內，沒有一個power值可以使無標度網絡圖譜結構R^2達到0.8或平均連接度降到100以下，可能是由于部分樣品與其他樣品差別太大造成的。這可能由批次效應、樣品異質性或實驗條件對表達影響太大等造成, 可以通過繪制樣品聚類查看分組信息、關聯(lián)批次信息、處理信息和有無異常樣品 (可以使用之前講過的熱圖簡化，增加行或列屬性)。如果這確實是由有意義的生物變化引起的，也可以使用后面程序中的經驗power值。

安裝WGCNA

WGCNA依賴的包比較多，bioconductor上的包需要自己安裝，cran上依賴的包可以自動安裝。通常在R中運行下面4條語句就可以完成WGCNA的安裝。

建議在編譯安裝R時增加--with-blas --with-lapack提高矩陣運算的速度，具體見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實戰(zhàn)

實戰(zhàn)采用的是官方提供的清理后的矩陣，原矩陣信息太多，容易給人誤導，后臺回復WGCNA?獲取數據。

數據讀入

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ī)表達矩陣，log2轉換后或
# Deseq2的varianceStabilizingTransformation轉換的數據
# 如果有批次效應，需要事先移除，可使用removeBatchEffect
# 如果有系統(tǒng)偏移(可用boxplot查看基因表達分布是否一致)，
# 需要quantile normalization

exprMat <- "WGCNA/LiverFemaleClean.txt"

# 官方推薦 "signed" 或 "signed hybrid"
# 為與原文檔一致，故未修改
type = "unsigned"

# 相關性計算
# 官方推薦 biweight mid-correlation & bicor
# corType: pearson or bicor
# 為與原文檔一致，故未修改
corType = "pearson"

corFnc = ifelse(corType=="pearson", cor, bicor)
# 對二元變量，如樣本性狀信息計算相關性時，
# 或基因表達嚴重依賴于疾病狀態(tài)時，需設置下面參數
maxPOutliers = ifelse(corType=="pearson",1,0.05)

# 關聯(lián)樣品性狀的二元變量時，設置
robustY = ifelse(corType=="pearson",T,F)

##導入數據##
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

數據篩選

## 篩選中位絕對偏差前75%的基因，至少MAD大于0.01
## 篩選后會降低運算量，也會失去部分信息
## 也可不做篩選，使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)),]

## 轉換為樣品在行，基因在列的矩陣
dataExpr <- as.data.frame(t(dataExprVar))

## 檢測缺失值
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="")

軟閾值的篩選原則是使構建的網絡更符合無標度網絡特征。

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)，縱軸是無標度網絡的評估參數，數值越高，
# 網絡越符合無標度特征 (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")
# 篩選標準。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

經驗power?(無滿足條件的power時選用)

# 無向網絡在power小于15或有向網絡power小于30內，沒有一個power值可以使
# 無標度網絡圖譜結構R^2達到0.8，平均連接度較高如在100以上，可能是由于
# 部分樣品與其他樣品差別太大。這可能由批次效應、樣品異質性或實驗條件對
# 表達影響太大等造成?？梢酝ㄟ^繪制樣品聚類查看分組信息和有無異常樣品。
# 如果這確實是由有意義的生物變化引起的，也可以使用下面的經驗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)) ? ? ?
 ? ? ? ? ?)
 ? ? ? ? ?)
}

網絡構建

##一步法網絡構建：One-step network construction and module detection##
# power: 上一步計算的軟閾值
# maxBlockSize: 計算機能處理的最大模塊的基因數量 (默認5000)；
#  4G內存電腦可處理8000-10000個，16G內存電腦可以處理2萬個，32G內存電腦可
# ?以處理3萬個
# ?計算資源允許的情況下最好放在一個block里面。
# corType: pearson or bicor
# numericLabels: 返回數字而不是顏色作為模塊的名字，后面可以再轉換為顏色
# saveTOMs：最耗費時間的計算，存儲起來，供后續(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, loadTOMs=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...

# 根據模塊中基因數目的多少，降序排列，依次編號為 `1-最大模塊數`。
# **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

層級聚類樹展示各個模塊

## 灰色的為**未分類**到模塊的基因。
# Convert labels to colors for plotting
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
# Plot the dendrogram and the module colors underneath
# 如果對結果不滿意，還可以recutBlockwiseTrees，節(jié)省計算時間
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
 ? ? ? ? ? ? ? ? ? ?"Module colors",
 ? ? ? ? ? ? ? ? ?  dendroLabels = FALSE, hang = 0.03,
 ? ? ? ? ? ? ? ? ?  addGuide = TRUE, guideHang = 0.05)

繪制模塊之間相關性圖

# module eigengene, 可以繪制線圖，作為每個模塊的基因表達趨勢的展示
MEs = net$MEs

### 不需要重新計算，改下列名字就好
### 官方教程是重新計算的，起始可以不用這么麻煩
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(
 ?as.numeric(str_replace_all(colnames(MEs),"ME",""))))
MEs_col = orderMEs(MEs_col)

# 根據基因間表達量進行聚類所得到的各模塊間的相關性圖
# marDendro/marHeatmap 設置下、左、上、右的邊距
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap",
 ? ? ? ? ? ? ? ? ? ? ?marDendro = c(3,3,2,4),
 ? ? ? ? ? ? ? ? ? ? ?marHeatmap = c(3,4,2,2), plotDendrograms = T,
 ? ? ? ? ? ? ? ? ? ? ?xLabelsAngle = 90)

## 如果有表型數據，也可以跟ME數據放一起，一起出圖
#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)

可視化基因網絡 (TOM plot)

# 如果采用分步計算，或設置的blocksize>=總基因數，直接load計算好的TOM結果
# 否則需要再計算一遍，比較耗費時間
# 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

# 這一部分特別耗時，行列同時做層級聚類
TOMplot(plotTOM, net$dendrograms, moduleColors,
 ? ? ? ?main = "Network heatmap plot, all genes")

導出網絡用于Cytoscape

Cytoscape繪制網絡圖見我們更新版的視頻教程或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 默認為0.5, 可以根據自己的需要調整，也可以都導出后在
# cytoscape中再調整
cyt = exportNetworkToCytoscape(TOM,
 ? ? ? ? ? ? edgeFile = paste(exprMat, ".edges.txt", sep=""),
 ? ? ? ? ? ? nodeFile = paste(exprMat, ".nodes.txt", sep=""),
 ? ? ? ? ? ? weighted = TRUE, threshold = 0,
 ? ? ? ? ? ? nodeNames = probes, nodeAttr = moduleColors)

關聯(lián)表型數據

trait <- "WGCNA/TraitsClean.txt"
# 讀入表型數據，不是必須的
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)), ]
}

### 模塊與表型數據關聯(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表示保留幾位小數
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"))

模塊內基因與表型數據關聯(lián), 從上圖可以看到MEmagenta與Insulin_ug_l相關，選取這兩部分進行分析。

## 從上圖可以看到MEmagenta與Insulin_ug_l相關

## 模塊內基因與表型數據關聯(lián)

# 性狀跟模塊雖然求出了相關性，可以挑選最相關的那些模塊來分析，
# 但是模塊本身仍然包含非常多的基因，還需進一步的尋找最重要的基因。
# 所有的模塊都可以跟基因算出相關系數，所有的連續(xù)型性狀也可以跟基因的表達
# 值算出相關系數。
# 如果跟性狀顯著相關基因也跟某個模塊顯著相關，那么這些基因可能就非常重要
# 。

### 計算模塊與基因的相關性矩陣

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
}

# 計算性狀與基因的相關性矩陣

## 只有連續(xù)型性狀才能進行計算，如果是離散變量，在構建樣品表時就轉為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.

# 最后把兩個相關性矩陣聯(lián)合起來,指定感興趣模塊進行分析
module = "magenta"
pheno = "Insulin_ug_l"
modNames = substring(colnames(MEs_col), 3)
# 獲取關注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 獲取模塊內的基因
moduleGenes = moduleColors == module

sizeGrWindow(7, 7)
par(mfrow = c(1,1))
# 與性狀高度相關的基因，也是與性狀相關的模型的關鍵基因
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)

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

### 計算鄰接矩陣
adjacency = adjacency(dataExpr, power = power)

### 把鄰接矩陣轉換為拓撲重疊矩陣，以降低噪音和假相關，獲得距離矩陣。
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM

### 層級聚類計算基因之間的距離樹
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)

### 通過計算模塊的代表性模式和模塊之間的定量相似性評估，合并表達圖譜相似
的模塊
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

## 分步法完結

Reference:

1. 官網：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/

2. 術語解釋：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/Simulated-00-Background.pdf

3. FAQ：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html

4. 生信博客：http://blog.genesino.com

5. 高顏值免費在線繪圖工具升級版來了~~~ （可以在線做 WGCNA 分析）
   

你可能還想看

• 自從用了這個神器，大規(guī)模RNA-seq數據挖掘我也可以

• NGS基礎 - FASTQ格式解釋和質量評估

• NGS基礎 - 高通量測序原理

• NGS基礎 - 參考基因組和基因注釋文件

• NGS基礎 - GTF/GFF文件格式解讀和轉換

• NGS基礎 - 測序原始數據下載

• WGCNA教程流程圖：
  

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集