一、前言

1.1 概念

加权基因共表达网络分析 (WGCNA, Weighted correlation network analysis)是用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集, 并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。

该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。

2.2 原理

从方法上来讲，WGCNA分为表达量聚类分析和表型关联两部分，主要包括基因之间相关系数计算、基因模块的确定、共表达网络、模块与性状关联四个步骤。

第一步计算任意两个基因之间的相关系数（Person Coefficient）。为了衡量两个基因是否具有相似表达模式，一般需要设置阈值来筛选，高于阈值的则认为是相似的。但是这样如果将阈值设为0.8，那么很难说明0.8和0.79两个是有显著差别的。因此，WGCNA分析时采用相关系数加权值，即对基因相关系数取N次幂，使得网络中的基因之间的连接服从无尺度网络分布(scale-freenetworks)，这种算法更具生物学意义。

第二步通过基因之间的相关系数构建分层聚类树，聚类树的不同分支代表不同的基因模块，不同颜色代表不同的模块。基于基因的加权相关系数，将基因按照表达模式进行分类，将模式相似的基因归为一个模块。这样就可以将几万个基因通过基因表达模式被分成了几十个模块，是一个提取归纳信息的过程。

•共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算 (冥次的值也就是软阈值 (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

•Module(模块)：高度內连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块。

•Eigengene （eigen + gene）：基因和样本构成的矩阵，https://en.wiktionary.org/wiki/eigengene

•Connectivity (连接度)：类似于网络中 "度" (degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。

•Module eigengene E: 给定模型的第一主成分，代表整个模型的基因表达谱。这个是个很巧妙的梳理，我们之前讲过PCA分析的降维作用，之前主要是拿来做可视化，现在用到这个地方，很好的用一个向量代替了一个矩阵，方便后期计算。(降维除了PCA，还可以看看tSNE)

•Intramodular connectivity: 给定基因与给定模型内其他基因的关联度，判断基因所属关系。

•Module membership: 给定基因表达谱与给定模型的eigengene的相关性。

•Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。

•Adjacency matrix (邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。

•TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

以上解释来自：WGCNA分析，简单全面的最新教程；https://www.jianshu.com/p/e9cc3f43441d

分析案例

WGCNA分析，如果你的数据样本较大，建议在服务器中进行运行！！

二、WGCNA 代码教程

2.1 加载所需的R包

加载WGCNAR包

#install.packages("WGCNA") 
#BiocManager::install('WGCNA') 
library(WGCNA)

进行基础设置，（此步骤可以省略）

options(stringsAsFactors = FALSE) 
enableWGCNAThreads() ## 打开多线程

2.2 加载数据矩阵

WGCNA.fpkm = read.table("WGCNA.inputdata.txt",header=T, 
                        comment.char = "", 
                        check.names=F)

数据矩阵格式

> WGCNA.fpkm[1:10,1:10] 
sample MT_0_1   MT_0_2   MT_0_3   MT_3_1   MT_3_2   MT_3_3   MT_6_1   MT_6_2   MT_6_3 
1  gene10      0 0.000000 0.000000 0.000000 0.000000 0.000000 0.225467 0.000000 0.116728 
2  gene12      0 0.000000 0.000000 0.000000 0.140796 0.000000 0.051667 0.076227 0.109181 
3  gene13      0 0.000000 0.000000 0.000000 0.000000 0.000000 0.090672 0.000000 0.000000 
4  gene14      0 0.000000 0.247270 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
5  gene16      0 0.000000 0.213587 0.051469 0.433966 0.786561 0.000000 0.246451 0.972179 
6  gene19      0 0.113876 0.060553 0.000000 0.331847 0.000000 0.099846 0.224528 0.191205 
7  gene20      0 0.000000 0.000000 0.183355 0.531899 0.383033 0.000000 2.460091 2.441143 
8  gene21      0 0.470667 0.000000 1.016470 0.583894 0.000000 1.074058 2.020643 1.007870 
9  gene24      0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
10 gene25      0 0.054502 0.026962 0.026076 0.025938 0.000000 0.028589 0.023683 0.000000

查看数据名称

names(WGCNA.fpkm) 
 [1] "sample"  "MT_0_1"  "MT_0_2"  "MT_0_3"  "MT_3_1"  "MT_3_2"  "MT_3_3"  "MT_6_1"  "MT_6_2"  "MT_6_3"  "MT_12_1" "MT_12_2" 
[13] "MT_12_3" "HG_0_1"  "HG_0_2"  "HG_0_3"  "HG_3_1"  "HG_3_2"  "HG_3_3"  "HG_6_1"  "HG_6_2"  "HG_6_3"  "HG_12_1" "HG_12_2" 
[25] "HG_12_3"

对数据进行提取

names(WGCNA.fpkm) 
datExpr0 = as.data.frame(t(WGCNA.fpkm[,-1])) 
names(datExpr0) = WGCNA.fpkm$sample; 
rownames(datExpr0) = names(WGCNA.fpkm[,-1])
check missing value
gsg = goodSamplesGenes(datExpr0, verbose = 3) 
gsg$allOK 
if (!gsg$allOK) 
{ 
  # Optionally, print the gene and sample names that were removed: 
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", "))) 
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", "))) 
  # Remove the offending genes and samples from the data: 
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes] 
}

