上期我们完成了WGCNA
输入数据的清洗,然后进行了样本的聚类与异常值的剔除,总体来说是非常简单的。😘 这期我们继续完成WGCNA
分析的第二步,网络构建和模块识别。🤒
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)
这个时候我们可以看到最佳的soft thresholding power β
为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")
这里我们只需要一个函数就可以完成网络构建和模块识别,即blockwiseModules
,不过参数比较多,大家需要去好好理解,这里不做太多介绍了,大家看帮助文档吧。🧐
这里会有一个maxBlockSize
的参数,大家可以根据自己的电脑配置进行选择,16 GB
内存最大值为20000
,32 GB
内存最大值为30000
。🤨
net <- blockwiseModules(datExpr, power = 6,
TOMType = "unsigned",
minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = T, pamRespectsDendro = F,
saveTOMs = T,
saveTOMFileBase = "femaleMouseTOM",
verbose = 3)
table(net$colors)
sizeGrWindow(12,9)
mergedColors <- labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors", dendroLabels = F, hang = 0.03, addGuide = T,
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")
这里我们把之前算好的power
应用一下。😚
softPower <- 6
adjacency <- adjacency(datExpr, power = softPower)
这里我们需要把邻接矩阵
转换为Topological Overlap Matrix
(TOM矩阵
),最大限度地减少噪音和假相关性。🧐
TOM <- TOMsimilarity(adjacency);
dissTOM <- 1-TOM
geneTree <- hclust(as.dist(dissTOM), method = "average")
plot(geneTree,
xlab="", sub="",
main = "Gene clustering on TOM-based dissimilarity",
labels = F, hang = 0.04)
我们一般喜欢基因多一点的模块,这里将最小的模块设置为30
。😏
补充一下 !这里0
模块是指未能归类到任意模块的基因。😷
minModuleSize <- 30
dynamicMods <- cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = F,
minClusterSize = minModuleSize)
table(dynamicMods)
上面的模块为数字,我们需要把它转成颜色进行可视化。🤨
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)
sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors,
"Dynamic Tree Cut",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05,
main = "Gene dendrogram and module colors")
MEList <- moduleEigengenes (datExpr, colors = dynamicColors)
MEs <- MEList$eigengenes
MEDiss = 1-cor(MEs);
METree <- hclust(as.dist(MEDiss), method = "average")
这里我们将相关性在0.75
以上的模块合并在一起,当然你也可以按照你的要求来调整。🫶
sizeGrWindow(7,6)
plot(METree, main = "Clustering of module eigengenes",xlab = "", sub = "")
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)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = F, hang = 0.03, addGuide = T, guideHang = 0.05)
moduleColors <- mergedColors
colorOrder <- c("grey", standardColors(50))
moduleLabels <- match(moduleColors, colorOrder)-1
MEs <- mergedMEs
我们save
一下这个data
吧,后面会再用到。😉
save(MEs, moduleLabels, moduleColors, geneTree,
file = "FemaleLiver-02-networkConstruction-stepByStep.RData")
📍
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
最后祝大家早日不卷!~