生信人

找回密码
立即注册
搜索
热搜: 活动 交友 discuz
生信人»论坛 › 生信基础 › 可视化 › 全代码WGCNA分析流程报告
发新帖
生信喵

0

收听

12

听众

309

主题

全代码WGCNA分析流程报告

发表于 2021-10-20 08:48:04 | 查看: 4803| 回复: 3
本帖最后由 生信喵 于 2021-10-20 08:52 编辑

第一步：安装所需的R包
  1. rm(list = ls()) #### 魔幻操作，一键清空~
  2. getwd()
  3. setwd('/home/xhs/helix/level2/L64_Codes/6_4_WGCNA_Analysis/')
  4. # 设置国内镜像，安装时运行一次即可
  5. options("repos"="https://mirrors.ustc.edu.cn/CRAN/")

  7. if(!"BiocManager"%in%installed.packages())
  8. install.packages("BiocManager",update = F,ask = F)

  10. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")

  12. # 如果使用R4.0版本需要安装Rtools40
  13. # 下载网站https://cran.r-project.org/bin/windows/Rtools/

  15. # 安装GO.db
  16. if(!"GO.db"%in%installed.packages())
  17. BiocManager::install("GO.db")

  19. # 安装flashClust
  20. if(!"flashClust"%in%installed.packages())
  21. BiocManager::install("flashClust")

  23. # 安装WGCNA
  24. if(!"WGCNA"%in%installed.packages())
  25. BiocManager::install("WGCNA")

  27. # 安装org.Hs.eg.db
  28. if(!"org.Hs.eg.db"%in%installed.packages())
  29. BiocManager::install("org.Hs.eg.db")

  31. # 安装clusterProfiler
  32. if(!"clusterProfiler"%in%installed.packages())
  33. BiocManager::install("clusterProfiler")

  35. # 安装最新版rlang
  36. BiocManager::install("rlang")

  38. # 安装富集分析包
  39. source("https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/GeneAnnotation/installAnRichment.R");
  40. # 上面无法加载成功时，用以下代码本地加载
  41. source("installAnRichment.R")
  42. installAnRichment()

  44. install.packages("anRichmentMethods_0.91-94.tar.gz", repos = NULL, type = "source")
  45. install.packages("anRichment_1.10-1.tar.gz", repos = NULL, type = "source")

  47. # 判断当前目录下是否存在data文件夹
  48. # 存在时，忽略
  49. # 不存在，创建data文件夹储存输入文件和结果文件
  50. if(!dir.exists("data")) dir.create("data")
  51. # 测试figures文件夹是否存在，
  52. # 存在时忽略，不存在时创建
  53. # figures文件夹储存所有的结果图片
  54. if(!dir.exists("figures")) dir.create("figures")

