2. 共表达网络构建及模块获取

基于前面处理的数据，获得共表达网络及模块主要有三种方式：

• 自动法：采用自动的网络构建和模块检测功能，快速获得结果
• 逐步法：采用分布的网络构建和模块检测功能，以便使用者在每一步中可以选择自定义的方法或者可替代的方法
• 分块法：一个自动的分块的网络建设和模块检测方法，解决分析数据集太大，不能在一次运行的情况

关于这三个方法的关系，我们可以这样理解：

1. 如果不像太折腾，就使用自动法，快速构建共表达网络并识别模块
2. 如果需要根据自己数据集的特征调整参数，要么调整自动法的函数的参数，或者采用逐步法，在每一步的处理中设置合适的参数
3. 在遇到要处理的基因太多的时候，就采用分块法，这是在自动法的基础上，先将数据集分成几块，然后分别计算，并在最后合并。

在实践操作中，大家选择三个方法中的一个进行计算。我们一般常用的方法是自动法。

2.1. 导入处理的数据：

首先，导入准备好的样本表达数据datExpr和临床信息datTraits：

# 设定工作目录
setwd("/home/cxt/data/R/")
# 导入包
library(WGCNA)
# 设置将strings不要转成factors
options(stringsAsFactors=F)

# 允许使用最大线程
allowWGCNAThreads()
# 或者直接指定线程数
# enableWGCNAThreads(nThreads=16)

# 导入数据
lnames = load(file = "dataInput.RData");
lnames

2.2. 自动法：网络构建和模块检测

2.2.1 确定软阈值β

当前的数据是一个基因表达矩阵,然后接着需要构建相关性矩阵。相关性矩阵有两种，分别是Signed相关性矩阵（区分正相关和负相关）和Unsigned相关性矩阵（不区分正负）。接着需要将相关性矩阵构建成邻接矩阵。由于设置阈值时具有人为干扰，WGCNA推荐使用软阈值β（soft-thresholding power）。原则上，需要尽可能接近无尺度网络（Scale independence第一次达到0.8以上），尽可能保留连通性信息。 这里我们使用pickSoftThreshold函数，来获取合适的软阈值β：

# 软阈值的预设范围
powers <- c(c(1:10), seq(from=12, to=20, by=2))
# 自动计算推荐的软阈值
sft <- pickSoftThreshold(datExpr, powerVector=powers, verbose=5, networkType="unsigned")
# 推荐值。如果是NA，就需要画图来自己挑选
sft$powerEstimate

# 作图
sizeGrWindow(9, 5)
par(mfrow=c(1, 2))
cex1 <- 0.9

# Scale-free topology fit index 作为评价软阈值指标
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")
# h值作为筛选指标
abline(h=0.80, col="red")

# Mean connectivity 作为评价软阈值的指标
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")

# 如果是NA，此时就需要根据图来自己指定数值了
sft$powerEstimate <- 6

2.2.2 构建网络并识别模块

直接调用blockwiseModules函数构建网络并确定模块

# 构建网络
# deepSplit是生成模块的梯度，从[0:4]中选择，越大模块越多
# minModuleSize最小模块的基因数目
# mergeCutHeight是合并相似的模块的合并系数，是通过主成分分析出来的
# mumericLabels 已数字命名模块
# nThreads 线程数，当设置为0时使用最大线程
net <- blockwiseModules(datExpr, power = sft$powerEstimate,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "femaleMouseTOM", 
                       verbose = 3)

# 查看每个模块的基因数，其中0模块下为没有计算进入模块的基因数
table(net$colors)

2.2.3. 结果可视化

# 打开新的绘图窗口
sizeGrWindow(12, 9)
# 把模块编号转成颜色
mergedColors <- labels2colors(net$colors)
# 绘制树状图和下面的模块颜色
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

2.2.4. 保存结果

# 计算模块特征向量MEs
moduleLabels <- net$colors
moduleColors <- labels2colors(net$colors)
MEs <- net$MEs
geneTree <- net$dendrograms[[1]]
# 保存数据
save(MEs, moduleLabels, moduleColors, geneTree,
     file="networkConstruction-auto.RData")

2.3. 逐步法：网络构建和模块检测

2.3.1. 确定软阈值β

与2.2.1 确定软阈值β相同

最终，我们确定了软阈值β为sft$powerEstimate中的值

2.3.2. 计算共表达相似性和邻接矩阵

在这步，根据基因-样本表达数据 ，计算基因之间的共表达相似性，并构建相应的基因-基因邻接矩阵

# 软阈值β为sft$powerEstimate
adjacency = adjacency(datExpr, power = sft$powerEstimate);