filter

meanFPKM= 0.5  ####--过滤标准，可以修改 
n=nrow(datExpr0) 
datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean) 
datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM] 
# for meanFpkm in row n+1 and it must be above what you set--select meanFpkm>opt$meanFpkm(by rp) 
filtered_fpkm=t(datExpr0) 
filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm) 
names(filtered_fpkm)[1]="sample" 
head(filtered_fpkm) 
write.table(filtered_fpkm, file="mRNA.symbol.uniq.filter.txt", 
            row.names=F, col.names=T,quote=FALSE,sep="\t")

2.3 对input data进行聚类

#############Sample cluster########### 
sampleTree = hclust(dist(datExpr0), method = "average") 
pdf(file = "1.sampleClustering.pdf", width = 15, height = 8) 
par(cex = 0.6) 
par(mar = c(0,6,6,0)) 
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1.5, cex.main = 2) 
### Plot a line to show the cut 
#abline(h = 180, col = "red")##剪切高度不确定，故无红线 
dev.off()

如果数据中有个别数据与其他数据有较大的差异性，为了让后续分析更精准，可以将其进行适当的删除即可。如上图所示.

2.4 加载性状数据

allTraits = read.table("type.txt",row.names=1,header=T,comment.char = "",check.names=F) 
dim(allTraits) 
names(allTraits) 
head(allTraits)

> head(allTraits) 
       Cold_3h Cold_6h Cold_12 
MT_0_1       0       0       0 
MT_0_2       0       0       0 
MT_0_3       0       0       0 
MT_3_1       3       0       0 
MT_3_2       3       0       0 
MT_3_3       3       0       0

表达量与性状数据进行匹配

fpkmSamples = rownames(datExpr0) 
traitSamples =rownames(allTraits) 
traitRows = match(fpkmSamples, traitSamples) 
datTraits = allTraits[traitRows,] 
rownames(datTraits) 
collectGarbage()

pdf(file="2.Sample_dendrogram_and_trait_heatmap.pdf",width=20,height=12) 
plotDendroAndColors(sampleTree2, traitColors, 
                    groupLabels = names(datTraits), 
                    main = "Sample dendrogram and trait heatmap",cex.colorLabels = 1.5, cex.dendroLabels = 1, cex.rowText = 2) 
dev.off()

2.5 选择最优的sftpower值

# 选择最优的sftpower值 
# softPower = sft$powerEstimate 
# NA 
softPower = 9

经验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))        
          ) 
          ) 
}

2.6 网络构建

adjacency = adjacency(datExpr0, power = softPower) 
TOM = TOMsimilarity(adjacency); 
dissTOM = 1-TOM 

# Call the hierarchical clustering function 
geneTree = hclust(as.dist(dissTOM), method = "average"); 
# Plot the resulting clustering tree (dendrogram) 
 
#sizeGrWindow(12,9) 
pdf(file="4_Gene clustering on TOM-based dissimilarity.pdf",width=24,height=18) 
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", 
     labels = FALSE, hang = 0.04) 
dev.off()

层级聚类树展示各个模块

# 选择模块的数量，可选大小 
minModuleSize = 30 
# Module identification using dynamic tree cut: 
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, 
                            deepSplit = 2, pamRespectsDendro = FALSE, 
                            minClusterSize = minModuleSize); 
table(dynamicMods) 
 
# Convert numeric lables into colors 
dynamicColors = labels2colors(dynamicMods) 
table(dynamicColors) 
# Plot the dendrogram and colors underneath 
#sizeGrWindow(8,6) 
pdf(file="5_Dynamic Tree Cut.pdf",width=8,height=6) 
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", 
                    dendroLabels = FALSE, hang = 0.03, 
                    addGuide = TRUE, guideHang = 0.05, 
                    main = "Gene dendrogram and module colors") 