复制代码
第二步：整理数据
  1. rm(list = ls())
  2. library(WGCNA)

  4. # 读取基因表达矩阵数据
  5. fpkm = read.table("data/fpkm.txt", header = T, row.names = 1, check.names = F)
  6. head(fpkm)

  8. ### 选取基因方法 ###

  10. ## 第一种，通过标准差选择
  11. ## 计算每个基因的标准差
  12. fpkm_sd = apply(fpkm,1,sd)#1是对每一行，2是对每一列
  13. head(WGCNA_matrix)
  14. ## 使用标准差对基因进行降序排序
  15. fpkm_sd_sorted = order(fpkm_sd, decreasing = T)
  16. ## 选择前5000个标准差较大的基因
  17. fpkm_num = fpkm_sd_sorted[1:5000]
  18. ## 从表达矩阵提取基因
  19. fpkm_filter = fpkm[fpkm_num,]
  20. ## 对表达矩阵进行转置
  21. WGCNA_matrix = t(fpkm_filter)#变成行名是样本，列名是基因
  22. ## 保存过滤后的数据
  23. save(WGCNA_matrix, file = "data/Step01-fpkm_sd_filter.Rda")

  25. ## 第二种，使用绝对中位差选择，推荐使用绝对中位差
  26. WGCNA_matrix = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:5000],])#mad代表绝对中位差
  27. save(WGCNA_matrix, file = "data/Step01-fpkm_mad_filter.Rda")

  29. ## 第三种，使用全基因
  30. WGCNA_matrix = t(fpkm)
  31. save(WGCNA_matrix, file = "data/Step01-fpkm_allgene.Rda")


  34. ### 去除缺失值较多的基因/样本 ###
  35. rm(list = ls())
  36. # 加载表达矩阵
  37. load(file = "data/Step01-fpkm_mad_filter.Rda")
  38. # 加载WGCNA包
  39. library(WGCNA)
  40. # 判断是否缺失较多的样本和基因
  41. datExpr0 = WGCNA_matrix
  42. gsg = goodSamplesGenes(datExpr0, verbose = 3)

  44. # 是否需要过滤，TRUE表示不需要，FALSE表示需要
  45. gsg$allOK
  46. # 当gsg$allOK为TRUE，以下代码不会运行，为FALSE时，运行以下代码过滤样本
  47. if (!gsg$allOK)
  48. {
  49. # 打印移除的基因名和样本名
  50. if (sum(!gsg$goodGenes)>0)
  51. printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
  52. if (sum(!gsg$goodSamples)>0)
  53. printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
  54. # 提取保留的基因和样本
  55. datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
  56. }

  58. ### 通过样本聚类识别离群样本，去除离群样本 ###
  59. sampleTree = hclust(dist(datExpr0), method = "average");#使用hclust函数进行均值聚类
  60. # 绘制样本聚类图确定离群样本
  61. sizeGrWindow(30,9)
  62. pdf(file = "figures/Step01-sampleClustering.pdf", width = 30, height = 9);
  63. par(cex = 0.6);
  64. par(mar = c(0,4,2,0))
  65. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
  66. cex.axis = 1.5, cex.main = 2)
  67. # 根据上图判断，需要截取的高度参数h
  68. abline(h = 120, col = "red")#在120的地方画条线
  69. dev.off()

  71. # 去除离群得聚类样本，cutHeight参数要与上述得h参数值一致
  72. clust = cutreeStatic(sampleTree, cutHeight = 120, minSize = 10)
  73. table(clust)
  74. # clust
  75. # 0 1 0就是要去除的，1就是要保存的
  76. # 15 162

  78. # clust 1聚类中包含着我们想要得样本，将其提取出来
  79. keepSamples = (clust==1)
  80. datExpr = datExpr0[keepSamples, ]

  82. # 记录基因和样本数，方便后续可视化
  83. nGenes = ncol(datExpr)#基因数
  84. nSamples = nrow(datExpr)#样本数
  85. save(datExpr, nGenes, nSamples,file = "data/Step01-WGCNA_input.Rda")
