1准备阶段：加载包、数据集
library(WGCNA)
# Load the data saved in the first part
lnames = load(file = "FemaleLiver-01-dataInput.RData")
2阈值选取
based on the criterion of approximate scale-free topology。使用pickSoftThreshold()函数进行网络拓扑的分析，得到备选软阈值对应的相关数值，如signed R^2
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
sizeGrWindow(9, 5)
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");
# 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")
得到下图的结果，此处设置的高度为0.9，达到这个高度的最小候选阈值为6，因此，我们选择软阈值为6.
3网络构建+模块识别
本文介绍两种方法构建一个加权共表达网络，第一种方法自动构建网络及模块识别（也就是说，这个方法使用一个函数就能实现），第二种方法通过分步骤的方式来构建网络及模块识别。
（1）自动的方式构建网络+模块识别
网络构建及模块识别：使用blockwiseModules()函数，将网络构建和模块识别的多步整合为一步。即自动构建一个加权网络，创建一个聚类树，将树形图的分枝定义为模块，将相近的模块合并成一个模块。这个函数输出：模块颜色、模块特征基因（module eigengenes）。这些量可以用于后续的分析。同时，我们可以将模块识别的结果进行可视化，即在聚类的树形图下绘制每个基因所属的分组。
net = blockwiseModules(datExpr, power = 6,
TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = TRUE,
saveTOMFileBase = "femaleMouseTOM",
verbose = 3)
此处选择的软阈值为6，设置模块中包含的基因个数最小为30（因为我们喜欢大的模块，因此这个值应该设置的尽可能大一些），一个中等的敏感度（deepSplit=2）。参数mergeCutHeight是模块融合的阈值。这个函数返回的是数值，不是模块的颜色标签。
net＄colors：显示了数据集中的每个基因被分配到什么组中，不同的组用不同的数字表示，比如1表示组1，2表示组2。。。
net＄MEs：网络中每个模块的特征值。
使用table(net＄colors)可以显示网络总共有多少个模块（分组），每个模块（组）中包含有多少个基因。
table(net$colors)
#
根据模块的大小（模块中包含的变量个数），对模块进行排序，分别对模块使用数字1-18进行标示。将不属于1-18模块的基因放到标签为0的模块中，也就是说，标签为0的模块中的基因，是没有被分配到任何模块的基因。
绘制基因聚类的树形图，在树形图下显示基因所属的分组。
# open a graphics window
sizeGrWindow(12, 9)
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
结果如下：
将模块相关信息进行保存，用于后续的分析。
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
file = "FemaleLiver-02-networkConstruction-auto.RData")
保存了每个基因属于哪个组（前两行，只不过第一行代码中组用数字表示，第二行代码中组用颜色表示，labels2colors()函数含义很明显，label_to_colors。）；第三行代码是保存每个组的特征向量，第四行是保存聚类的树形图。需要保存的变量有①模块的特征基因②模块的数字标签③模块的颜色标签④基因的树形图。
（2）逐步构建网络+模块识别
Step_1:Co-expression similarity and adjacency
softPower = 6
adjacency = adjacency(datExpr, power = softPower)
Step_2:计算拓扑重叠矩阵（TOM）
# Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
Step_3:使用TOM(拓扑重叠矩阵)进行聚类，绘制聚类得到的树形图。
# Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average")
# Plot the resulting clustering tree (dendrogram)
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity", labels = FALSE, hang = 0.04)
Step_4:使用dynamic tree cut来识别模块。
需要事前设置好最小的模块允许包含的变量个数。因为，我们喜欢大的模块，此处设置模块的最小值为30.使用cutreeDynamic()函数对树形图进行切割。
# 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)
table(dynamicMods)
在树形图下方绘制模块，使用plotDendroAndColors()函数
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
Step_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")
# Plot the result
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
We choose a height cut of 0.25, corresponding to correlation of 0.75, to merge
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# 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 modules:
mergedMEs = merge$newMEs
为了查看合并模块后，网络中模块有什么变化。我们再次绘制基因树形图，并在树形图下面绘制原始模块颜色和合并后的模块颜色.
sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)
#dev.off()
保存模块相关变量，用于后续的分析.需要保存的变量有①模块的特征基因②模块的数字标签③模块的颜色标签④基因的树形图。
# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50))
moduleLabels = match(moduleColors, colorOrder)-1
MEs = mergedMEs
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree, file = "FemaleLiver-02-networkConstruction-stepByStep.RData")
*注：本文中谈到的组或者是模块是一个概念。
