GEO数据挖掘系列，第21篇学习笔记：WGCNA代码实操之正式阶段
为了无缝衔接上一篇学习笔记，该文中的序号将接着上一篇来标注
3.WGCNA正式开始
以下代码将不用再做太多修改
3.1 软阈值的筛选
情况一：能得到推荐值
设置一系列软阈值，范围是1-30之间，后面的数没必要全部画，就seq一下
powers = c(1:10, seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
sft$powerEstimate#如果自己没有合适的阈值，此行代码将算出一个值
sft$powerEstimate
的结果就是推荐的软阈值，拿到了可以直接用，无视下面的图4
情况二：没能得到推荐值
有的数据
sft$powerEstimate
的结果是1或NA，极少情况用1作为阈值，而NA就是没有推荐的意思
那就要运行以下代码，看画出来的图4，选拐点了。以下代码不需要改动，唯一可能要改动的就是0.9这个数值
cex1 = 0.9
png(file = "4.Soft threshold.png", width = 2000, height = 1500,res = 300)
par(mfrow = c(1,2))
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=cex1,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")
dev.off()
图4 Soft threshold
上一篇学习笔记也提到了，cex1一般设置为0.9，不太合适（就是大多数软阈值对应的纵坐标都达不到0.9）时，可以设置为0.8或者0.85。一般不能再低了
3.2 一步构建网络
如果前面没有得到推荐的软阈值，这里的
power
就要手动输入一个值。根据上面的图4，选择左图纵坐标第一个达到上面设置的cex1值的软阈值
如果从图4中都找不到一个合适的软阈值，根据经验，这个WGCNA就做不了了；或者就得跑到前面重新调整表达矩阵里纳入的基因了
power = sft$powerEstimate
power
net = blockwiseModules(datExpr, power = power,
TOMType = "unsigned",
minModuleSize = 30,
reassignThreshold = 0,
mergeCutHeight = 0.25,
deepSplit = 2 ,
numericLabels = TRUE,
pamRespectsDendro = FALSE,
saveTOMs = TRUE,
saveTOMFileBase = "testTOM",
verbose = 3)
后面的结果如果不理想，比如划分的模块太少，或者青色、灰色的模块占到大多数，其他颜色都很少，或者颜色模块聚类是凌乱的，可以倒回来调整几个参数：
deepSplit
默认2，调整划分模块的敏感度，值越大，越敏感，得到的模块就越多
minModuleSize
默认30，参数设置最小模块的基因数，值越小，小的模块就会被保留下来
mergeCutHeight
默认0.25，设置合并相似性模块的距离，值越小，就越不容易被合并，保留下来的模块就越多
更详细设置可参考https://zhuanlan.zhihu.com/p/34697561
可以调整试试，但是主要还得看选择的基因是否给力，不给力的话调整了也只能相对好一点点
不是很有必要去尝试分步法构建网络，得到的结果一样，可以调整的参数上面也都有
table(net$colors)#此处展示得到了多少模块，每个模块里面有多少基因
mergedColors = labels2colors(net$colors)
png(file = "5.DendroAndColors.png", width = 2000, height = 1200,res = 300)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
图5 DendroAndColors
好的结果就是每个颜色都差不多在一起（青色和灰色不在考虑范围内），以及青色和灰色的基因不要太多
因为灰色代表没有合适的聚类，青色是基因数量的模块，比如你输入5000个基因，其中3000个都属于青色，剩下的模块基因数量太少，就很难受了
3.3 保存每个模块对应的基因
这里把每个模块对应的基因存为了Rdata，用于数据挖掘下一步需求，提取基因。比如你需要某模块的基因与差异基因取交集等
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]
gm = data.frame(net$colors)
gm$color = moduleColors
head(gm)
genes = split(rownames(gm),gm$color)
save(genes,file = "genes.Rdata")
3.4 模块与表型的相关性
计算基因与表型的相关性矩阵
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
#热图
png(file = "6.labeledHeatmap.png", width = 2000, height = 2000,res = 300)
# 设置热图上的文字（两行数字：第一行是模块与各种表型的相关系数；
# 第二行是p值）
# signif 取有效数字
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# 然后对moduleTraitCor画热图
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = blueWhiteRed (50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
dev.off()
我们希望找到和每个表型相关性较强的模块，正负相关都可。相关系数越大越好，如果能有个0.8，那结论就比较稳。
没有的话0.6或0.7几也行，但再小就不能用了
图6 labeledHeatmap
3.5. GS与MM
GS代表模块里的每个基因与形状的相关性
MM代表每个基因和所在模块之间的相关性，表示是否与模块的趋势一致
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
下面的代码中，第几列的表型是最关心的，下面的i就设置为几
与关心的表型相关性最高的模块赋值给下面的module
# i和module是对应的，第2列它的模块颜色即青色的
i = 2
#module = "pink"
module = "turquoise"
assign(colnames(traitData)[i],traitData[i])
instrait = eval(parse(text = colnames(traitData)[i]))
geneTraitSignificance = as.data.frame(cor(datExpr, instrait, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(instrait), sep="")
names(GSPvalue) = paste("p.GS.", names(instrait), sep="")
png(file = paste0("7.MM-GS-scatterplot.png"), width = 2000, height = 2000,res = 300)
column = match(module, modNames) #找到目标模块所在列
moduleGenes = moduleColors==module #找到模块基因所在行
par(mfrow = c(1,1))
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()
图7 MM-GS-scatterplot
我们希望看到强悍的正相关。数据挖掘里面，可以把整个模块的基因拿出来和别的步骤结果取交集了
也可以通过这里找到GS和MM都大的基因
f = data.frame(MM = abs(geneModuleMembership[moduleGenes, column]),
GS = abs(geneTraitSignificance[moduleGenes, 1]))