2.3.3. 构建拓扑重叠矩阵

为了使共表达关系中噪声（noise）和伪关联（spurious associations）的影响最小化，我们将邻接矩阵变换为拓扑重叠矩阵（Topological Overlap Matrix, TOM），计算相应的差异：

# Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM

2.3.4. 基因聚类分析

根据TOM矩阵，我们使用层次聚类来生成基因的层次聚类树(树状图)。

注意，我们使用的函数hclust提供了比标准hclust函数更快的层次集群例程。

# 调用层次聚类函数
geneTree = hclust(as.dist(dissTOM), method = "average");
# 画出层次聚类树
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04);

2.3.5. 模块识别

在得到的聚类树（树状图）中，每个叶子（即一条短垂直线）对应于一个基因。 树状图的分支一起密集互连，表示高度共表达的基因。 模块识别相当于识别各个分支（“从树状图中切割分支”）。 分枝切割有几种方法，我们的标准方法是使用dynamicTreeCut包中动态树剪切（Dynamic Tree Cut）方法：

# 如果想要大的模块，就把最小模块的的大小设置的高一点
minModuleSize = 30;
# 使用动态树剪切（Dynamic Tree Cut）确定模块
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                deepSplit = 2, pamRespectsDendro = FALSE,
                minClusterSize = minModuleSize);
table(dynamicMods)

上述函数返回了22个模块，并每个模块的标签是1-22。标签0表示未分配的基因。现在，我们将确认的模块绘制在距离树的下面：

# 将模块的数字标签转换为颜色标签
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# 绘制聚类树及其模块
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")

3.4.6. 模块合并

动态树剪切（Dynamic Tree Cut）方法虽然能确定表达数据高度相似的模块，但是模块之间的基因是高度共表达的，所以在合并模块的过程中是谨慎的。所以，为例确定模块间的共表达相似度，我们根据模块之间的连接性，确定他们之间的特征基因（Eigengenes）和聚类：

# 计算模块特征基因
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# 计算模块特征基因的不同
MEDiss = 1-cor(MEs);
# 聚类模块特征基因
METree = hclust(as.dist(MEDiss), method = "average");
# 绘制结果
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")

根据聚类结果，我们选择高度（height）为0.25的切线(对应于相关性为0.75)，来合并模块：

MEDissThres = 0.25
# 在聚类树中绘制切线
abline(h=MEDissThres, col = "red")
# 使用一个自动合并函数
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# 合并模块的颜色
mergedColors = merge$colors;
# 新的合并模块的特征基因
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()

3.4.7. 保存数据

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

# moduleColors重命名为合并模块的颜色名称
moduleColors = mergedColors
# 构造与颜色相对应的数字标签
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# 保存模块颜色和标签，以供后续部分使用
save(MEs, moduleLabels, moduleColors, geneTree, file = "networkConstruction-stepByStep.RData")

2.4. 分块法：网络构建和模块检测

在本教程中，我们使用了3600个测量探针的相对较小的数据集。然而，现代微阵列一次可测量高达50,000个探针表达水平。即使在大型服务器上，构建和分析具有如此大量节点的网络也具有计算上的挑战性。在WGCNA包中，提出了一个解决这种情况的方法。该方法允许用户使用如此大量的基因进行网络分析。

我们不是实际使用非常大的数据集，而是简单地假设硬件限制限制了可以一次分析到2000的基因数量。基本思想是使用两级聚类。首先，我们使用快速，计算上廉价且相对粗略的聚类方法将基因预聚类成接近并且不超过2000个基因的大小的块。然后，我们分别在每个块中执行完整的网络分析。最后，合并了特征基因高度相关的模块。

分块方法的优点是内存占用空间小得多（这是标准台式计算机上大数据集的主要问题），并且显着加快了计算速度。权衡是由于使用更简单的聚类来获得块，块可能不是最佳的，导致一些外围基因被分配到不同于完整网络分析中的模块。

我们现在假装即使我们在这里使用的相对较少数量的基因3600也太大了，我们运行分析的计算机不能在一个区块中处理超过2000个基因。 自动网络构建和模块检测功能blockwiseModules可以自动处理分块; 用户只需指定可容纳块的最大数量的基因：

2.4.1. 确定软阈值β

与2.2.1 确定软阈值β相同

最终，我们确定了软阈值β为为sft$powerEstimate中的值

2.4.2. 分块构建网络并识别模块

这里，我们在分块的基础上，采用2.2. 自动法构建网络并识别模块