复制代码
第三步：一步法WGCNA
  1. # 清空所有变量
  2. rm(list = ls())
  3. # 加载包
  4. library(WGCNA)
  5. # 允许多线程运行
  6. enableWGCNAThreads()
  7. # 加载表达矩阵
  8. load("data/Step01-WGCNA_input.Rda")

  10. ### 选择软阈值 ###
  11. powers = c(c(1:10), seq(from = 12, to=20, by=2))
  12. # 进行网络拓扑分析
  13. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)#β=power,就是软阈值

  15. # 可视化结果
  16. sizeGrWindow(9, 5)
  17. pdf(file = "figures/Step02-SoftThreshold.pdf", width = 9, height = 5);

  19. par(mfrow = c(1,2))#一个画板上，画两个图，一行两列
  20. cex1 = 0.9;
  21. # 无尺度网络阈值得选择
  22. plot(sft$fitIndices[,1],
  23. -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  24. xlab="Soft Threshold (power)",#x轴
  25. ylab="Scale Free Topology Model Fit,signed R^2",type="n",#y轴
  26. main = paste("Scale independence"));#标题
  27. text(sft$fitIndices[,1],
  28. -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
  29. labels=powers,
  30. cex=cex1,
  31. col="red");

  33. # 用红线标出R^2的参考值
  34. abline(h=0.90,col="red")
  35. # 平均连接度
  36. plot(sft$fitIndices[,1],
  37. sft$fitIndices[,5],
  38. xlab="Soft Threshold (power)",
  39. ylab="Mean Connectivity",
  40. type="n",
  41. main = paste("Mean connectivity"))
  42. text(sft$fitIndices[,1],
  43. sft$fitIndices[,5],
  44. labels=powers,
  45. cex=cex1,
  46. col="red")

  48. dev.off()

  50. # 无尺度网络检验，验证构建的网络是否是无尺度网络
  51. softpower=sft$powerEstimate

  53. ADJ = abs(cor(datExpr,use="p"))^softpower#相关性取绝对值再幂次
  54. k = as.vector(apply(ADJ,2,sum,na.rm=T))#对ADJ的每一列取和，也就是频次

  56. pdf(file = "figures/Step02-scaleFree.pdf",width = 14)
  57. par(mfrow = c(1,2))

  59. hist(k)#直方图
  60. scaleFreePlot(k,main="Check scale free topology")

  62. dev.off()


  65. ### 一步构建网络 ###
  66. net = blockwiseModules(datExpr, #处理好的表达矩阵
  67. power = sft$powerEstimate,#选择的软阈值
  68. TOMType = "unsigned", #拓扑矩阵类型，none表示邻接矩阵聚类，unsigned最常用，构建无方向
  69. minModuleSize = 30,#网络模块包含的最少基因数
  70. reassignThreshold = 0, #模块间基因重分类的阈值
  71. mergeCutHeight = 0.25,#合并相异度低于0.25的模块
  72. numericLabels = TRUE, #true，返回模块的数字标签 false返回模块的颜色标签
  73. pamRespectsDendro = FALSE,#调用动态剪切树算法识别网络模块后，进行第二次的模块比较，合并相关性高的模块
  74. saveTOMs = TRUE,#保存拓扑矩阵
  75. saveTOMFileBase = "data/Step02-fpkmTOM",
  76. verbose = 3)#0，不反回任何信息，＞0返回计算过程
  77. # 保存网络构建结果
  78. save(net, file = "data/Step02-One_step_net.Rda")

  80. # 加载网络构建结果
  81. load(file = "data/Step02-One_step_net.Rda")
  82. # 打开绘图窗口
  83. sizeGrWindow(12, 9)
  84. pdf(file = "figures/Step02-moduleCluster.pdf", width = 12, height = 9);
  85. # 将标签转化为颜色
  86. mergedColors = labels2colors(net$colors)
  87. # 绘制聚类和网络模块对应图
  88. plotDendroAndColors(dendro = net$dendrograms[[1]], #hclust函数生成的聚类结果
  89. colors = mergedColors[net$blockGenes[[1]]],#基因对应的模块颜色
  90. groupLabels = "Module colors",#分组标签
  91. dendroLabels = FALSE, #false,不显示聚类图的每个分支名称
  92. hang = 0.03,#调整聚类图分支所占的高度
  93. addGuide = TRUE, #为聚类图添加辅助线
  94. guideHang = 0.05,#辅助线所在高度
  95. main = "Gene dendrogram and module colors")
  96. dev.off()

  98. # 加载TOM矩阵
  99. load("data/Step02-fpkmTOM-block.1.RData")
  100. # 网络特征向量
  101. MEs = moduleEigengenes(datExpr, mergedColors)$eigengenes
  102. # 对特征向量排序
  103. MEs = orderMEs(MEs)
  104. # 可视化模块间的相关性
  105. sizeGrWindow(5,7.5);
  106. pdf(file = "figures/Step02-moduleCor.pdf", width = 5, height = 7.5);

  108. par(cex = 0.9)
  109. plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2),
  110. marHeatmap = c(3,4,1,2), cex.lab = 0.8,
  111. xLabelsAngle = 90)
  112. dev.off()
