[富集分析] 3、clusterprofile富集分析--下游可视化

1、富集分析的背景知识 - 简书 (jianshu.com)
2、clusterprofile富集分析--上游分析 - 简书 (jianshu.com)
3、clusterprofile富集分析--下游可视化 - 简书 (jianshu.com)
4、clusterprofile富集分析--多组设计的富集及可视化 - 简书 (jianshu.com)

富集分析结果的可视化主要依赖enrichplot包,可分为三类:(1)单纯展示富集最显著的通路;(2)富集通路与差异基因的网络关系;(3)富集通路之间的网络关系。由于GSEA方法的独特性,可以对Enrichment score(正负),Running score进行可视化。

0、示例数据

library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
library(DOSE)
library(ReactomePA)
data(geneList)
gene <- names(geneList)[abs(geneList) > 2]
#ORA
res1 <- enrichPathway(gene)
res1 <- setReadable(res1, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
res2 <- enrichGO(gene, OrgDb = org.Hs.eg.db, ont = "ALL")
res2 <- setReadable(res2, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
#GSEA
res3 <- gseKEGG(geneList     = geneList,
                organism     = 'hsa')
res3 <- setReadable(res3, OrgDb = org.Hs.eg.db, keyType="ENTREZID")

1、展示被富集的通路

1.1 barplot

(1)经典出图
  • 其实这里的barplot()并非enrichplot包的函数,就是一个基础绘图函数。
  • 如下图:从高到低,按p值升序排列;横轴标示差异基因与通路基因的重复基因数。
  • barplot只支持ORA富集分析结果的可视化。
barplot(res1, showCategory=20) 

如果barplot()结果返回如下结果报错,很有可能是因为根据既有的过滤标准没有被富集到的通路,放宽限定的阈值即可。例如pvalueCutoff = 1, qvalueCutoff = 1
此外可通过设置label_format=参数,设置轴标签文字的每行字符数长度,过长则会自动换行。

Error in ans[ypos] <- rep(yes, length.out = len)[ypos] : 
  replacement has length zero
In addition: Warning message:
In rep(yes, length.out = len) : 'x' is NULL so the result will be NULL
(2)仅展示感兴趣的相关通路
df = res1@result
# 根据dataframe结果,选择感兴趣的通路
interesting_set = c("R-HSA-141444","R-HSA-174048", "R-HSA-176412")
# 找到Description对应的ID
index= which(rownames(res1@result) %in% interesting_set)
Description=res1@result$Description[index]
barplot(res1, showCategory=Description) 
(3)针对分类关系的基因集,可分别展示
table(res2@result$ONTOLOGY)
# BP  CC  MF 
# 178  31  21
barplot(res2, split = "ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale = "free")
(4)针对GSEA富集分析结果,存在Enrichment score正负(上下调)的关系,可以通过ggplot2包可视化这些关系
library(ggplot2)
library(forcats)
library(ggstance)
ggplot(res3, aes(NES, fct_reorder(Description, NES), fill=qvalues), showCategory=20) + 
  geom_bar(stat='identity') + 
  scale_fill_continuous(low='red', high='blue', guide=guide_colorbar(reverse=TRUE)) + 
  theme_minimal() + ylab(NULL)
(5)简洁图--only P值
  • 通路 + 显著程度
df = res1@result[1:10,]
library(ggplot2)
library(ggstatsplot)
df$logP = -log10(df$p.adjust)
df=df[order(df$logP,decreasing = T),]
#限定通路名的每行长度
df$nice_description = stringr::str_wrap(gsub("_"," ",df$Description), width = 30)
# reorder(Description, logP)绘图时,对Description按照logP升序排列
# reorder(Description, logP)按照logP降序排列
ggplot(df, aes(x=reorder(nice_description, logP), y=logP)) + 
  geom_bar(stat="identity", fill="#7fc97f") + 
  xlab("Pathway names") +
  ylab("-log10(P.adj)") +
  scale_y_continuous(expand=c(0,0)) + coord_flip() + 
  theme_ggstatsplot() +
  theme(plot.title = element_text(size = 15,hjust = 0.5), #调整title的大小、位置  
        axis.text = element_text(size = 12,face = 'bold'), #调整轴标题的大小、位置
        axis.ticks.y=element_blank(), #不显示刻度
        panel.grid = element_blank()) + #空白背景
  ggtitle("Enrichment Barplot") 
  • 通路 + 显著程度 + 上下调信息(GSEA)
#选择感兴趣的上下调通路,如下选择的是上下调的Top5(ES)
df = res3@result
df = df[order(df$enrichmentScore),]
#设定分组:1--ES>0;   -1--ES<0
up = head(subset(df, enrichmentScore>0),5);up$group=1
down = tail(subset(df, enrichmentScore<0),5);down$group=-1
dat=rbind(up,down)
#dat$group = factor(dat$group);str(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group 
dat=dat[order(dat$pvalue,decreasing = F),]
ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) + 
  geom_bar(stat="identity", aes(fill=factor(group, levels = c(1,-1),  
                                            labels = c("ES>0","ES<0")))) + 
  xlab("Pathway names") +
  ylab("-log10(P.adj)") +
  coord_flip() + 
  theme_ggstatsplot() +
  scale_y_continuous(breaks=c(-4, -2, 0, 2, 4),
                     labels=c("4", "2", "0","2","4")) +
  scale_fill_manual(values = c("#ff6633","#34bfb5")) + 
  theme(plot.title = element_text(size = 15,hjust = 0.5),  
        axis.text = element_text(size = 12,face = 'bold'),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        legend.title = element_blank()) +
  ggtitle("Pathway Enrichment") 

1.2 dotplot

  • dotplot就是enrichplot的函数了,适合展示ORA的富集结果
  • 如下图所示,dotplot相比barplot增添了size,可增添一个维度的信息。
  • 默认color→p值,dot size→counts(差异基因与通路基因集重复数), GeneRatio→比例(counts/通路基因集大小)。我们可以设置参数修改这些映射关系。
dotplot(res3, split = "ONTOLOGY") 
  • 类似dotplot的棒棒糖图
library(ggplot2)
library(forcats)
ggplot(res1, showCategory = 20, 
       aes(GeneRatio, fct_reorder(Description, GeneRatio))) + 
  geom_segment(aes(xend=0, yend = Description)) +
  geom_point(aes(color=p.adjust, size = Count)) +
  scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
  scale_size_continuous(range=c(2, 10)) +
  theme_minimal() + 
  ylab(NULL) 

2、富集通路与差异基因的网络关系图

2.1 cnetplot

可设置的参数很多,常用的有

  • showCategory :设置展示几个通路及相关差异基因,也可以为字符串,指定感兴趣的通路
  • node_label:设置是否展示通路名/基因名 (“category,” “gene,” “all” and “none”)
  • layout:设置排布样式,e.g. 'star', 'circle', 'gem', 'dh', 'graphopt', 'grid', 'mds', 'randomly', 'fr', 'kk', 'drl' or 'lgl'
library(patchwork)
p1 = cnetplot(res1, foldChange=geneList,showCategory = 3)
p2 = cnetplot(res1, foldChange=geneList,showCategory = 3, circular = T)
p3 = cnetplot(res1, foldChange=geneList,showCategory = 3, colorEdge = T, node_label="gene")
p4 = cnetplot(res1, foldChange=geneList,showCategory = 3, layout = "gem")
(p1 + p2) / (p3 + p4)

由于需要具体展示基因,所以不太适合GSEA,因为往往GSEA会包含过多的基因,不适合展示。

2.2 heatplot

  • 类似上面的cnetplot,只是变成了热图的形式展示。
heatplot(res1, foldChange=geneList)

3、富集通路之间的网络关系图

  • 对富集的通路,两两间计算相似性,绘制通路网络图
res1 <- enrichplot::pairwise_termsim(res1)
emapplot(res1, showCategory = 15, layout="kk", cex_category=1.5,min_edge = 0.8) 
res3 <- enrichplot::pairwise_termsim(res3)
emapplot(res3, showCategory = 15, layout="kk", min_edge = 0.1, color = "NES") 

4、GSEA running score

4.1 gseaplot

如下图结果,上半部分表示位于Cell cycle基因集中的差异基因的差异倍数;下半部分表示根据这些基因的打分过程(简单理解:从左到右遍历时,遇到差异基因就加分,不是就减分)

  • 设置参数by = "preranked"只展示上半部分,by = "runningScore"只展示下半部分。
gseaplot(res3, geneSetID = 1, title = res3$Description[1])

4.2 gseaplot2

相比gseaplot,gseaplot2可一次展示多个基因集的富集结果

gseaplot2(res3, geneSetID = 1:3, pvalue_table = T)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,511评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,495评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,595评论 0 225
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,558评论 0 190
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,715评论 3 270
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,672评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,112评论 2 291
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,837评论 0 181
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,417评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,928评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,316评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,773评论 2 234
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,253评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,827评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,440评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,523评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,583评论 2 249

推荐阅读更多精彩内容