WGCNA分析，簡單全面的最新教程

05-20

WGCNA分析，簡單全面的最新教程

來自專欄生信寶典

WGCNA分析，簡單全面的最新教程

Jump to...

WGCNA基本概念
基本分析流程
WGCNA包實戰
輸入數據和參數選擇

安裝WGCNA
WGCNA實戰

數據讀入
軟閾值篩選
經驗power (無滿足條件的power時選用)
網路構建
層級聚類樹展示各個模塊
繪製模塊之間相關性圖
可視化基因網路 (TOM plot)
導出網路用於Cytoscape
關聯表型數據
分步法展示每一步都做了什麼

Reference:

WGCNA基本概念

加權基因共表達網路分析 (WGCNA, Weighted correlation network

analysis)是用來描述不同樣品之間基因關聯模式的系統生物學方法，可以用來鑒定高度協同變化的基因集,

並根據基因集的內連性和基因集與表型之間的關聯鑒定候補生物標記基因或治療靶點。

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

理解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. 模塊與性狀進行關聯分析，找出與關注性狀相關度最高的模塊；3. 模塊與樣本進行關聯分析，找到樣品特異高表達的模塊。

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

Connectivity (連接度)：類似於網路中「度」
(degree)的概念。每個基因的連接度是與其相連的基因的邊屬性之和。

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

Intramodular connectivity:

給定基因與給定模型內其他基因的關聯度，判斷基因所屬關係。

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

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

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

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

基本分析流程

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

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

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

研究模型之間的關係，從系統層面查看不同模型的互作網路。

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

導出TOM矩陣，繪製相關性圖。

WGCNA包實戰

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

輸入數據和參數選擇

WGCNA本質是基於相關係數的網路分析方法，適用於多樣品數據模式，一般要求樣本數多於15個。樣本數多於20時效果更好，樣本越多，結果越穩定。
基因表達矩陣:
常規表達矩陣即可，即基因在行，樣品在列，進入分析前做一個轉置。RPKM、FPKM或其它標準化方法影響不大，推薦使用Deseq2的varianceStabilizingTransformation或log2(x+1)對標準化後的數據做個轉換。如果數據來自不同的批次，需要先移除批次效應
(記得上次轉錄組培訓課講過如何操作)。如果數據存在系統偏移，需要做下quantile normalization。
性狀矩陣：用於關聯分析的性狀必須是數值型特徵
(如下面示例中的Height, Weight,
Diameter)。如果是區域或分類變數，需要轉換為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)。(這個根據自己的需要，看看上面寫的每個網路怎麼計算的，更知道怎麼選擇)

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

實戰採用的是官方提供的清理後的矩陣，原矩陣信息太多，容易給人誤導，後台回復

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=<number_of_processors>

## *

## * 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.

# 常規表達矩陣，log2轉換後或

# Deseq2的varianceStabilizingTransformation轉換的數據

# 如果有批次效應，需要事先移除，可使用removeBatchEffect

# 如果有系統偏移(可用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)

# 對二元變數，如樣本性狀信息計算相關性時，

# 或基因表達嚴重依賴於疾病狀態時，需設置下面參數

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

# 關聯樣品性狀的二元變數時，設置

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

##導入數據##

dataExpr <- read.table(exprMat, sep= , 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：最耗費時間的計算，存儲起來，供後續使用

# 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，節省計算時間

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)

關聯表型數據

trait <- "WGCNA/TraitsClean.txt"

# 讀入表型數據，不是必須的

if(trait != "") {

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

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

sampleName = rownames(dataExpr)

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

}

### 模塊與表型數據關聯

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), "
(", 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"))

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

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

## 模塊內基因與表型數據關聯

# 性狀跟模塊雖然求出了相關性，可以挑選最相關的那些模塊來分析，

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

# 所有的模塊都可以跟基因算出相關係數，所有的連續型性狀也可以跟基因的表達

# 值算出相關係數。

# 如果跟性狀顯著相關基因也跟某個模塊顯著相關，那麼這些基因可能就非常重要

# 。

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

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

}

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

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

# 最後把兩個相關性矩陣聯合起來,指定感興趣模塊進行分析

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
"),

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:

官網
https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
術語解釋
https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/Simulated-00-Background.pdf
FAQ
https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/faq.html
生信博客 http://blog.genesino.com

WGCNA分析，簡單全面的最新教程

WGCNA分析，簡單全面的最新教程

WGCNA分析，簡單全面的最新教程

WGCNA基本概念

基本分析流程

WGCNA包實戰

輸入數據和參數選擇

安裝WGCNA

WGCNA實戰

數據讀入

軟閾值篩選

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

網路構建

層級聚類樹展示各個模塊

繪製模塊之間相關性圖

可視化基因網路 (TOM plot)

導出網路用於Cytoscape

關聯表型數據

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

Reference:

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