bwnet = blockwiseModules(datExpr, maxBlockSize = 2000,
                         power = sft$powerEstimate, TOMType = "unsigned", 
                         minModuleSize = 30,
                         reassignThreshold = 0, mergeCutHeight = 0.25,
                         numericLabels = TRUE,
                         saveTOMs = TRUE,
                         saveTOMFileBase = "femaleMouseTOM-blockwise",
                         verbose = 3)

该功能的输出可能看起来有点神秘，但它易于使用。 例如，bwnet$colors包含模块赋值，bwnet$MEs包含模块的模块特征基因。

注意：

1. 函数blockwiseModules有许多参数，目前提供的参数可能不适合用户希望分析的特定数据集，对于想要根据自己的数据调整代码的，可以阅读R环境中包中提供的帮助文件，并尝试调整网络结构和模块检测参数，以获得更好的（生物学上更相关）分析结果。
2. 关于块大小的第二个警告。特别是，参数maxBlockSize告诉函数读者计算机可以处理的最大块的大小。在这个例子中，我们将最大块大小设置为2000来说明块分析及其结果，但对于大多数现代计算机来说，这个值是不必要的小。默认值为5000，适用于大多数现代桌面。如果服务器可以访问内存超过4 GB的大型工作站，则可以增加参数maxBlockSize。 1）16GB工作站最多可处理20000个探测器；2）32GB工作站应该可以处理30000个；3）4GB标准台式机或笔记本电脑可以处理多达8000-10000个探测器。通常，最好以尽可能少的块分析数据集。

2.4.3. 结果比较：与自动法

下面我们将该分析的结果与自动法的2.2.2. 构建网络并识别模块的结果进行比较，其中所有基因都在一个区块中进行分析。为了使比较更容易，我们重新标记了逐块模块标签，以便与单块模块有重大重叠的模块具有相同的标签：

2.4.3.1. 导入自动法分析的结果

# 导入自动法最终保存的结果
load(file = "networkConstruction-auto.RData");
# 模块标签重命名
bwLabels = matchLabels(bwnet$colors, moduleLabels);
# 将标签转换为用于绘图的颜色
bwModuleColors = labels2colors(bwLabels)
# 查看模块信息：
table(bwLabels)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
142 472 470 479 271 327 130 209 153 121 100 100 104 77 73 81 40 42 34 91 84
# 表示有20个模块，标记为1到20，标记为0的是所有模块外的基因。

2.4.3.2. 绘制分块法获得的结果

在每个分块中用于确定模块的的分层聚类树(tree)保存在bwnet$dendrograms[[1]]和bwnet$dendrograms[[2]]中。树形图可以与颜色分配一起显示，使用以下代码:

# 打开新绘图窗口
sizeGrWindow(6,6)
# 绘制block 1中的聚类树和模块
plotDendroAndColors(bwnet$dendrograms[[1]], bwModuleColors[bwnet$blockGenes[[1]]],
                    "Module colors", main = "Gene dendrogram and module colors in block 1", 
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
# 绘制block 2中的聚类树和模块
plotDendroAndColors(bwnet$dendrograms[[2]], bwModuleColors[bwnet$blockGenes[[2]]],
                    "Module colors", main = "Gene dendrogram and module colors in block 2", 
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

如果用户希望更改一些树的切割、模块成员和模块合并标准，WGCNA提供了recutBlockwiseTrees函数，该函数可以应用修改后的标准，而无需重新计算网络和集群树图，从而节省了大量时间。

2.4.3.3. 比较分块法和自动法结果

1）目视检查

我们绘制自动法的聚类树结果，并在下面分别绘制自动法的模块识别结果和分块法的模块识别结果，以比较两个方法：

sizeGrWindow(12,9)
plotDendroAndColors(geneTree,
                    cbind(moduleColors, bwModuleColors),
                    c("Single block", "2 blocks"),
                    main = "Single block gene dendrogram and module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

**结果：**目视检查证实，在分块法和自动法识别的模块有很好的一致性

2）模块特征基因比较

后续要分块法和自动法中相互对应的模块的模块特征基因是否极其相似的。我们首先计算两种方法中对应模块的特征基因：

# 自动法：模块特征基因
singleBlockMEs = moduleEigengenes(datExpr, moduleColors)$eigengenes;
# 分块法：模块特征基因
blockwiseMEs = moduleEigengenes(datExpr, bwModuleColors)$eigengenes;

接下来，我们将分块法和自动法中模块按名称进行匹配，并计算相应特征基因的相关性：

# 匹配两种方法的模块
single2blockwise = match(names(singleBlockMEs), names(blockwiseMEs))
# 比较模块中的特征基因相关性
signif(diag(cor(blockwiseMEs[, single2blockwise], singleBlockMEs)), 3)

结果：每个模块的相关性都非常接近于1，再次表明分块法和自动法分析结果非常相似。