查看原文
其他

RNA-seq入门实战(六):GO、KEGG富集分析与enrichplot超全可视化攻略

生信技能树 生信技能树 2022-08-10


连续两次求贤令:曾经我给你带来了十万用户,但现在祝你倒闭,以及 生信技能树知识整理实习生招募,让我走大运结识了几位优秀小伙伴!大家开始根据我的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,多一点数据认知,让他们的科研上一个台阶:


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存