在前面的帖子中介绍了数据的导入和清洗,网络构建的两种方法,模块与性状的关联,这篇文章将介绍如果进行模块可视化。
WGCNA(1):R包安装及数据导入清洗 - 简书 (jianshu.com)
WGCNA(2a):一步法完成网络构建和模块检测 - 简书 (jianshu.com)
WGCNA(2b):分步法完成网络构建和模块检测 - 简书 (jianshu.com)
WGCNA(3):基因模块与性状关联识别重要基因 - 简书 (jianshu.com)
1. 准备工作
导入前期数据,这里我选择了分步法构建网络的结果,大家可以根据自己的数据选择使用哪种方法构建网络。
# 设置工作目录
> setwd("D:/RNA-seq/WGCNA/mad0.3/cor 0.25")
# 载入WGCNA包
> library('WGCNA')
# 允许R语言以最大线程运行
> options(stringsAsFactors = FALSE)
> allowWGCNAThreads()
Allowing multi-threading with up to 4 threads.
# 载入第一步中的表达量和表型值
> lnames = load(file = "WGCNA0.3-dataInput.RData")
> lnames
# 载入第二步的网络数据
> lnames = load(file = "networkConstruction-stepByStep.RData")
> lnames
> nGenes = ncol(datExpr)
> nSamples = nrow(datExpr)
2. 利用R进行网络可视化
2.1 可视化TOM矩阵
以下代码只适用于一步法或分步法构建的网络,如果是 block-wise方法,则需要各个部分分开绘图。
# 模块检测时的计算,重新算一次
> dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6)
# 转换dissTOM,方便在热图中显示
> plotTOM = dissTOM^7
# 将对角线设置为NA
> diag(plotTOM) = NA
# 绘图
> pdf("Network heatmap plot.pdf", width=9, height=9)
> TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot_all genes")
> dev.off()
部分可视化TOM矩阵
上面这张图画起来非常非常慢,又占时间,又占内存,而且意义也不是很大,感觉是在凑图,4w个基因,服务器上要画半个多小时,所以,如果实在想要这张图凑数据,其实选部分基因就可以了(Fig.2)。
> nSelect = 1000
# For reproducibility, we set the random seed
> set.seed(10)
> select = sample(nGenes, size = nSelect)
> selectTOM = dissTOM[select, select]
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
> selectTree = hclust(as.dist(selectTOM), method = "average")
> selectColors = moduleColors[select]
# Open a graphical window
> pdf("Network heatmap plot_selected genes.pdf", width = 9, height = 9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
> plotDiss = selectTOM^7
> diag(plotDiss) = NA
> TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot_selected genes")
> dev.off()
这个热图的配色可以改,详情见:https://cloud.tencent.com/developer/article/1628251?from=article.detail.1516749
2.2 eigengenes可视化
在eigengenes中加入表型
# 提取性状weight
> weight = as.data.frame(datTraits$weight_g)
> names(weight) = "weight"
# 在eigengenes模块中加入性状
> MET = orderMEs(cbind(MEs, weight))
# 绘制eigengenes和性状之间的关系图
> pdf("Visualization of the eigengene network representing the relationships among the modules and the clinical trait weight.pdf",width = 7, height=6)
#一页多图,一行2列
> par(mfrow = c(1,2))
#字号
> cex1 = 0.9
> plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE)
> plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)
> dev.off()
图中可以看出weight性状和blue、brown还有red模块聚类在一起。