全代码WGCNA分析流程报告

生信喵

第一步：安装所需的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/")
   
6. 
   
7. if(!"BiocManager"%in%installed.packages())
   
8. install.packages("BiocManager",update = F,ask = F)
   
9. 
   
10. options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
    
11. 
    
12. # 如果使用R4.0版本需要安装Rtools40
    
13. # 下载网站https://cran.r-project.org/bin/windows/Rtools/
    
14. 
    
15. # 安装GO.db
    
16. if(!"GO.db"%in%installed.packages())
    
17. BiocManager::install("GO.db")
    
18. 
    
19. # 安装flashClust
    
20. if(!"flashClust"%in%installed.packages())
    
21. BiocManager::install("flashClust")
    
22. 
    
23. # 安装WGCNA
    
24. if(!"WGCNA"%in%installed.packages())
    
25. BiocManager::install("WGCNA")
    
26. 
    
27. # 安装org.Hs.eg.db
    
28. if(!"org.Hs.eg.db"%in%installed.packages())
    
29. BiocManager::install("org.Hs.eg.db")
    
30. 
    
31. # 安装clusterProfiler
    
32. if(!"clusterProfiler"%in%installed.packages())
    
33. BiocManager::install("clusterProfiler")
    
34. 
    
35. # 安装最新版rlang
    
36. BiocManager::install("rlang")
    
37. 
    
38. # 安装富集分析包
    
39. source("https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/GeneAnnotation/installAnRichment.R");
    
40. # 上面无法加载成功时，用以下代码本地加载
    
41. source("installAnRichment.R")
    
42. installAnRichment()
    
43. 
    
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")
    
46. 
    
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")
    
55. 
    

复制代码

第二步：整理数据

1. rm(list = ls())
   
2. library(WGCNA)
   
3. 
   
4. # 读取基因表达矩阵数据
   
5. fpkm = read.table("data/fpkm.txt", header = T, row.names = 1, check.names = F)
   
6. head(fpkm)
   
7. 
   
8. ### 选取基因方法 ###
   
9. 
   
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")
    
24. 
    
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")
    
28. 
    
29. ## 第三种，使用全基因
    
30. WGCNA_matrix = t(fpkm)
    
31. save(WGCNA_matrix, file = "data/Step01-fpkm_allgene.Rda")
    
32. 
    
33. 
    
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)
    
43. 
    
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. }
    
57. 
    
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()
    
70. 
    
71. # 去除离群得聚类样本，cutHeight参数要与上述得h参数值一致
    
72. clust = cutreeStatic(sampleTree, cutHeight = 120, minSize = 10)
    
73. table(clust)
    
74. # clust
    
75. # 0 1 0就是要去除的，1就是要保存的
    
76. # 15 162
    
77. 
    
78. # clust 1聚类中包含着我们想要得样本，将其提取出来
    
79. keepSamples = (clust==1)
    
80. datExpr = datExpr0[keepSamples, ]
    
81. 
    
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")
   
9. 
   
10. ### 选择软阈值 ###
    
11. powers = c(c(1:10), seq(from = 12, to=20, by=2))
    
12. # 进行网络拓扑分析
    
13. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)#β=power,就是软阈值
    
14. 
    
15. # 可视化结果
    
16. sizeGrWindow(9, 5)
    
17. pdf(file = "figures/Step02-SoftThreshold.pdf", width = 9, height = 5);
    
18. 
    
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");
    
32. 
    
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")
    
47. 
    
48. dev.off()
    
49. 
    
50. # 无尺度网络检验，验证构建的网络是否是无尺度网络
    
51. softpower=sft$powerEstimate
    
52. 
    
53. ADJ = abs(cor(datExpr,use="p"))^softpower#相关性取绝对值再幂次
    
54. k = as.vector(apply(ADJ,2,sum,na.rm=T))#对ADJ的每一列取和，也就是频次
    
55. 
    
56. pdf(file = "figures/Step02-scaleFree.pdf",width = 14)
    
57. par(mfrow = c(1,2))
    
58. 
    
59. hist(k)#直方图
    
60. scaleFreePlot(k,main="Check scale free topology")
    
61. 
    
62. dev.off()
    
63. 
    
64. 
    
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")
    
79. 
    
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()
    
97. 
    
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);
     
107. 
     
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()
     

复制代码

生信喵

