[富集分析] 3、clusterprofile富集分析--下游可
2021-08-15 本文已影响0人
小贝学生信
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)
