1.导入数据
导入第一步生成的网络关系和第二步的性状之间的关系。
# 加载 R 包
library(WGCNA)
# 导入数据
load(file = "wgcna-network.Rdata")
load(file = "trait_analysis.RData")
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
2. 基因共表达网络可视化
2.1. 使用所有基因绘制
{
# 计算 dissTOM, dissTOM = 1 - TOM
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate);
# 取 7 次方,仅为展示更显著
plotTOM = dissTOM^7
# 绘图
diag(plotTOM) = NA;
sizeGrWindow(9,9)
geneTree = net$dendrograms[[1]]
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
官网推荐了用随机400条基因作图,全部用于作图的话没必要,而且时间很慢。
2.2 随机选择 400 条绘图
nSelect = 400
set.seed(10)
select = sample(nGenes, size = nSelect)
selectTOM = dissTOM[select, select]
# 重新聚类
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select]
# Open a graphical window
sizeGrWindow(9,9)
plotDiss = selectTOM^7
# 绘图
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
}
TOMplot
3. 模块特征向量网络可视化
{
# 重新计算MEs
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5)
par(cex = 0.9)
plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
## 当然也可以单独绘制
# Plot the dendrogram
sizeGrWindow(6,6)
par(cex = 1.0)
plotEigengeneNetworks(MEs, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)
}
Eigengene adjacency heatmap
4. 生成 Cytoscape 的输入文件
{
# 重新计算 TOM
TOM = TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate);
# 要可视化的模块
modules = c("blue")
# 要可视化的基因
probes = colnames(datExpr)
inModule = is.finite(match(moduleColors, modules));
modProbes = probes[inModule]
# 候选基因的 TOM
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
weighted = TRUE,
threshold = 0.02,
nodeNames = modProbes,
nodeAttr = moduleColors[inModule])
}
网友评论