RNA-seq入门实战（六）：GO、KEGG富集分析与enrichplot超全可视化攻略
连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中几个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！
他前面的分享是：
下面是他对我们b站转录组视频课程的详细笔记
承接上节：RNA-seq入门实战（四）：差异分析前的准备——数据检查，以及 RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较
本节概览：1.获取DEG结果的上下调差异基因2.bitr()函数转化基因名为entrez ID3.利用clusterProfiler进行KEGG与GO富集4.用enrichplot进行富集结果可视化：pathview goplot barplot dotplot cnetplot emapplot treeplot heatplot upsetplot
1. 获取DEG结果的上下调差异基因
载入上节RNA-seq入门的简单实战（五）中保存的三种差异分析结果数据，这里示范选取DESeq2的结果数据，进行筛选条件设置后获取上下调基因名
rm(list=ls())options(stringsAsFactors = F) library(clusterProfiler)library(enrichplot)library(tidyverse)# library(org.Hs.eg.db)library(org.Mm.eg.db)library(DOSE)library(pathview) #BiocManager::install("pathview",ask = F,update = F)
##数据载入和目录设置setwd("C:/Users/Lenovo/Desktop/test")source('FUNCTIONS.R')load(file = '1.counts.Rdata')load(list.files(path = "./3.DEG",pattern = 'DEG_results.Rdata',full.names = T))dir.create("4.Enrichment_KEGG_GO")setwd("4.Enrichment_KEGG_GO")
##### 筛选条件设置 #######log2FC_cutoff = log2(10)pvalue_cutoff = 0.001padj_cutoff = 0.001
####获取DEG结果的上下调基因 #####DEG_DESeq2[,c(2,5,6)] DEG_limma_voom[,c(1,4,5)] DEG_edgeR[,c(1,4,5)] need_DEG <- DEG_DESeq2[,c(2,5,6)] head(need_DEG)colnames(need_DEG) <- c('log2FoldChange','pvalue','padj')
gene_up=rownames(need_DEG[with(need_DEG,log2FoldChange>log2FC_cutoff & pvalue<pvalue_cutoff & padj<padj_cutoff),])gene_down=rownames(need_DEG[with(need_DEG,log2FoldChange < -log2FC_cutoff & pvalue<pvalue_cutoff & padj<padj_cutoff),])2. 转化基因名为entrez ID
在进行差异基因富集分析前，需要先将基因名为entrez ID。
转化ID前要载入org.Hs.eg.db\org.Mm.eg.db，其包含着各大主流数据库的数据，如entrez ID和ensembl等等，使用keytypes(org.Mm.eg.db) 可查看所有支持及可转化类型。
一般利用clusterProfiler包的bitr()函数或者mapIds()函数进行转换，用法如下：
#### 转化基因名为entrez ID #####org.Hs.eg.db\org.Mm.eg.db包含着各大主流数据库的数据，如entrez ID和ensembl,#keytypes(org.Hs.eg.db) #查看所有支持及可转化类型 常用 "ENTREZID","ENSEMBL","SYMBOL"gene_up_entrez <- as.character(na.omit(bitr(gene_up, #数据集 fromType="SYMBOL", #输入格式 toType="ENTREZID", # 转为ENTERZID格式 OrgDb="org.Mm.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"gene_down_entrez <- as.character(na.omit(bitr(gene_down, #数据集 fromType="SYMBOL", #输入格式 toType="ENTREZID", # 转为ENTERZID格式 OrgDb="org.Mm.eg.db")[,2])) #"org.Hs.eg.db" "org.Mm.eg.db"gene_diff_entrez <- unique(c(gene_up_entrez ,gene_down_entrez ))##或者使用mapIds# gene_up_ENTREZID <- as.character(na.omit(mapIds(x = org.Mm.eg.db,# keys = gene_up,# keytype = "SYMBOL",# column = "ENTREZID")))3. 利用clusterProfiler进行KEGG与GO富集
有了上一步所得上下调差异基因的entrez ID，我们就可以利用clusterProfiler包的enrichKEGG()和enrichGO()函数进行富集分析了。在这里仅示范对上调基因进行富集，实际应用中可以将上下调和合并的基因都分别进行富集查看结果。需要注意以下事项：
函数的参数pvalueCutoff 默认为 0.05，qvalueCutoff 默认为 0.2，可根据实际情况自行调整大一些
enrichGO()中要设置OrgDb = "org.Mm.eg.db"；readable = TRUE表示将entrez ID转化为gene Symbol；ont 参数可以选择"BP", "MF" "CC" 或 "ALL"，表示对细胞组分（cellular component, CC）、分子功能（molecular function, MF）、生物过程（biological process, BP)或全部的基因通路进行富集，可以根据实际需求选择
enrichKEGG()要设置organism = "mmu"，由于其没有readable参数，因此之后还要用DOSE包的setReadable进行转化entrez ID为gene Symbol
#### KEGG、GO富集 ####kegg_enrich_results <- enrichKEGG(gene = gene_up_entrez, organism = "mmu", #物种人类 hsa 小鼠mmu pvalueCutoff = 0.05, qvalueCutoff = 0.2)kegg_enrich_results <- DOSE::setReadable(kegg_enrich_results, OrgDb="org.Mm.eg.db", keyType='ENTREZID')#ENTREZID to gene Symbolwrite.csv(kk_read@result,'KEGG_gene_up_enrichresults.csv') save(kegg_enrich_results, file ='KEGG_gene_up_enrichresults.Rdata')
go_enrich_results <- enrichGO(gene = gene_up_entrez, OrgDb = "org.Mm.eg.db", ont = "ALL" , #One of "BP", "MF" "CC" "ALL" pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE)write.csv(go_enrich_results@result, 'GO_gene_up_BP_enrichresults.csv') save(go_enrich_results, file ='GO_gene_up_enrichresults.Rdata')4. 富集结果可视化：pathview goplot barplot dotplot cnetplot emapplot treeplot heatplot upsetplot
在富集到通路后就要进行可视化展示了，参见clusterprofiler说明书📖 Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)，其中enrichplot包可以对富集结果进行超级丰富的可视化。
下面就来总结介绍一下clusterprofiler说明书介绍的各种富集结果可视化方法，其中4.1pathview 只针对KEGG， 4.2goplot只针对GO结果，其他可视化方法则对两者都通用
4.1 pathview ——KEGG通路可视化
将KEGG通路进行可视化一般有以下三种方法：
*使用函数 browseKEGG(enrich_results, select_pathway)进行网页查看相关通路，如
browseKEGG(kegg_enrich_results, 'mmu04310') #网页查看通路 网页版的pathview https://pathview.uncc.edu/guest-home，可以上传数据进行在线可视化
pathview包：pathview()函数中需要输入含log2FC信息的gene.data、pathway.id 和species物种信息，会生成含有基因上下调信息的基因通路图。
注意函数中有一个重要参数kegg.native ：若TRUE则输出完整gene pathway的png文件，若为FALSE则输出只含输入基因信息的pdf文件 。
参数limit可以对图例color bar范围（即log2FC范围）进行调整。
其他参数使用详见官方说明：pathview.pdf (bioconductor.org)，再推荐一篇简要中文教程：可视化kegg通路-pathview包 | KeepNotes blog (bioinfo-scrounger.com)下面我们选取富集排名第3（富集结果是按pvalue由小到大排序的）的"Wnt signaling pathway" 进行示范。
#构建含log2FC信息的genelist genelist <- as.numeric(DEG_DESeq2[,2]) names(genelist) <- row.names(DEG_DESeq2)# genelist ID转化genelist_entrez <- genelistnames(genelist_entrez) <- bitr(names(genelist), fromType="SYMBOL",toType="ENTREZID", OrgDb=OrgDb)[,2] genelist_entrez <- genelist_entrez[!is.na(names(genelist_entrez))]
##查看与选择所需通路kk_read@result$Description[1:10] #查看前10通路i=3 select_pathway <- kk_read@result$ID[i] #选择所需通路的ID号
pathview(gene.data = genelist_entrez, pathway.id = select_pathway, species = 'mmu' , # 人类hsa 小鼠mmu kegg.native = T,# TRUE输出完整pathway的png文件，F输出基因列表的pdf文件 new.signature = F, #pdf是否显示pathway标注 limit = list(gene=2.5, cpd=1) #图例color bar范围调整 )
pathview(gene.data = genelist_entrez, pathway.id = select_pathway, species = 'mmu' , # 人类hsa 小鼠mmu kegg.native = F,# TRUE输出完整pathway的png文件，F输出基因列表的pdf文件 new.signature = F, #pdf是否显示pathway标注 limit = list(gene=2.5, cpd=1) #图例color bar范围调整 参数kegg.native = T，所得如下：
参数kegg.native = F，所得如下：
4.2 goplot—— GO富集结果的有向无环图 directed acyclic graph
注意当enrichGO()的ont为'ALL'时不能运行该图，会报错。以下 goplot展示的是enrichGO()的ont='BP'时的go_enrich_results
### goplot : Gene Ontology is organized as a directed acyclic graph.有向无环图gop <- goplot(go_enrich_results, showCategory = 10)ggsave(gop, filename = "goplot.pdf", width=10, height=10)4.3 barplot与dotplot——最常用的可视化图形
barplot与dotplot展示富集通路的p.adj，generatio，count信息。
如果enrichGO的ont为'ALL'，barplot与dotplot还可以设置使不同ont项目通路分隔开展示
### barplotbarp <- barplot(go_enrich_results, font.size=14, showCategory=10)+ theme(plot.margin=unit(c(1,1,1,1),'lines')) #如果enrichGO的ont为'ALL'if (F) { barp <- barplot(go_enrich_results, split= "ONTOLOGY", font.size =14)+ facet_grid(ONTOLOGY~., scale="free")+ theme(plot.margin=unit(c(1,1,1,1),'lines')) #调整图形四周留白大小 }ggsave(barp,filename = paste0(fn,'.pdf'), width=10, height=10)
### dotplot dotp <- enrichplot::dotplot(go_enrich_results,font.size =14)+ theme(legend.key.size = unit(10, "pt"),#调整图例大小 plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小#如果enrichGO的ont为'ALL'if (F) { dotp <- enrichplot::dotplot(go_enrich_results,font.size =8,split = 'ONTOLOGY')+ facet_grid(ONTOLOGY~., scale="free")+ theme(legend.key.size = unit(10, "pt"),#调整图例大小 plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小 }ggsave(dotp,filename = paste0(fn,'.pdf'),width =10,height =10)4.4 cnetplot——Gene-Concept Network
cnetplot展示了各通路下的基因网络。绘制cnetplot有两种展现方式， 更改参数circular 为 F（默认）或T可以分别得到散布状和圈状分布的cnetplot；cnetplot还可以输入含log2FC信息的genelist ，会将log2FC信息展现在图中
### cnetplot: Gene-Concept Network#构建含log2FC信息的genelist genelist <- as.numeric(DEG_DESeq2[,2]) names(genelist) <- row.names(DEG_DESeq2)
cnetp1 <- cnetplot(go_enrich_results, foldChange = genelist, showCategory = 6, colorEdge = T, node_label = 'all', color_category ='steelblue')cnetp2 <- cnetplot(go_enrich_results, foldChange = genelist, showCategory = 6, node_label = 'gene', circular = T, colorEdge = T)ggsave(cnetp1,filename ='cnetplot.pdf', width =12,height =10)ggsave(cnetp2,filename = 'cnetplot_cir.pdf', width =15,height =10)4.5 emapplot ——Enrichment Map
emapplot富集图将被富集的术语组织成一个边缘连接重叠基因集的网络，相互重叠的基因集往往会聚集在一起，因此有助于我们识别功能模块。
注意使用emapplot前还需要用pairwise_termsim()处理富集结果
### emapplot :Enrichment Mappt <- pairwise_termsim(go_enrich_results)emapp <- emapplot(pt, showCategory = 30, node_label = 'category') # 'category', 'group', 'all' and 'none'.)ggsave(emapp,filename = 'emapplot.pdf',width =10,height =10)4.6 heatplot: Heatmap-like functional classification
与cnetplot展示内容类似，但是将不同富集通路的相同基因聚在了一起，有助于我们更好识别基因模块
## heatplot: Heatmap-like functional classification #构建含log2FC信息的genelist genelist <- as.numeric(DEG_DESeq2[,2]) names(genelist) <- row.names(DEG_DESeq2)
heatp <- heatplot(go_enrich_results, foldChange = genelist, showCategory = 5)ggsave(heatp, filename ='heatplot.pdf', width=20, height=10)4.7 upsetplot——不同富集通路间的重叠基因数量
upsetplot展示不同富集通路间的重叠基因数量。
## upsetplot # emphasizes the gene overlapping among different gene setsupsetp <- upsetplot(go_enrich_results, n = 10)+ theme(plot.margin=unit(c(1,1,1,1),'lines')) #调整图形四周留白大小ggsave(upsetp, filename = 'upsetplot.pdf', width=15, height=10)4.7 treeplot——集结果术语的层次聚类与高频词标记
treeplot对富集结果术语进行层次聚类, 并使用高频词标记，有助于我们从繁多的富集结果中快速提取有用关键信息。
## Tree plot #pt <- pairwise_termsim(go_enrich_results)treep <- treeplot(pt, showCategory = 30)ggsave(treep, filename = 'treeplot.pdf', width=18, height=10)GO、KEGG富集与可视化到此就结束啦
如果所得显著差异基因很少，或无法富集到有生物学意义的通路时又该怎么办呢，下节将介绍一种不依赖于差异基因筛选的富集分析方法——GSEA
参考资料
📖 Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)pathview.pdf (bioconductor.org)
可视化kegg通路-pathview包 | KeepNotes blog (bioinfo-scrounger.com)
GitHub - jmzeng1314/GEO
【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili
【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表 生信入门课-2022年5~6月场，你的生物信息学第一课