之前我们完成了WGCNA
输入数据的清洗,网络构建和模块识别。😘
但这里我们的示例数据内所含有的基因其实是很少的,而在实际情况中,一个简单的测序可能就要包含上万个基因,这对大家的电脑无疑是不小的压力。🤒
在WGCNA
的包内其实也提供了解决方案,基本思想是分级聚类
。🧐
1️⃣ 首先,我们使用快速但相对粗糙的聚类方法,用于将基因预聚类成大小接近的模块
,且不超过你所设定的基因最大值
。😂
2️⃣ 然后我们分别在每个模块
中执行完整的网络分析。🤠
3️⃣ 最后,合并特征基因高度相关的模块。😏
rm(list = ls())
library(WGCNA)
library(tidyverse)
load("FemaleLiver-01-dataInput.RData")
首先我们还是要和之前一样进行soft thresholding power β
的计算。🤒
powers <- c(c(1:10), seq(from = 12, to=20, by=2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
显然,我们的结果和之前是一样的,6
。😜
sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.9
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")
abline(h=0.90,col="red")
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")
这里我们就要设置每个block
的最大size
是多少了,分次计算,再合并,以此减少内存的负担。🤒
bwnet <- blockwiseModules(
datExpr, maxBlockSize = 2000,
power = 6, TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = T,
saveTOMs = T,
saveTOMFileBase = "femaleMouseTOM-blockwise",
verbose = 3)
table(bwnet$colors)
这里我们再对比一下结果,看看2种方法得出结果的区别。🤪
这里我们把之前不分次计算的结果载入进来,后面会用到label
和colors
。😷
load(file = "FemaleLiver-02-networkConstruction-auto.RData")
把colors
和label
匹配起来,嘿嘿。🤨
bwLabels <- matchLabels(bwnet$colors, moduleLabels);
bwModuleColors <- labels2colors(bwLabels)
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 = F, hang = 0.03,
addGuide = T, 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 = F, hang = 0.03,
addGuide = T, guideHang = 0.05)
我们以可视化的形式对比一下,分割出来的模块差异不大。😂
sizeGrWindow(12,9)
plotDendroAndColors(geneTree,
cbind(moduleColors, bwModuleColors),
c("Single block", "2 blocks"),
main = "Single block gene dendrogram and module colors",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05)
这里我们提取一下2种方法得到的module eigengenes
。🥳
singleBlockMEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
blockwiseMEs <- moduleEigengenes(datExpr, bwModuleColors)$eigengenes
match
之后看一下结果,嘿嘿。🫶
高度一致,所以这种blockwise
的方法,请放心食用吧,各位。😉
single2blockwise <- match(names(singleBlockMEs), names(blockwiseMEs))
signif(diag(cor(blockwiseMEs[, single2blockwise], singleBlockMEs)), 3)
📍
Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559
最后祝大家早日不卷!~