复制代码
发表于 2021-10-20 08:49:16
  1. ## TOMplot

  3. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate); #1-相关性=相异性

  5. nSelect = 400
  6. # 随机选取400个基因进行可视化，设置seed值，保证结果的可重复性
  7. set.seed(10);#设置随机种子数
  8. select = sample(nGenes, size = nSelect);#从5000个基因选择400个
  9. selectTOM = dissTOM[select, select];#选择这400*400的矩阵
  10. # 对选取的基因进行重新聚类
  11. selectTree = hclust(as.dist(selectTOM), method = "average")#用hclust重新聚类
  12. selectColors = mergedColors[select];#提取相应的颜色模块
  13. # 打开一个绘图窗口
  14. sizeGrWindow(9,9)
  15. pdf(file = "figures/Step02-TOMplot.pdf", width = 9, height = 9);
  16. # 美化图形的设置
  17. plotDiss = selectTOM^7;
  18. diag(plotDiss) = NA;
  19. # 绘制TOM图
  20. TOMplot(plotDiss, #拓扑矩阵，该矩阵记录了每个节点之间的相关性
  21. selectTree, #基因的聚类结果
  22. selectColors, #基因对应的模块颜色
  23. main = "Network heatmap plot, selected genes")
  24. dev.off()
复制代码
第四步：分步法WGCNA
  1. # 清空环境变量
  2. rm(list = ls())
  3. #加载R包
  4. library(WGCNA)
  5. # 允许多线程运行
  6. enableWGCNAThreads()
  7. # 加载表达矩阵
  8. load("data/Step01-WGCNA_input.Rda")

  10. ## 选择软阈值
  11. sizeGrWindow(9,5);
  12. par(mfrow = c(1,2));
  13. powers = c(c(1:10), seq(from = 12, to=20, by=2))
  14. RpowerTable=pickSoftThreshold(datExpr, powerVector=powers)[[2]]
  15. cex1=0.7
  16. pdf(file="figures/Step03-softThresholding.pdf",width = 14)
  17. par(mfrow=c(1,2))
  18. plot(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",main = paste("Scale independence"))
  19. text(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2], labels=powers,cex=cex1,col="red")

  21. abline(h=0.9,col="red")
  22. plot(RpowerTable[,1], RpowerTable[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))
  23. text(RpowerTable[,1], RpowerTable[,5], labels=powers, cex=cex1,col="red")
  24. dev.off()

  26. ## 无尺度网络验证
  27. softpower=7

  29. ADJ = abs(cor(datExpr,use="p"))^softpower
  30. k = as.vector(apply(ADJ,2,sum,na.rm=T))

  32. pdf(file = "figures/Step03-scaleFree.pdf",width = 14)
  33. par(mfrow = c(1,2))
  34. hist(k)
  35. scaleFreePlot(k,main="Check scale free topology")
  36. dev.off()

  38. ## 计算邻接矩阵
  39. adjacency = adjacency(datExpr,power=softpower)
  40. ## 计算TOM拓扑矩阵
  41. TOM = TOMsimilarity(adjacency)
  42. ## 计算相异度
  43. dissTOM = 1- TOM

  45. #模块初步聚类分析
  46. library(flashClust)
  47. geneTree = flashClust(as.dist(dissTOM),method="average")
  48. #绘制层次聚类树
  49. pdf(file = "figures/Step03-GeneClusterTOM-based.pdf")
  50. plot(geneTree,xlab="",sub="",main="Gene clustering on TOM-based",labels=FALSE,hang=0.04)
  51. dev.off()

  53. #构建初步基因模块
  54. #设定基因模块中至少30个基因
  55. minModuleSize=30
  56. # 动态剪切树识别网络模块
  57. dynamicMods = cutreeDynamic(dendro = geneTree,#hclust函数的聚类结果
  58. distM = dissTOM,#
  59. deepSplit = 2,
  60. pamRespectsDendro = FALSE,
  61. minClusterSize = minModuleSize)#设定基因模块中至少30个基因
  62. # 将标签转换为颜色
  63. dynamicColors = labels2colors(dynamicMods)
  64. table(dynamicColors)#看聚类到哪些模块，哪些颜色

  66. pdf(file="figures/Step03-DynamicTreeCut.pdf",width=9,height=5)
  67. plotDendroAndColors(dendro = geneTree,
  68. colors = dynamicColors,
  69. groupLabels = "Dynamic Tree Cut",
  70. dendroLabels = FALSE, hang = 0.03,
  71. addGuide = TRUE,
  72. main = "Gene dendrogram and module colors")
  73. dev.off()


  76. #前面用动态剪切树聚类了一些模块，现在要对这步结果进一步合并，合并相似度大于0.75的模块，降低网络的复杂度
  77. #计算基因模块特征向量
  78. MEList = moduleEigengenes(datExpr,colors=dynamicColors)#计算特征向量
  79. MEs = MEList$eigengenes;#提取特征向量
  80. MEDiss1 = 1-cor(MEs);#计算相异度
  81. METree1 = flashClust(as.dist(MEDiss1),method="average")#对相异度进行flashClust聚类
  82. #设置特征向量相关系数大于0.75
  83. MEDissThres = 0.25;#相异度在0.25以下，也就是相似度大于0.75，对这些模块合并
  84. #合并模块
  85. merge = mergeCloseModules(datExpr, #合并相似度大于0.75的模块
  86. dynamicColors,
  87. cutHeight = MEDissThres,
  88. verbose=3)
  89. mergedColors = merge$colors

  91. table(dynamicColors)#动态剪切树的模块颜色
  92. table(mergedColors)#合并后的模块颜色，可以看到从18个模块变成了14个模块

  94. mergedMEs = merge$newMEs#合并后的14个模块

  96. #重新命名合并后的模块
  97. moduleColors = mergedColors;
  98. colorOrder = c("grey",standardColors(50));
  99. moduleLabels = match(moduleColors,colorOrder)-1;
  100. MEs = mergedMEs;
  101. MEDiss2 = 1-cor(MEs);#计算相异度
  102. METree2 = flashClust(as.dist(MEDiss2),method="average");#对合并后的模块进行聚类

  104. #绘制聚类结果图
  105. pdf(file="figures/Step03-MECombined.pdf",width=12,height=5)
  106. par(mfrow=c(1,2))
  107. plot(METree1,xlab="",sub="",main="Clustering of ME before combined")# METree1是动态剪切树形成的模块
  108. abline(h=MEDissThres,col="red")#相异度为0.25
  109. plot(METree2,xlab="",sub="",main="Clustering of ME after combined")# METree2是合并后的模块
  110. dev.off()

  112. pdf(file="figures/Step03-MergedDynamics.pdf",width=10,height=4)
  113. plotDendroAndColors(dendro = geneTree,#剪切树
  114. colors = cbind(dynamicColors,mergedColors),#将两种方法形成的模块颜色合并在一起
  115. groupLabels = c("Dynamic Tree Cut","Merged Dynamics"),
  116. dendroLabels = FALSE,
  117. hang = 0.03,
  118. addGuide=TRUE,
  119. guideHang=0.05,
  120. main="Gene Dendrogram and module colors")
  121. dev.off()

  123. # 模块中基因数
  124. write.table(table(moduleColors),"data/Step03-MEgeneCount.txt",quote = F,row.names = F)

  126. # 保存构建的网络信息
  127. moduleColors=mergedColors
  128. colorOrder=c("grey", standardColors(50))
  129. moduleLabels=match(moduleColors, colorOrder)-1
  130. MEs=mergedMEs
  131. save(MEs, moduleLabels, moduleColors, geneTree, file="data/Step03-Step_by_step_buildnetwork.rda")




  136. ## 绘制样本间的相关性
  137. load("data/Step01-WGCNA_input.Rda")
  138. load(file="data/Step03-Step_by_step_buildnetwork.rda")


  141. MEs = orderMEs(MEs)

  143. sizeGrWindow(5,7.5);
  144. pdf(file = "figures/Step03-moduleCor.pdf", width = 5, height = 7.5);

  146. par(cex = 0.9)
  147. plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
  148. = 90)
  149. dev.off()

  151. ## TOMplot

  153. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = softpower);

  155. nSelect = 400
  156. # 随机选取400个基因进行可视化，s设置seed值，保证结果的可重复性
  157. set.seed(10);
  158. select = sample(nGenes, size = nSelect);
  159. selectTOM = dissTOM[select, select];
  160. # 对选取的基因进行重新聚类
  161. selectTree = hclust(as.dist(selectTOM), method = "average")
  162. selectColors = mergedColors[select];
  163. # 打开一个绘图窗口
  164. sizeGrWindow(9,9)
  165. pdf(file = "figures/Step03-TOMplot.pdf", width = 9, height = 9);
  166. # 美化图形的设置
  167. plotDiss = selectTOM^7;
  168. diag(plotDiss) = NA;
  169. TOMplot(plotDiss, selectTree, selectColors,
  170. main = "Network heatmap plot, selected genes")
  171. dev.off()
