1.導(dǎo)入數(shù)據(jù)

導(dǎo)入第一步生成的網(wǎng)絡(luò)關(guān)系和第二步的性狀之間的關(guān)系。

# 加載 R 包
library(WGCNA)

# 導(dǎo)入數(shù)據(jù)
load(file = "wgcna-network.Rdata")
load(file = "trait_analysis.RData")
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

2. 基因共表達(dá)網(wǎng)絡(luò)可視化

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")

官網(wǎng)推薦了用隨機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. 模塊特征向量網(wǎng)絡(luò)可視化

{
  # 重新計算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)
 
## 當(dāng)然也可以單獨繪制 
  # 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])
}