RNA-seq入门实战（六）：GO、KEGG富集分析与enric
本节概览：
1.获取DEG结果的上下调差异基因
2.bitr()函数转化基因名为entrez ID
3.利用clusterProfiler进行KEGG与GO富集
4.用enrichplot进行富集结果可视化：pathview goplot barplot dotplot cnetplot emapplot treeplot heatplot upsetplot
承接上期文章：RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较，本节介绍使用差异基因进行GO、KEGG富集分析与enrichplot可视化
1. 获取DEG结果的上下调差异基因
载入上节RNA-seq入门的简单实战（五）中保存的三种差异分析结果数据，这里示范选取DESeq2的结果数据，进行筛选条件设置后获取上下调基因名
rm(list=ls())
options(stringsAsFactors = F)
library(clusterProfiler)
library(enrichplot)
library(tidyverse)
library(ggstatsplot)
library(ggnewscale)
# 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.001
padj_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 Symbol
write.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') #网页查看通路
image.png
-
网页版的pathview https://pathview.uncc.edu/guest-home，可以上传数据进行在线可视化
image.png
-
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 <- genelist
names(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，所得如下：
image.png
参数kegg.native = F，所得如下：
image.png
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)
image.png
4.3 barplot与dotplot——最常用的可视化图形
barplot与dotplot展示富集通路的p.adj，generatio，count信息。
如果enrichGO的ont为'ALL'，barplot与dotplot还可以设置使不同ont项目通路分隔开展示
### barplot
barp <- 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)
image.png
image.png
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)
image.png
image.png
4.5 emapplot ——Enrichment Map
emapplot富集图将被富集的术语组织成一个边缘连接重叠基因集的网络，相互重叠的基因集往往会聚集在一起，因此有助于我们识别功能模块。
注意使用emapplot前还需要用pairwise_termsim()处理富集结果
### emapplot :Enrichment Map
pt <- 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)
image.png
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)
image.png
4.7 upsetplot——不同富集通路间的重叠基因数量
upsetplot展示不同富集通路间的重叠基因数量。
## upsetplot # emphasizes the gene overlapping among different gene sets
upsetp <- 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)
image.png
4.7 treeplot——集结果术语的层次聚类与高频词标记
treeplot对富集结果术语进行层次聚类, 并使用高频词标记，有助于我们从繁多的富集结果中快速提取有用关键信息。
使用 treeplot也需要用pairwise_termsim()处理富集结果。
## Tree plot #
pt <- pairwise_termsim(go_enrich_results)
treep <- treeplot(pt,
showCategory = 30)
ggsave(treep, filename = 'treeplot.pdf', width=18, height=10)
image.png
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
RNA-seq实战系列文章：( 2022.5.21 更新)
RNA-seq入门实战（零）：RNA-seq流程前的准备——Linux与R的环境创建 - 简书 (jianshu.com)
RNA-seq入门实战（一）：上游数据下载、格式转化和质控清洗
RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon
RNA-seq入门实战（三）：从featureCounts与Salmon输出文件获取counts矩阵
RNA-seq入门实战（四）：差异分析前的准备——数据检查
RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较
RNA-seq入门实战（六）：GO、KEGG富集分析与enrichplot超全可视化攻略
RNA-seq入门实战（七）：GSEA——基因集富集分析
RNA-seq入门实战（八）：GSVA——基因集变异分析
未完待续......
