WGCNA（2b）：分步法完成网络构建和模块检测

在前两篇帖子中介绍了数据的导入清洗，以及一步法构建网络，这篇文章将介绍网络构建的第二种方法，分步法。
WGCNA（1）：R包安装及数据导入清洗 - 简书 (jianshu.com)
WGCNA（2a）：一步法完成网络构建和模块检测 - 简书 (jianshu.com)

参考材料还是官方说明书：
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-man.pdf

1. 准备工作

$ cd ~/WGCNA

$ R
> library(WGCNA)
> enableWGCNAThreads(15)

#载入第一步中的表达量和表型值
> lnames = load(file = "WGCNA0.3-dataInput.RData")
> lnames
[1] "datExpr"  "datTraits"

2. 分步法构建基因网络和模块识别

2.1 确定合适的软阈值：网络拓扑分析

# 设置网络构建参数选择范围
> powers = c(c(1:10), seq(from = 12, to=20, by=2)) 
> powers
 [1]  1  2  3  4  5  6  7  8  9 10 12 14 16 18 20
# 计算无尺度分布拓扑矩阵
> sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

软阈值计算结果可视化

> pdf(" Analysis of network topology for various soft-thresholding powers.pdf",width = 9, height=5)
#一页多图，一行2列
> par(mfrow = c(1,2))
#字号
> cex1 = 0.9
# Scale-free topology fit index as a function of the soft-thresholding power（左图）
> 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") 
#查看最佳软阈值
> sft$powerEstimate
[1] 6 
# this line corresponds to using an R^2 cut-off of h
> abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power（右图）
> 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")
> dev.off()

Figure 1: Analysis of network topology for various soft-thresholding powers. The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis).

2.2 共表达相似性和邻接性

通过上述计算得知，最佳软阈值为6。

> softPower = 6
> adjacency = adjacency(datExpr, power = softPower)

2.3 拓扑覆盖矩阵TOM

为了减小噪声和虚假关联的影响，我们将邻接变换为拓扑重叠矩阵，并计算相应的不相似度。

> TOM = TOMsimilarity(adjacency)
> dissTOM = 1-TOM

2.4 通过TOM进行聚类

现在使用层次聚类来产生基因的层次聚类树(树状图)。注意，我们使用了hclust函数，它提供了比标准hclust函数更快的分层集群例程。

# 调用层次聚类函数
> geneTree = hclust(as.dist(dissTOM), method = "average")
# 聚类结果绘图
> pdf(" Gene clustering on TOM-based dissimilarity.pdf",width = 12, height=9)
> plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)
> dev.off()

在下图中，每个短分枝代表一个基因，相邻的基因表示他们高度共表达。聚类方法是dynamicTreeCut包中的Dynamic TreeCut。

Figure 2: Clustering dendrogram of genes, with dissimilarity based on topological overlap

# 由于基因数目比较多，我们喜欢大的模块，所以我们把最小模块的大小设置得比较高，这个数值要根据自己的数据进行修改
> minModuleSize = 200
# 模块识别使用动态树切割，这里敏感度deepSplit我设置了1，需要根据自己数据进行调试，默认是2
> dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 1, pamRespectsDendro = FALSE, minClusterSize = minModuleSize) 
# 查看模块数
> table(dynamicMods)
# 转换数字标签为颜色
> dynamicColors = labels2colors(dynamicMods)
# 查看颜色
> table(dynamicColors)
# 绘图
> pdf("Gene dendrogram and module colors.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()

Figure 3: Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned module colors.

2.5 合并表达模式相似的模块

量化各个模块的共表达相似性，根据模块间的相关性合并表达模式相似的模块。

# Calculate eigengenes
> 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")
# 画图
> pdf("Clustering of module eigengenes.pdf",width = 7, height=6)
> plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

这里选择切割高度0.3，即模块之间相关性达到0.7进行合并，这个系数也要根据自己的数据进行更改，常见的是0.25（Fig.4）。

> MEDissThres = 0.3
# 在树状图中加入切割线
> abline(h=MEDissThres, col = "red")
# 调用自动归并函数
> merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# 合并模块的颜色
> mergedColors = merge$colors
# 合并后新模块
> mergedMEs = merge$newMEs

Figure 4: Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned merged module colors and the original module colors.

为了查看合并对模块颜色的影响，我们再次绘制了基因树状图，下面依次是原始和合并后的模块颜色（Fig.5）。

> pdf("PlotsgeneDendro3.pdf",width = 12, height=9)
> plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)
> dev.off()

在随后的分析中，我们将在mergedColors中使用合并的模块颜色，保存相关的变量，以便在后续分析中使用。

# 对模块颜色重命名
> moduleColors = mergedColors
# 将颜色转换为数据标签
> colorOrder = c("grey", standardColors(50))
> moduleLabels = match(moduleColors, colorOrder)-1
> MEs = mergedMEs
# 保存模块颜色和标签，以供后续部分使用
> save(MEs, moduleLabels, moduleColors, geneTree, file = "networkConstruction-stepByStep.RData")

Figure 5: Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned merged module colors and the original module colors.