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

本节概览：
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(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")
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)
kk_read <- 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') #网页查看通路

网页版的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 <- genelist
names(genelist_entrez) <- bitr(names(genelist),
                               fromType="SYMBOL",toType="ENTREZID", 
                               OrgDb="org.Mm.eg.db")[,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项目通路分隔开展示

### 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)

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 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)

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 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)

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)

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

最后编辑于：2022.07.26 23:19:31

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 194,319评论 5赞 459
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 81,801评论 2赞 371
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 141,567评论 0赞 319
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,156评论 1赞 263
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,019评论 4赞 355
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,090评论 1赞 272
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,500评论 3赞 381
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,192评论 0赞 253
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,474评论 1赞 290
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,566评论 2赞 309
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,338评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,212评论 3赞 312
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,572评论 3赞 298
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,890评论 0赞 17
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,169评论 1赞 250
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,478评论 2赞 341
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,661评论 2赞 335