1. ## TOMplot
   
2. 
   
3. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate); #1-相关性=相异性
   
4. 
   
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")
   
9. 
   
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")
    
20. 
    
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()
    
25. 
    
26. ## 无尺度网络验证
    
27. softpower=7
    
28. 
    
29. ADJ = abs(cor(datExpr,use="p"))^softpower
    
30. k = as.vector(apply(ADJ,2,sum,na.rm=T))
    
31. 
    
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()
    
37. 
    
38. ## 计算邻接矩阵
    
39. adjacency = adjacency(datExpr,power=softpower)
    
40. ## 计算TOM拓扑矩阵
    
41. TOM = TOMsimilarity(adjacency)
    
42. ## 计算相异度
    
43. dissTOM = 1- TOM
    
44. 
    
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()
    
52. 
    
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)#看聚类到哪些模块，哪些颜色
    
65. 
    
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()
    
74. 
    
75. 
    
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
    
90. 
    
91. table(dynamicColors)#动态剪切树的模块颜色
    
92. table(mergedColors)#合并后的模块颜色，可以看到从18个模块变成了14个模块
    
93. 
    
94. mergedMEs = merge$newMEs#合并后的14个模块
    
95. 
    
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");#对合并后的模块进行聚类
     
103. 
     
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()
     
111. 
     
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()
     
122. 
     
123. # 模块中基因数
     
124. write.table(table(moduleColors),"data/Step03-MEgeneCount.txt",quote = F,row.names = F)
     
125. 
     
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")
     
132. 
     
133. 
     
134. 
     
135. 
     
136. ## 绘制样本间的相关性
     
137. load("data/Step01-WGCNA_input.Rda")
     
138. load(file="data/Step03-Step_by_step_buildnetwork.rda")
     
139. 
     
140. 
     
141. MEs = orderMEs(MEs)
     
142. 
     
143. sizeGrWindow(5,7.5);
     
144. pdf(file = "figures/Step03-moduleCor.pdf", width = 5, height = 7.5);
     
145. 
     
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()
     
150. 
     
151. ## TOMplot
     
152. 
     
153. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = softpower);
     
154. 
     
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()

复制代码

生信喵

第五步：模块功能分析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")
   
9. 
   
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)
    
17. 
    
18. # 对样本进行聚类
    
19. sampleTree2 = hclust(dist(datExpr), method = "average")
    
20. 
    
21. # 将临床信息转换为颜色，白色表示低，红色表示高，灰色表示缺失
    
22. traitColors = numbers2colors(datTraits, signed = FALSE)
    
23. 
    
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()
    
31. 
    
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)
    
43. 
    
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()
    
62. 
    
63. 
    
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)
    
70. 
    
71. # 计算基因模块特征
    
72. geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
    
73. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
    
74. 
    
75. # 对结果进行命名
    
76. names(geneModuleMembership) = paste("MM", modNames, sep="");
    
77. names(MMPvalue) = paste("p.MM", modNames, sep="");
    
78. 
    
79. # 计算M分期基因特征显著性
    
80. geneTraitSignificance = as.data.frame(cor(datExpr, M_stage, use = "p"));
    
81. GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
    
82. 
    
83. # 对结果进行命名
    
84. names(geneTraitSignificance) = paste("GS.", names(M_stage), sep="")
    
85. names(GSPvalue) = paste("p.GS.", names(M_stage), sep="")
    
86. 
    
87. # 设置需要分析的模块名称，此处为brown模块
    
88. module = "brown"
    
89. # 提取brown模块数据
    
90. column = match(module, modNames);
    
91. moduleGenes = moduleColors==module;
    
92. 
    
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()
     
104. 
     
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()
     
111. 
     
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)
     
118. 
     
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')
     
125. 
     
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)
     
153. 
     
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)
     
160. 
     
161. # 给结果命名
     
162. colnames(screenTab) = c("module", "GOID", "term name", "p-val", "Bonf", "FDR", "size")
     
163. rownames(screenTab) = NULL
     
164. 
     
165. # 查看结果
     
166. head(screenTab)
     
167. # 写入文件中
     
168. write.table(screenTab, file = "data/Step04-GOEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)
     
169. 
     
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)
     
190. 
     
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)
     
197. 
     
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)
     
205. 
     
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]
     
225. 
     
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可视化

生信喵

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