学习WGCNA总结

来源：tiramisutes211,014

学习WGCNA总结-图片1
在转录组数据处理过程中我们经常会用到差异表达分析这一概念，通过比较不同处理或不同组织间基因表达量(FPKM)差异来寻找特异基因，但这前提是你的不同处理或不同组织样本较少，当不同处理或组织有较多样本，如40个，此时的两两比较有780组比较^_^,这根本不是我们想要的结果；

此时就需要WGCNA(weighted gene co-expression network analysis)将复杂的数据进行归纳整理。除了这种最常见的比较差异表达，我们还想知道在不同处理或不同组织间是否有些基因的表达存在内在的联系或相关性？WGCNA同样可以帮助我们预测基因间的相互作用关系。

WGCNA术语

权重(weghted)

学习WGCNA总结-图片2

Module

模块(module)：表达模式相似的基因分为一类，这样的一类基因成为模块；

Eigengene

Eigengene（eigen- +‎ gene）：基因和样本构成的矩阵，https://en.wiktionary.org/wiki/eigengene；

Adjacency Matrix

邻近矩阵：是图的一种存储形式，用一个一维数组存放图中所有顶点数据；用一个二维数组存放顶点间关系（边或弧）的数据，这个二维数组称为邻接矩阵；
学习WGCNA总结-图片3

Topological Overlap Matrix(TOM)

学习WGCNA总结-图片4

整体思路

先对数据进行处理→分层聚类→表达模式相似的基因组成模块→研究某一个模块中相关基因的功能富集(GO,KEGG)，各个模块与样本表型数据间的相关性，各个模块与样本本身间的相关性(没有表型数据的情况，如不同组织)→具体到特定模块后分析其所包含基因间的相互作用网络关系，并找出其中的关键基因。

分析构建的网络寻找以下有用信息

 

  • 这类处于调控网络中心的基因称为核心基因（hub gene），这类基因通常是转录因子等关键的调控因子，是值得我们优先深入分析和挖掘的对象。
  • 在网络中，被调控线连接的基因，其表达模式是相似的。那么它们潜在有相似的功能。所以，在这个网络中，如果线条一端的基因功能是已知的，那么就可以预测线条另一端的功能未知的基因也有相似的功能。

 

R脚本

