• GEO生信数据挖掘（九）WGCNA分析

    第六节，我们使用结核病基因数据，做了一个数据预处理的实操案例。例子中结核类型，包括结核，潜隐进展，对照和潜隐，四个类别。第七节延续上个数据，进行了差异分析。 第八节对差异基因进行富集分析。本节进行WGCNA分析。

    更新版本：GEO生信数据挖掘（九）肺结核数据-差异分析-WGCNA分析（900行代码整理注释更新版本）

    目录

    加载数据，进行聚类

    初次聚类观察

    自己定义红线位置，进行切割划分

    载入性状数据

    增加形状信息后，再次聚类

    网络构建

    选取soft-thresholding powers

    基于tom的差异的基因聚类，绘制聚类树

    根据聚类情况，设置颜色

    计算eigengenes

    模块的自动合并

    模块与临床形状的关系热图 （关键数据）

    红色模块样本表达情况（相关性大）

    产生了很多数据（各个模块的和临床性状）

    后续挖掘核心基因时，需要用到Cytoscape，生成绘图所需要的数据

    加载数据，进行聚类

    1. library(WGCNA)
    2. #读取目录名称，方便复制粘贴
    3. dir()
    4. #加载数据
    5. load('DEG_TB_LTBI_step13.Rdata')
    6.  
    7.  
    8. #这里行为样品名，列为基因名
    9. ################################样品聚类#################### 
    10. datExpr = t(dataset_TB_LTBI_DEG)
    11. #初次聚类
    12. sampleTree = hclust(dist(datExpr), method = "average")
    13. # Plot the sample tree: Open a graphic output window of size 20 by 15 inches
    14. # The user should change the dimensions if the window is too large or too small.
    15. sizeGrWindow(12,9)
    16. #pdf(file='sampleCluestering.pdf',width = 12,height = 9)
    17. par(cex=0.6)
    18. par(mar=c(0,4,2,0))
    19. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
    20.      cex.axis = 1.5, cex.main = 2)
    21.  
    22. #结果图片自己导出PDF，文件名=1_sampleClustering.pdf
    23.  
    24.  
    25.  
    26. ### Plot a line to show the cut
    27. abline(h = 87, col = "red")##剪切高度不确定，故无红线
    28.  
    29. dev.off()

    初次聚类观察

    自己定义红线位置，进行切割划分

    本例发现右侧有些样本孤立，适合被剔除，设置红线87切割。

    左侧也被切成两块，需要做处理，保留。

    1.  
    2. ### Determine cluster under the line
    3. clust = cutreeStatic(sampleTree, cutHeight = 87, minSize = 10)
    4. table(clust)
    5. #clust
    6. #0  1  2 
    7. #5 57 40
    8.  
    9. ### clust 1 contains the samples we want to keep.
    10. keepSamples = (clust==1|clust==2)
    11. datExpr0 = datExpr[keepSamples, ]
    12. dim(datExpr0) #[1]   97 2813
    13. #保存数据
    14. save(datExpr0,file='datExpr0_cluster_filter.Rdata')

    载入性状数据

    匹配样本名称，性状数据与表达数据保证一致

    1.  
    2. #################### 载入性状数据###########################
    3. #加载自己的性状数据
    4. load('design_TB_LTBI.Rdata')
    5. traitData=design
    6. #Loading clinical trait data
    7. #traitData = read.table("trait_D.txt",row.names=1,header=T,comment.char = "",check.names=F)########trait file name can be changed######性状数据文件名，根据实际修改，如果工作路径不是实际性状数据路径，需要添加正确的数据路径
    8. dim(traitData)
    9. #names(traitData)
    10. # remove columns that hold information we do not need.
    11. #allTraits = traitData
    12. dim(traitData)
    13. names(traitData)
    14.  
    15. # Form a data frame analogous to expression data that will hold the clinical traits.
    16. fpkmSamples = rownames(datExpr0)
    17. traitSamples =rownames(traitData)
    18. #匹配样本名称，性状数据与表达数据保证一致
    19. traitRows = match(fpkmSamples, traitSamples)
    20. datTraits = traitData[traitRows,]
    21. rownames(datTraits) 
    22. collectGarbage()

    增加形状信息后，再次聚类

    1.  
    2. # Re-cluster samples
    3. sampleTree2 = hclust(dist(datExpr0), method = "average")
    4. # Convert traits to a color representation: white means low, red means high, grey means missing entry
    5. traitColors = numbers2colors(datTraits, signed = FALSE)
    6. # Plot the sample dendrogram and the colors underneath.
    7.  
    8. #sizeGrWindow(20,20)
    9. ##pdf(file="2_Sample dendrogram and trait heatmap.pdf",width=12,height=12)
    10. plotDendroAndColors(sampleTree2, traitColors,
    11.                     groupLabels = names(datTraits),
    12.                     main = "Sample dendrogram and trait heatmap")
    13.  
    14. dev.off()

    下方红色，大致分成了两类，效果不错。

    网络构建

    1. #############################network constr########################################
    2.  
    3. # Allow multi-threading within WGCNA. At present this call is necessary.
    4. # Any error here may be ignored but you may want to update WGCNA if you see one.
    5. # Caution: skip this line if you run RStudio or other third-party R environments.
    6. # See note above.
    7. enableWGCNAThreads()
    8.  
    9.  
    10. # Choose a set of soft-thresholding powers
    11. powers = c(1:15)
    12.  
    13. # Call the network topology analysis function
    14. sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
    15.  
    16. # Plot the results:
    17. sizeGrWindow(15, 9)
    18. #pdf(file="3_Scale independence.pdf",width=9,height=5)
    19. #pdf(file="Rplot03.pdf",width=9,height=5)
    20. par(mfrow = c(1,2))
    21. cex1 = 0.9
    22. # Scale-free topology fit index as a function of the soft-thresholding power
    23. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    24.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    25.      main = paste("Scale independence"));
    26. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    27.      labels=powers,cex=cex1,col="red");
    28. # this line corresponds to using an R^2 cut-off of h
    29. abline(h=0.90,col="red")
    30. # Mean connectivity as a function of the soft-thresholding power
    31. plot(sft$fitIndices[,1], sft$fitIndices[,5],
    32.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    33.      main = paste("Mean connectivity"))
    34. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    35. dev.off()

    选取soft-thresholding powers

    测试阈值，注意观察，突破红线的附近时取值，下方代码时候的是自适应的方法选取 soft-thresholding powers

    1. ######chose the softPower
    2. #datExpr0= datExpr0[,-1]
    3. softPower =sft$powerEstimate
    4. adjacency = adjacency(datExpr0, power = softPower)
    5.  
    6. ##### Turn adjacency into topological overlap
    7. TOM = TOMsimilarity(adjacency);
    8. dissTOM = 1-TOM
    9.  
    10. # Call the hierarchical clustering function
    11. geneTree = hclust(as.dist(dissTOM), method = "average");
    12. # Plot the resulting clustering tree (dendrogram)
    13.  
    14. #sizeGrWindow(12,9)
    15. pdf(file="4_Gene clustering on TOM-based dissimilarity.pdf",width=12,height=9)
    16. plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
    17.      labels = FALSE, hang = 0.04)
    18. dev.off()

    基于tom的差异的基因聚类，绘制聚类树

    根据聚类情况，设置颜色

    1.  
    2. # We like large modules, so we set the minimum module size relatively high:
    3. minModuleSize = 30
    4. # Module identification using dynamic tree cut:
    5. dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
    6.                             deepSplit = 2, pamRespectsDendro = FALSE,
    7.                             minClusterSize = minModuleSize);
    8. table(dynamicMods)
    9.  
    10.  
    11.  
    12. # Convert numeric lables into colors
    13. dynamicColors = labels2colors(dynamicMods)
    14. table(dynamicColors)
    15. # Plot the dendrogram and colors underneath
    16. #sizeGrWindow(8,6)
    17. pdf(file="5_Dynamic Tree Cut.pdf",width=8,height=6)
    18. plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
    19.                     dendroLabels = FALSE, hang = 0.03,
    20.                     addGuide = TRUE, guideHang = 0.05,
    21.                     main = "Gene dendrogram and module colors")
    22. dev.off()

    计算eigengenes

    1.  
    2. # Calculate eigengenes
    3. MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
    4. MEs = MEList$eigengenes
    5. # Calculate dissimilarity of module eigengenes
    6. MEDiss = 1-cor(MEs);
    7. # Cluster module eigengenes
    8. METree = hclust(as.dist(MEDiss), method = "average")
    9. # Plot the result
    10. #sizeGrWindow(7, 6)
    11. pdf(file="6_Clustering of module eigengenes.pdf",width=7,height=6)
    12. plot(METree, main = "Clustering of module eigengenes",
    13.      xlab = "", sub = "")
    14. MEDissThres = 0.25######剪切高度可修改
    15. # Plot the cut line into the dendrogram
    16. abline(h=MEDissThres, col = "red")
    17. dev.off()

    模块的自动合并

    1.  
    2. # Call an automatic merging function
    3. merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
    4. # The merged module colors
    5. mergedColors = merge$colors
    6. # Eigengenes of the new merged modules:
    7. mergedMEs = merge$newMEs
    8.  
    9. #sizeGrWindow(12, 9)
    10. pdf(file="7_merged dynamic.pdf", width = 9, height = 6)
    11. plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
    12.                     c("Dynamic Tree Cut", "Merged dynamic"),
    13.                     dendroLabels = FALSE, hang = 0.03,
    14.                     addGuide = TRUE, guideHang = 0.05)
    15. dev.off()
    16.  
    17.  
    18. # Rename to moduleColors
    19. moduleColors = mergedColors
    20. # Construct numerical labels corresponding to the colors
    21. colorOrder = c("grey", standardColors(50))
    22. moduleLabels = match(moduleColors, colorOrder)-1
    23. MEs = mergedMEs
    24.  
    25. # Save module colors and labels for use in subsequent parts
    26. save(MEs, TOM, dissTOM,  moduleLabels, moduleColors, geneTree, sft, file = "networkConstruction-stepByStep.RData")

    模块与临床形状的关系热图 （关键数据）

    1. #############################relate modules to external clinical triats######################################
    2. # Define numbers of genes and samples
    3. nGenes = ncol(datExpr0)
    4. nSamples = nrow(datExpr0)
    5.  
    6. moduleTraitCor = cor(MEs, datTraits, use = "p")
    7. moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
    8.  
    9. #sizeGrWindow(10,6)
    10. pdf(file="8_Module-trait relationships.pdf",width=10,height=6)
    11. # Will display correlations and their p-values
    12. textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
    13.                    signif(moduleTraitPvalue, 1), ")", sep = "")
    14.  
    15. dim(textMatrix) = dim(moduleTraitCor)
    16. par(mar = c(6, 8.5, 3, 3))
    17.  
    18. # Display the correlation values within a heatmap plot #修改性状类型 data.frame
    19. labeledHeatmap(Matrix = moduleTraitCor,
    20.                xLabels = names(data.frame(datTraits)),
    21.                yLabels = names(MEs),
    22.                ySymbols = names(MEs),
    23.                colorLabels = FALSE,
    24.                colors = greenWhiteRed(50),
    25.                textMatrix = textMatrix,
    26.                setStdMargins = FALSE,
    27.                cex.text = 0.5,
    28.                zlim = c(-1,1),
    29.                main = paste("Module-trait relationships"))
    30. dev.off()

    挑选相关性最高的，具有统计学意义的（p<0.05）,red模块最佳！

    红色模块样本表达情况（相关性大）

    产生了很多数据（各个模块的和临床性状）

    1. ######## Define variable weight containing all column of datTraits
    2.  
    3. ###MM and GS
    4.  
    5.  
    6. # names (colors) of the modules
    7. modNames = substring(names(MEs), 3)
    8.  
    9. geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"))
    10. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
    11.  
    12. names(geneModuleMembership) = paste("MM", modNames, sep="")
    13. names(MMPvalue) = paste("p.MM", modNames, sep="")
    14.  
    15. #names of those trait
    16. traitNames=names(data.frame(datTraits))
    17. class(datTraits)
    18.  
    19. geneTraitSignificance = as.data.frame(cor(datExpr0, datTraits, use = "p"))
    20. GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
    21.  
    22. names(geneTraitSignificance) = paste("GS.", traitNames, sep="")
    23. names(GSPvalue) = paste("p.GS.", traitNames, sep="")
    24.  
    25.  
    26. ####plot MM vs GS for each trait vs each module
    27.  
    28.  
    29. ##########example:royalblue and CK
    30. module="red"
    31. column = match(module, modNames)
    32. moduleGenes = moduleColors==module
    33.  
    34. trait="TB"
    35. traitColumn=match(trait,traitNames)
    36.  
    37. sizeGrWindow(7, 7)
    38.  
    39. #par(mfrow = c(1,1))
    40. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    41. abs(geneTraitSignificance[moduleGenes, traitColumn]),
    42. xlab = paste("Module Membership in", module, "module"),
    43. ylab = paste("Gene significance for ",trait),
    44. main = paste("Module membership vs. gene significance\n"),
    45. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    46. ######
    47.  
    48. for (trait in traitNames){
    49.   traitColumn=match(trait,traitNames)
    50.   
    51.   for (module in modNames){
    52.     column = match(module, modNames)
    53.     moduleGenes = moduleColors==module
    54.     
    55.     if (nrow(geneModuleMembership[moduleGenes,]) > 1){####进行这部分计算必须每个模块内基因数量大于2，由于前面设置了最小数量是30，这里可以不做这个判断，但是grey有可能会出现1个gene,它会导致代码运行的时候中断，故设置这一步
    56.       
    57.       #sizeGrWindow(7, 7)
    58.       pdf(file=paste("9_", trait, "_", module,"_Module membership vs gene significance.pdf",sep=""),width=7,height=7)
    59.       par(mfrow = c(1,1))
    60.       verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    61.                          abs(geneTraitSignificance[moduleGenes, traitColumn]),
    62.                          xlab = paste("Module Membership in", module, "module"),
    63.                          ylab = paste("Gene significance for ",trait),
    64.                          main = paste("Module membership vs. gene significance\n"),
    65.                          cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    66.       dev.off()
    67.     }
    68.   }
    69. }
    70.  
    71. #####
    72. names(datExpr0)
    73. probes = names(datExpr0)
    74.  
    75.  
    76. #################export GS and MM############### 
    77.  
    78. geneInfo0 = data.frame(probes= probes,
    79.                        moduleColor = moduleColors)
    80.  
    81. for (Tra in 1:ncol(geneTraitSignificance))
    82. {
    83.   oldNames = names(geneInfo0)
    84.   geneInfo0 = data.frame(geneInfo0, geneTraitSignificance[,Tra],
    85.                          GSPvalue[, Tra])
    86.   names(geneInfo0) = c(oldNames,names(geneTraitSignificance)[Tra],
    87.                        names(GSPvalue)[Tra])
    88. }
    89.  
    90. for (mod in 1:ncol(geneModuleMembership))
    91. {
    92.   oldNames = names(geneInfo0)
    93.   geneInfo0 = data.frame(geneInfo0, geneModuleMembership[,mod],
    94.                          MMPvalue[, mod])
    95.   names(geneInfo0) = c(oldNames,names(geneModuleMembership)[mod],
    96.                        names(MMPvalue)[mod])
    97. }
    98. geneOrder =order(geneInfo0$moduleColor)
    99. geneInfo = geneInfo0[geneOrder, ]
    100.  
    101. write.table(geneInfo, file = "10_GS_and_MM.xls",sep="\t",row.names=F)
    102.  
    1. ####################################################Visualizing the gene network#######################################################
    4. nGenes = ncol(datExpr0)
    5. nSamples = nrow(datExpr0)
    8. # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
    9. plotTOM = dissTOM^7
    10. # Set diagonal to NA for a nicer plot
    11. diag(plotTOM) = NA
    15. # Call the plot function
    17. sizeGrWindow(9,9)  #这个耗电脑内存
    18. pdf(file="12_Network heatmap plot_all gene.pdf",width=9, height=9)
    19. TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
    20. dev.off()
    23. nSelect = 400
    24. # For reproducibility, we set the random seed
    25. set.seed(10)
    26. select = sample(nGenes, size = nSelect)
    27. selectTOM = dissTOM[select, select]
    28. # There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
    29. selectTree = hclust(as.dist(selectTOM), method = "average")
    30. selectColors = moduleColors[select]
    32. # Open a graphical window
    33. #sizeGrWindow(9,9)
    34. # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
    35. # the color palette; setting the diagonal to NA also improves the clarity of the plot
    36. plotDiss = selectTOM^7
    37. diag(plotDiss) = NA
    39. pdf(file="13_Network heatmap plot_selected genes.pdf",width=9, height=9)
    40. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
    41. dev.off()
    45. ####################################################Visualizing the gene network of eigengenes####################################################
    46.  
    47.  
    48. #sizeGrWindow(5,7.5)
    49. pdf(file="14_Eigengene dendrogram and Eigengene adjacency heatmap.pdf", width=5, height=7.5)
    50. par(cex = 0.9)
    51. plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90)
    52. dev.off()
    53.  
    54. #or devide into two parts
    55. # Plot the dendrogram
    56. #sizeGrWindow(6,6);
    57. pdf(file="15_Eigengene dendrogram_2.pdf",width=6, height=6)
    58. par(cex = 1.0)
    59. plotEigengeneNetworks(MEs, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE)
    60. dev.off()
    61.  
    62. pdf(file="15_Eigengene adjacency heatmap_2.pdf",width=6, height=6)
    63. # Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
    64. par(cex = 1.0)
    65. plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)
    66. dev.off()
    67.  
    68.

    后续挖掘核心基因时，需要用到Cytoscape，生成绘图所需要的数据

    1. ###########################Exporting to Cytoscape all one by one ##########################
    2.  
    3.  
    4.  
    5.  
    6. # Select each module
    7. '''
    8. Error in exportNetworkToCytoscape(modTOM, edgeFile = paste("CytoscapeInput-edges-",  : 
    9.   Cannot determine node names: nodeNames is NULL and adjMat has no dimnames.
    10.   datExpr0 格式需要dataframe
    11. '''
    12. modules =module
    13. for (mod in 1:nrow(table(moduleColors)))
    14. {
    15.   
    16.   modules = names(table(moduleColors))[mod]
    17.   # Select module probes
    18.   probes = names(data.frame(datExpr0))  # 
    19.   inModule = (moduleColors == modules)
    20.   modProbes = probes[inModule]
    21.   modGenes = modProbes
    22.   # Select the corresponding Topological Overlap
    23.   modTOM = TOM[inModule, inModule]
    24.   
    25.   dimnames(modTOM) = list(modProbes, modProbes)
    26.   # Export the network into edge and node list files Cytoscape can read
    27.   cyt = exportNetworkToCytoscape(modTOM,
    28.                                  edgeFile = paste("CytoscapeInput-edges-", modules , ".txt", sep=""),
    29.                                  nodeFile = paste("CytoscapeInput-nodes-", modules, ".txt", sep=""),
    30.                                  weighted = TRUE,
    31.                                  threshold = 0.02,
    32.                                  nodeNames = modProbes,
    33.                                  altNodeNames = modGenes,
    34.                                  nodeAttr = moduleColors[inModule])
    35. }

    关系网络的构建完毕，绘图找核心基因，Cytoscape 到底怎么玩？

  • 相关阅读:
    MySQL基本操作
    MySQL案例详解 三：MMM高可用架构及其故障切换
    Jmeter使用
    docker 查看容器挂载情况
    阿里云分布式深度学习训练架构Whale
    [CVE-2021-45105] Apache Log4j2 漏洞复现与原理详细分析
    Benchmarking Detection Transfer Learning with Vision Transformers（2021-11）
    Java基础this关键字02
    如何使用oracle数据库的length()、lengthb()、replace()、regexp_substr()函数
    Ajax（2）
  • 原文地址：https://blog.csdn.net/zzh1464501547/article/details/133872767