dev.off()

模块间进行矫正合并
计算模块合并的高度

# Calculate eigengenes 
MEList = moduleEigengenes(datExpr0, colors = dynamicColors) 
MEs = MEList$eigengenes 
# Calculate dissimilarity of module eigengenes 
MEDiss = 1-cor(MEs); 
# Cluster module eigengenes 
METree = hclust(as.dist(MEDiss), method = "average") 
# Plot the result 
#sizeGrWindow(7, 6) 
pdf(file="6_Clustering of module eigengenes.pdf",width=7,height=6) 
plot(METree, main = "Clustering of module eigengenes", 
     xlab = "", sub = "") 
MEDissThres = 0.4######剪切高度可修改 
# Plot the cut line into the dendrogram 
abline(h=MEDissThres, col = "red") 
dev.off()

重新绘制合并后的模块层次聚类图

# Call an automatic merging function 
merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3) 
# The merged module colors 
mergedColors = merge$colors 
# Eigengenes of the new merged modules: 
mergedMEs = merge$newMEs 
table(mergedColors) 
 
#sizeGrWindow(12, 9) 
pdf(file="7_merged dynamic.pdf", width = 9, height = 6) 
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), 
                    c("Dynamic Tree Cut", "Merged dynamic"), 
                    dendroLabels = FALSE, hang = 0.03, 
                    addGuide = TRUE, guideHang = 0.05) 
dev.off()

2.7 各个模块与数量性状间的相关性计算

## 模块与数量性状间的相关性 
# 
moduleColors = mergedColors 
# Construct numerical labels corresponding to the colors 
colorOrder = c("grey", standardColors(50)) 
moduleLabels = match(moduleColors, colorOrder)-1 
MEs = mergedMEs 
 
nGenes = ncol(datExpr0) 
nSamples = nrow(datExpr0) 
moduleTraitCor = cor(MEs, datTraits, use = "p") 
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples) 
 
 
pdf(file="8_Module-trait relationships.pdf",width=10,height=10) 
 
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", 
                   signif(moduleTraitPvalue, 1), ")", sep = "") 
 
dim(textMatrix) = dim(moduleTraitCor) 
par(mar = c(6, 8.5, 3, 3)) 
 
labeledHeatmap(Matrix = moduleTraitCor, 
               xLabels = names(datTraits), 
               yLabels = names(MEs), 
               ySymbols = names(MEs), 
               colorLabels = FALSE, 
               colors = greenWhiteRed(50), 
               textMatrix = textMatrix, 
               setStdMargins = FALSE, 
               cex.text = 0.5, 
               zlim = c(-1,1), 
               main = paste("Module-trait relationships")) 
dev.off()

2.8 可视化基因网络 (TOM plot)

数据处理

nGenes = ncol(datExpr0) 
nSamples = nrow(datExpr0) 
nSelect = 400 
# For reproducibility, we set the random seed  
set.seed(10) 
select = sample(nGenes, size = nSelect) 
selectTOM = dissTOM[select, select] 
# 
selectTree = hclust(as.dist(selectTOM), method = "average") 
selectColors = moduleColors[select] 
 
 
plotDiss = selectTOM^7 
diag(plotDiss) = NA
library("gplots") 
pdf(file="Network heatmap plot_selected genes.pdf",width=9, height=9) 
mycol = colorpanel(250,'red','orange','lemonchiffon') 
TOMplot(plotDiss, selectTree, selectColors, col=mycol ,main = "Network heatmap plot, selected genes") 
dev.off()

2.9 绘制模块之间相关性图

pdf(file="Eigengene dendrogram and Eigengene adjacency heatmap.pdf", width=5, height=7.5) 
par(cex = 0.9) 
plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90) 
dev.off()

pdf(file="Eigengene dendrogram_2.pdf",width=6, height=6) 
par(cex = 1.0) 
plotEigengeneNetworks(MEs, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE) 
dev.off()

pdf(file="Eigengene adjacency heatmap_2.pdf",width=6, height=6) 
par(cex = 1.0) 
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90) 
dev.off()

参考：

1.一文看懂WGCNA 分析(2019更新版)
2.WGCNA分析，简单全面的最新教程

“小杜的生信筆記”公众号、知乎、简书平台，主要发表或收录生物信息学的教程，以及基于R的分析和可视化（包括数据分析，图形绘制等）；分享感兴趣的文献和学习资料！!