输入数据为RNA-seq不同处理或组织所有样本的FPKM值组成的矩阵，切记含有 0 的要去掉；

  1. setwd("F:/WGCNA")
  2. library(WGCNA)
  3. options(stringsAsFactors = FALSE)
  4. enableWGCNAThreads()
  5. #1. 数据读入，处理和保存
  6. fpkm<- read.csv("trans_counts.counts.matrix.TMM_normalized.FPKM.nozero.csv")
  7. head(fpkm)
  8. dim(fpkm)
  9. names(fpkm)
  10. datExpr0=as.data.frame(t(fpkm[,-c(1)]));
  11. names(datExpr0)=fpkm$trans;
  12. rownames(datExpr0)=names(fpkm)[-c(1)];
  13. #data<-log10(date[,-1]+0.01)
  14. gsg = goodSamplesGenes(datExpr0, verbose = 3);
  15. gsg$allOK
  16. sampleTree = hclust(dist(datExpr0), method = "average")
  17. #sizeGrWindow(12,9)
  18. par(cex = 0.6)
  19. par(mar = c(0,4,2,0))
  20. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
  21.      cex.axis = 1.5, cex.main = 2)
  22. abline(h = 80000, col = "red");
  23. clust = cutreeStatic(sampleTree, cutHeight = 80000, minSize = 10)
  24. table(clust)
  25. keepSamples = (clust==1)
  26. datExpr = datExpr0[keepSamples, ]
  27. nGenes = ncol(datExpr)
  28. nSamples = nrow(datExpr)
  29. save(datExpr, file = "AS-green-FPKM-01-dataInput.RData")
  30. #2. 选择合适的阀值
  31. powers = c(c(1:10), seq(from = 12, to=20, by=2))
  32. # Call the network topology analysis function
  33. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
  34. # Plot the results:
  35. ##sizeGrWindow(9, 5)
  36. par(mfrow = c(1,2));
  37. cex1 = 0.9;
  38. # Scale-free topology fit index as a function of the soft-thresholding power
  39. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  40.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
  41.      main = paste("Scale independence"));
  42. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  43.      labels=powers,cex=cex1,col="red");
  44. # this line corresponds to using an R^2 cut-off of h
  45. abline(h=0.90,col="red")
  46. # Mean connectivity as a function of the soft-thresholding power
  47. plot(sft$fitIndices[,1], sft$fitIndices[,5],
  48.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
  49.      main = paste("Mean connectivity"))
  50. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
  51.  
  52. #=====================================================================================
  53. # 网络构建有两种方法，One-step和Step-by-step；
  54. #  第一种：一步法进行网络构建
  55. #=====================================================================================
  56.  
  57. #3. 一步法网络构建：One-step network construction and module detection
  58. net = blockwiseModules(datExpr, power = 6, maxBlockSize = 6000,
  59.                        TOMType = "unsigned", minModuleSize = 30,
  60.                        reassignThreshold = 0, mergeCutHeight = 0.25,
  61.                        numericLabels = TRUE, pamRespectsDendro = FALSE,
  62.                        saveTOMs = TRUE,
  63.                        saveTOMFileBase = "AS-green-FPKM-TOM",
  64.                        verbose = 3)
  65. table(net$colors)
  66. #4. 绘画结果展示
  67. # open a graphics window
  68. #sizeGrWindow(12, 9)
  69. # Convert labels to colors for plotting
  70. mergedColors = labels2colors(net$colors)
  71. # Plot the dendrogram and the module colors underneath
  72. plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
  73.                     "Module colors",
  74.                     dendroLabels = FALSE, hang = 0.03,
  75.                     addGuide = TRUE, guideHang = 0.05)
  76. #5.结果保存
  77. moduleLabels = net$colors
  78. moduleColors = labels2colors(net$colors)
  79. table(moduleColors)
  80. MEs = net$MEs;
  81. geneTree = net$dendrograms[[1]];
  82. save(MEs, moduleLabels, moduleColors, geneTree,
  83.      file = "AS-green-FPKM-02-networkConstruction-auto.RData")
  84. #6. 导出网络到Cytoscape
  85. # Recalculate topological overlap if needed
  86. TOM = TOMsimilarityFromExpr(datExpr, power = 6);
  87. # Read in the annotation file
  88. # annot = read.csv(file = "GeneAnnotation.csv");
  89. # Select modules需要修改，选择需要导出的模块颜色
  90. modules = c("turquoise", "blue");
  91. # Select module probes选择模块探测
  92. probes = names(datExpr)
  93. inModule = is.finite(match(moduleColors, modules));
  94. modProbes = probes[inModule];
  95. #modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
  96. # Select the corresponding Topological Overlap
  97. modTOM = TOM[inModule, inModule];
  98. dimnames(modTOM) = list(modProbes, modProbes)
  99. # Export the network into edge and node list files Cytoscape can read
  100. cyt = exportNetworkToCytoscape(modTOM,
  101.                                edgeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
  102.                                nodeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
  103.                                weighted = TRUE,
  104.                                threshold = 0.02,
  105.                                nodeNames = modProbes,
  106.                                #altNodeNames = modGenes,
  107.                                nodeAttr = moduleColors[inModule]);
  108. #=====================================================================================
  109. #  分析网络可视化，用heatmap可视化权重网络，heatmap每一行或列对应一个基因，颜色越深表示有较高的邻近
  110. #=====================================================================================
  111. options(stringsAsFactors = FALSE);
  112. lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
  113. lnames
  114. lnames = load(file = "AS-green-FPKM-02-networkConstruction-auto.RData");
  115. lnames
  116. nGenes = ncol(datExpr)
  117. nSamples = nrow(datExpr)
  118. #1. 可视化全部基因网络
  119. # Calculate topological overlap anew: this could be done more efficiently by saving the TOM
  120. # calculated during module detection, but let us do it again here.
  121. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
  122. # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
  123. plotTOM = dissTOM^7;
  124. # Set diagonal to NA for a nicer plot
  125. diag(plotTOM) = NA;
  126. # Call the plot function
  127. #sizeGrWindow(9,9)
  128. TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
  129. #随便选取1000个基因来可视化
  130. nSelect = 1000
  131. # For reproducibility, we set the random seed
  132. set.seed(10);
  133. select = sample(nGenes, size = nSelect);
  134. selectTOM = dissTOM[select, select];
  135. # There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
  136. selectTree = hclust(as.dist(selectTOM), method = "average")
  137. selectColors = moduleColors[select];
  138. # Open a graphical window
  139. #sizeGrWindow(9,9)
  140. # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
  141. # the color palette; setting the diagonal to NA also improves the clarity of the plot
  142. plotDiss = selectTOM^7;
  143. diag(plotDiss) = NA;
  144. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
  145. #=====================================================================================
  146. #  第二种：一步步的进行网络构建
  147. #=====================================================================================
  148. ###################Step-by-step network construction and module detection
  149. #2.选择合适的阀值，同上
  150. #3. 网络构建：(1) Co-expression similarity and adjacency
  151. softPower = 6;
  152. adjacency = adjacency(datExpr, power = softPower);
  153. #(2) 邻近矩阵到拓扑矩阵的转换，Turn adjacency into topological overlap
  154. TOM = TOMsimilarity(adjacency);
  155. dissTOM = 1-TOM
  156. # (3) 聚类拓扑矩阵
  157. #Call the hierarchical clustering function
  158. geneTree = hclust(as.dist(dissTOM), method = "average");
  159. # Plot the resulting clustering tree (dendrogram)
  160. #sizeGrWindow(12,9)
  161. plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
  162.      labels = FALSE, hang = 0.04);
  163. #(4) 聚类分支的休整dynamicTreeCut
  164. # We like large modules, so we set the minimum module size relatively high:
  165. minModuleSize = 30;
  166. # Module identification using dynamic tree cut:
  167. dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
  168.                             deepSplit = 2, pamRespectsDendro = FALSE,
  169.                             minClusterSize = minModuleSize);
  170. table(dynamicMods)
  171. #4. 绘画结果展示
  172. # Convert numeric lables into colors
  173. dynamicColors = labels2colors(dynamicMods)
  174. table(dynamicColors)
  175. # Plot the dendrogram and colors underneath
  176. #sizeGrWindow(8,6)
  177. plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
  178.                     dendroLabels = FALSE, hang = 0.03,
  179.                     addGuide = TRUE, guideHang = 0.05,
  180.                     main = "Gene dendrogram and module colors")
  181. #5. 聚类结果相似模块的融合，Merging of modules whose expression profiles are very similar
  182. #在聚类树中每一leaf是一个短线，代表一个基因，
  183. #不同分之间靠的越近表示有高的共表达基因，将共表达极其相似的modules进行融合
  184. # Calculate eigengenes
  185. MEList = moduleEigengenes(datExpr, colors = dynamicColors)
  186. MEs = MEList$eigengenes
  187. # Calculate dissimilarity of module eigengenes
  188. MEDiss = 1-cor(MEs);
  189. # Cluster module eigengenes
  190. METree = hclust(as.dist(MEDiss), method = "average");
  191. # Plot the result
  192. #sizeGrWindow(7, 6)
  193. plot(METree, main = "Clustering of module eigengenes",
  194.      xlab = "", sub = "")
  195. #选择有75%相关性的进行融合
  196. MEDissThres = 0.25
  197. # Plot the cut line into the dendrogram
  198. abline(h=MEDissThres, col = "red")
  199. # Call an automatic merging function
  200. merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
  201. # The merged module colors
  202. mergedColors = merge$colors;
  203. # Eigengenes of the new merged modules:
  204. mergedMEs = merge$newMEs;
  205. #绘制融合前(Dynamic Tree Cut)和融合后(Merged dynamic)的聚类图
  206. #sizeGrWindow(12, 9)
  207. #pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
  208. plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
  209.                     c("Dynamic Tree Cut", "Merged dynamic"),
  210.                     dendroLabels = FALSE, hang = 0.03,
  211.                     addGuide = TRUE, guideHang = 0.05)
  212. #dev.off()
  213. # 只是绘制融合后聚类图
  214. plotDendroAndColors(geneTree,mergedColors,"Merged dynamic",
  215.                     dendroLabels = FALSE, hang = 0.03,
  216.                     addGuide = TRUE, guideHang = 0.05)
  217. #5.结果保存
  218. # Rename to moduleColors
  219. moduleColors = mergedColors
  220. # Construct numerical labels corresponding to the colors
  221. colorOrder = c("grey", standardColors(50));
  222. moduleLabels = match(moduleColors, colorOrder)-1;
  223. MEs = mergedMEs;
  224. # Save module colors and labels for use in subsequent parts
  225. save(MEs, moduleLabels, moduleColors, geneTree, file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData")
  226. #6. 导出网络到Cytoscape
  227. # Recalculate topological overlap if needed
  228. TOM = TOMsimilarityFromExpr(datExpr, power = 6);
  229. # Read in the annotation file
  230. # annot = read.csv(file = "GeneAnnotation.csv");
  231. # Select modules需要修改
  232. modules = c("brown", "red");
  233. # Select module probes
  234. probes = names(datExpr)
  235. inModule = is.finite(match(moduleColors, modules));
  236. modProbes = probes[inModule];
  237. #modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
  238. # Select the corresponding Topological Overlap
  239. modTOM = TOM[inModule, inModule];
  240. dimnames(modTOM) = list(modProbes, modProbes)
  241. # Export the network into edge and node list files Cytoscape can read
  242. cyt = exportNetworkToCytoscape(modTOM,
  243.                                edgeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
  244.                                nodeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
  245.                                weighted = TRUE,
  246.                                threshold = 0.02,
  247.                                nodeNames = modProbes,
  248.                                #altNodeNames = modGenes,
  249.                                nodeAttr = moduleColors[inModule]);
  250. #=====================================================================================
  251. #  分析网络可视化，用heatmap可视化权重网络，heatmap每一行或列对应一个基因，颜色越深表示有较高的邻近
  252. #=====================================================================================
  253. options(stringsAsFactors = FALSE);
  254. lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
  255. lnames
  256. lnames = load(file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData");
  257. lnames
  258. nGenes = ncol(datExpr)
  259. nSamples = nrow(datExpr)
  260. #1. 可视化全部基因网络
  261. # Calculate topological overlap anew: this could be done more efficiently by saving the TOM
  262. # calculated during module detection, but let us do it again here.
  263. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
  264. # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
  265. plotTOM = dissTOM^7;
  266. # Set diagonal to NA for a nicer plot
  267. diag(plotTOM) = NA;
  268. # Call the plot function
  269. #sizeGrWindow(9,9)
  270. TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
  271. #随便选取1000个基因来可视化
  272. nSelect = 1000
  273. # For reproducibility, we set the random seed
  274. set.seed(10);
  275. select = sample(nGenes, size = nSelect);
  276. selectTOM = dissTOM[select, select];
  277. # There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
  278. selectTree = hclust(as.dist(selectTOM), method = "average")
  279. selectColors = moduleColors[select];
  280. # Open a graphical window
  281. #sizeGrWindow(9,9)
  282. # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
  283. # the color palette; setting the diagonal to NA also improves the clarity of the plot
  284. plotDiss = selectTOM^7;
  285. diag(plotDiss) = NA;
  286. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
  287. #此处画的是根据基因间表达量进行聚类所得到的各模块间的相关性图
  288. MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
  289. MET = orderMEs(MEs)
  290. sizeGrWindow(7, 6) 
  291. plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)

 

部分结果图简单解释

学习WGCNA总结-图片5

Cytoscape生成网络图

只需要第二个edges文件就能构建网络图。导入该文件后，在软件的导入设置中，将第一列设置为fromNode，第二列设置为toNode，最后把第三列设为网络关系属性，完成设置，便可生成网络图了。

WGCNA样本要求

由于WGCNA是基于相关系数的表达调控网络分析方法。当样本数过低的时候，相关系数的计算是不可靠的，得到的调控网络价值不大。所以，我们推荐的样本数如下：

 

  • 当独立样本数≥8（非重复样本）时，可以考虑基于Pearson相关系数的WGCNA共表达网络的方法（效果看实际情况而定）；
  • 当样本数≥15（可以包含生物学重复）时，WGCNA方法会有更好的效果。
  • 当样品数＜8时，不建议进行该项分析。

 

报错暨解决办法

按照WGCNA手册第五步Network visualization using WGCNA functions时报错如下：

  1. > TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
  2.  
  3. Error in .heatmap(as.matrix(dissim), Rowv = as.dendrogram(dendro, hang = 0.1),  :
  4. row dendrogram ordering gave index of wrong length

看到row dendrogram ordering gave index of wrong length这句报错内容，分别察看plotTOM, geneTree, moduleColors这三个变量length;

  1. > dim(plotTOM)
  2. > geneTree
  3. > moduleColors

果然，三者的length不同，发现geneTree少了一些，往回找geneTree来源 geneTree = net$dendrograms[[1]]，net来源于网络构建过程：

  1. net = blockwiseModules(datExpr, power = 6,
  2. TOMType = "unsigned", minModuleSize = 30,
  3. reassignThreshold = 0, mergeCutHeight = 0.25,
  4. numericLabels = TRUE, pamRespectsDendro = FALSE,
  5. saveTOMs = TRUE,
  6. saveTOMFileBase = "femaleMouseTOM",
  7. verbose = 3)

所以，这是问题所在，继续察看文档发现blockwiseModules函数默认最大maxBlockSize=5000，而我们的数据超过了这个值，所以函数自动做了拆分处理，而解决办法也很简单，设置maxBlockSize参数大于我们的值即可。

评论  2  访客  2
    • fengxj 0
      1F
      回复

      非常有用，谢谢！

      • hope 1
        2F
        回复

        原文出处：http://tiramisutes.github.io/2016/09/14/WGCNA.html

      发表评论

      匿名网友

      拖动滑块以完成验证
      加载失败