复制代码


回复 显示全部楼层 道具 举报

发表于 2021-10-20 08:50:50
第五步：模块功能分析GO和KEGG
  1. # 清空环境变量
  2. rm(list = ls())
  3. getwd()
  4. setwd('/home/xhs/helix/level2/L64_Codes/6_4_WGCNA_Analysis/')
  5. # 加载包
  6. library(WGCNA)
  7. # 加载表达矩阵
  8. load("data/Step01-WGCNA_input.Rda")

  10. # 读入临床信息
  11. clinical = read.table("data/clinical.txt",stringsAsFactors = TRUE, header = T,row.names = 1,na.strings = "",sep = "\t")
  12. # 查看临床信息
  13. head(clinical)
  14. # 对表达矩阵进行预处理
  15. datTraits = as.data.frame(do.call(cbind,lapply(clinical, as.numeric)))
  16. rownames(datTraits) = rownames(clinical)

  18. # 对样本进行聚类
  19. sampleTree2 = hclust(dist(datExpr), method = "average")

  21. # 将临床信息转换为颜色，白色表示低，红色表示高，灰色表示缺失
  22. traitColors = numbers2colors(datTraits, signed = FALSE)

  24. pdf(file = "figures/Step04-Sample_dendrogram_and_trait_heatmap.pdf", width = 24);
  25. # 样本聚类图与样本性状热图
  26. plotDendroAndColors(sampleTree2,
  27. traitColors,
  28. groupLabels = names(datTraits),
  29. main = "Sample dendrogram and trait heatmap")
  30. dev.off()

  32. #### 网络的分析
  33. ###### 基因模块与临床信息的关系
  34. # 加载构建的网络
  35. load(file = "data/Step03-Step_by_step_buildnetwork.rda")
  36. # 对模块特征矩阵进行排序
  37. MEs=orderMEs(MEs)
  38. #计算模型特征矩阵和样本信息矩阵的相关度。
  39. moduleTraitCor=cor(MEs, datTraits, use="p")
  40. write.table(file="data/Step04-modPhysiological.cor.xls",moduleTraitCor,sep="\t",quote=F)
  41. moduleTraitPvalue=corPvalueStudent(moduleTraitCor, nSamples)
  42. write.table(file="data/Step04-modPhysiological.p.xls",moduleTraitPvalue,sep="\t",quote=F)

  44. #使用labeledHeatmap()将上述相关矩阵和p值可视化。
  45. pdf(file="figures/Step04-Module_trait_relationships.pdf",width=9,height=7)
  46. textMatrix=paste(signif(moduleTraitCor,2),"\n(",signif(moduleTraitPvalue,1),")",sep="")
  47. dim(textMatrix)=dim(moduleTraitCor)
  48. # 基因模块与临床信息相关性图
  49. labeledHeatmap(Matrix=moduleTraitCor,#模块和表型的相关性矩阵，这个参数最重要，其他可以不变
  50. xLabels=colnames(datTraits),
  51. yLabels=names(MEs),
  52. ySymbols=names(MEs),
  53. colorLabels=FALSE,
  54. colors=blueWhiteRed(50),
  55. textMatrix=textMatrix,
  56. setStdMargins=FALSE,
  57. cex.text=0.5,
  58. cex.lab=0.5,
  59. zlim=c(-1,1),
  60. main=paste("Module-trait relationships"))
  61. dev.off()


  64. ## 单一模块与某一表型相关性
  65. M_stage = as.data.frame(datTraits$M_stage)
  66. # 分析自己感兴趣的临床信息，此处以M_stage为示例
  67. names(M_stage) = "M_stage"
  68. # 模块对应的颜色
  69. modNames = substring(names(MEs), 3)

  71. # 计算基因模块特征
  72. geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
  73. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

  75. # 对结果进行命名
  76. names(geneModuleMembership) = paste("MM", modNames, sep="");
  77. names(MMPvalue) = paste("p.MM", modNames, sep="");

  79. # 计算M分期基因特征显著性
  80. geneTraitSignificance = as.data.frame(cor(datExpr, M_stage, use = "p"));
  81. GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));

  83. # 对结果进行命名
  84. names(geneTraitSignificance) = paste("GS.", names(M_stage), sep="")
  85. names(GSPvalue) = paste("p.GS.", names(M_stage), sep="")

  87. # 设置需要分析的模块名称，此处为brown模块
  88. module = "brown"
  89. # 提取brown模块数据
  90. column = match(module, modNames);
  91. moduleGenes = moduleColors==module;

  93. # 可视化brown模块与M分期的相关性分析结果
  94. sizeGrWindow(7, 7);
  95. pdf(file="figures/Step04-Module_membership_vs_gene_significance.pdf")
  96. par(mfrow = c(1,1));
  97. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
  98. abs(geneTraitSignificance[moduleGenes, 1]),
  99. xlab = paste("Module Membership in", module, "module"),
  100. ylab = "Gene significance for M Stage",
  101. main = paste("Module membership vs. gene significance\n"),
  102. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  103. dev.off()

  105. ## 单一特征与所有模块关联分析
  106. GS = as.numeric(cor(datTraits$M_stage,datExpr, use="p"))
  107. GeneSignificance = abs(GS);
  108. pdf(file="figures/Step04-Gene_significance_for_M_stage_across_module.pdf",width=9,height=5)
  109. plotModuleSignificance(GeneSignificance, moduleColors, ylim=c(0,0.2), main="Gene significance for M stage across module");
  110. dev.off()

  112. ## 模块中的hub基因
  113. #### 为每一个模块寻找hub基因
  114. HubGenes <- chooseTopHubInEachModule(datExpr,#WGCNA分析输入的表达矩阵
  115. moduleColors)#模块颜色信息
  116. # 保存hub基因结果
  117. write.table (HubGenes,file = "data/Step04-HubGenes_of_each_module.xls",quote=F,sep='\t',col.names = F)

  119. #### 与某种特征相关的hub基因
  120. NS = networkScreening(datTraits$M_stage,#M分期
  121. MEs,#
  122. datExpr)#WGCNA分析输入的表达矩阵
  123. # 将结果写入到文件
  124. write.table(NS,file="data/Step04-Genes_for_M_stage.xls",quote=F,sep='\t')

  126. ## 模块GO/KEGG分析
  127. # 加载R包
  128. library(anRichment)
  129. library(clusterProfiler)
  130. ##### GO分析
  131. # 构建GO背景基因集
  132. GOcollection = buildGOcollection(organism = "human")
  133. geneNames = colnames(datExpr)
  134. # 将基因SYMBOL转换为ENTREZID基因名
  135. geneID = bitr(geneNames,fromType = "SYMBOL", toType = "ENTREZID",
  136. OrgDb = "org.Hs.eg.db", drop = FALSE)
  137. # 将基因名对应结果写入文件中
  138. write.table(geneID, file = "data/Step04-geneID_map.xls", sep = "\t", quote = TRUE, row.names = FALSE)
  139. # 进行GO富集分析
  140. GOenr = enrichmentAnalysis(classLabels = moduleColors,#基因所在的模块信息
  141. identifiers = geneID$ENTREZID,
  142. refCollection = GOcollection,
  143. useBackground = "given",
  144. threshold = 1e-4,
  145. thresholdType = "Bonferroni",
  146. getOverlapEntrez = TRUE,
  147. getOverlapSymbols = TRUE,
  148. ignoreLabels = "grey");
  149. # 提取结果，并写入结果到文件
  150. tab = GOenr$enrichmentTable
  151. names(tab)
  152. write.table(tab, file = "data/Step04-GOEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)

  154. # 提取主要结果，并写入文件
  155. keepCols = c(1, 3, 4, 6, 7, 8, 13)
  156. screenTab = tab[, keepCols]
  157. # 小数位为2位
  158. numCols = c(4, 5, 6)
  159. screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)

  161. # 给结果命名
  162. colnames(screenTab) = c("module", "GOID", "term name", "p-val", "Bonf", "FDR", "size")
  163. rownames(screenTab) = NULL

  165. # 查看结果
  166. head(screenTab)
  167. # 写入文件中
  168. write.table(screenTab, file = "data/Step04-GOEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)

  170. ##### KEGG富集分析
  171. # AnRichment没有直接提供KEGG数据的背景集
  172. # 这里使用MSigDBCollection构建C2通路数据集
  173. KEGGCollection = MSigDBCollection("data/msigdb_v7.1.xml", MSDBVersion = "7.1",
  174. organism = "human",
  175. excludeCategories = c("h","C1","C3","C4","C5","C6","C7"))
  176. # KEGG分析
  177. KEGGenr = enrichmentAnalysis(classLabels = moduleColors,
  178. identifiers = geneID$ENTREZID,
  179. refCollection = KEGGCollection,
  180. useBackground = "given",
  181. threshold = 1e-4,
  182. thresholdType = "Bonferroni",
  183. getOverlapEntrez = TRUE,
  184. getOverlapSymbols = TRUE,
  185. ignoreLabels = "grey")
  186. # 提取KEGG结果，并写入文件
  187. tab = KEGGenr$enrichmentTable
  188. names(tab)
  189. write.table(tab, file = "data/Step04-KEGGEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)

  191. # 提取主要结果并写入文件
  192. keepCols = c(1, 3, 4, 6, 7, 8, 13)
  193. screenTab = tab[, keepCols]
  194. # 取两位有效数字
  195. numCols = c(4, 5, 6)
  196. screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)

  198. # 对结果表格进行重命名
  199. colnames(screenTab) = c("module", "ID", "term name", "p-val", "Bonf", "FDR", "size")
  200. rownames(screenTab) = NULL
  201. # 查看结果
  202. head(screenTab)
  203. # 写入文件中
  204. write.table(screenTab, file = "data/Step04-KEGGEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)

  206. ### 输出cytoscape可视化
  207. # 重新计算TOM，power值设置为前面选择好的
  208. TOM = TOMsimilarityFromExpr(datExpr, power = 7)
  209. # 输出全部网络模块
  210. cyt = exportNetworkToCytoscape(TOM,
  211. edgeFile = "data/Step04-CytoscapeInput-edges-all.txt",#基因间的共表达关系
  212. nodeFile = "data/Step04-CytoscapeInput-nodes-all.txt",#
  213. weighted = TRUE,
  214. threshold = 0.1,
  215. nodeNames = geneID$SYMBOL,
  216. altNodeNames = geneID$ENTREZID,
  217. nodeAttr = moduleColors)
  218. # 输出感兴趣网络模块
  219. modules = c("brown", "red")
  220. # 选择上面模块中包含的基因
  221. inModule = is.finite(match(moduleColors, modules))
  222. modGenes = geneID[inModule,]
  223. # 选择指定模块的TOM矩阵
  224. modTOM = TOM[inModule, inModule]

  226. # 输出为Cytoscape软件可识别格式
  227. cyt = exportNetworkToCytoscape(modTOM,
  228. edgeFile = paste("data/Step04-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
  229. nodeFile = paste("data/Step04-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
  230. weighted = TRUE,
  231. threshold = 0.2,
  232. nodeNames = modGenes$SYMBOL,
  233. altNodeNames = modGenes$ENTREZID,
  234. nodeAttr = moduleColors[inModule])
复制代码
第六步：Cytoscape可视化

回复 显示全部楼层 道具 举报

发表于 2021-10-20 08:51:40
本帖最后由 生信喵 于 2021-11-1 11:57 编辑

代码中需要用到的输入数据：临床信息和转录组数据。获取示例数据请在公众号"生信喵实验柴"后台回复“20211020”。

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有账号？立即注册

回复 显示全部楼层 道具 举报

返回列表
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

使用Markdown编辑器编辑 使用富文本编辑器编辑

QQ|Archiver|手机版|小黑屋|生信人 ( 萌ICP备20244422号 )

GMT+8, 2024-10-30 21:20 , Processed in 0.075320 second(s), 30 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表