【原】复现Nature图表：分组富集分析条形图展示通路及基因

TS的美梦 2024-04-03

展开全文

❝
1、《KS科研分享与服务》公众号有QQ交流群，进入门槛是20元（完全是为了防止白嫖党，请理解），请考虑清楚。群里有免费推文的注释代码和示例数据（终身拥有），没有付费内容，群成员福利是购买单个付费内容半价！
2、《KS科研分享与服务》微信VIP群只针对购买打包代码的小伙伴（公众号所有付费内容合集）！微信群不是单独的，是对于打包的人答疑解惑和交流的平台、群成员专享视频教程，帖子提前发布，以及其他更多福利！
点击：→ 加入微信vip群：2024-2025《KS科研分享与服务》付费内容打包集合
3、需进QQ群或者打包代码入微信VIP的小伙伴请添加作者微信了解，请备注目的，除此之外请勿添加，谢谢！
详情请联系作者：
❞

最近很多博主捅了这篇Nature文章的窝了，它的富集分析展示很好，我们这里也复现一下（群主真容像彭于晏）。这个图展示的很有特色，不仅展示了通路，就和我们之前展示的一样（富集分析柱状图大集合：通路展示在柱子上）。重要的一点是还展示了通路基因，那么在实际应用中，可以展示自己关注的通路，基因可以展示那些重要的基因即可。不仅仅是单细胞数据富集分析的展示，其他的富集分析也是可以这样展示的，只需要整理成相应的作图数据格式即可！

（reference：CHIT1-positive microglia drive motor neuron ageing in the primate spinal cord）

我们的复现效果如下：

完整版代码已发布2024年微信VIP，微信群_聊天记录_文件，自行下载！

我们这里使用上下调基因进行演示，首先构建一下富集分析的数据：

Macrophage <- subset(human_data, celltype=='Macrophage')diff <- FindMarkers(Macrophage, ident.1 = "GM", ident.2 = "BM",                    group.by = "group", logfc.threshold = 0.25,min.pct = 0.25)
diff$gene <- rownames(diff)diff$group <- ""diff$group <- ifelse(diff$avg_log2FC>0,"up",'down')

使用clusterProfiler进行富集分析，其他工具或者网站富集分析的，整理数据作图即可：这里我们挑选上下调前几个terms展示，实际中可展示自己需要的terms。

#富集分析,我们这里就以KEGG为例子group <- data.frame(gene=diff$gene,group=diff$group)#分组情况#gene转化为IDGene_ID <- bitr(diff$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")#构建文件并分析data  <- merge(Gene_ID,group,by.x='SYMBOL',by.y='gene')diff_KEGG <- compareCluster(ENTREZID~group,                            data=data,                            fun = "enrichKEGG",#函数选择什么定义什么分析                            pAdjustMethod = "BH",                            pvalueCutoff = 0.01,                            qvalueCutoff = 0.01,                            organism= "hsa")#物种
#将gene ID转化为gene symboldiff_KEGG = setReadable(diff_KEGG,OrgDb = "org.Hs.eg.db", keyType = "ENTREZID")#获取富集分析表格文件diff_KEGG <- diff_KEGG@compareClusterResult
diff_KEGG <- diff_KEGG %>%   group_by(group) %>%   top_n(n = 5, wt = -qvalue)

ggplot作图：

#排序diff_KEGG$group <- factor(diff_KEGG$group, levels = c("up","down"))# 使用排序索引重新排列数据框diff_KEGG <- diff_KEGG[order(diff_KEGG$group), ]#terms因子顺序diff_KEGG$Description <- factor(diff_KEGG$Description, levels = diff_KEGG$Description)
#展示的基因，我们选择每个terms展示5个基因，实际情况可以展示自己关注的基因diff_KEGG$geneID  <- sapply(strsplit(diff_KEGG$geneID , "/"), function(x) paste(x[1:5], collapse = "/"))

ggplot(diff_KEGG, aes(x = -log10(qvalue), y = rev(Description), fill = group))+  geom_bar(stat = "identity", width = 0.5)+  geom_text(aes(x=0.1,y=rev(Description),label = Description),size=3.5, hjust =0)+  theme_classic()+  theme(axis.text.y = element_blank(),        axis.ticks.y = element_blank(),        axis.title.y = element_text(colour = 'black', size = 12),        axis.line = element_line(colour = 'black', linewidth =0.5),        axis.text.x = element_text(colour = 'black', size = 10),        axis.ticks.x = element_line(colour = 'black'),        axis.title.x = element_text(colour = 'black', size = 12),        legend.position = "none")+  scale_x_continuous(expand = c(0,0))+  scale_fill_manual(values = c("#CB5640","#65B0C6"))+  geom_text(data = diff_KEGG,            aes(x = 0.1, y = rev(Description), label = geneID, color = group),            size = 4,            fontface = 'italic',             hjust = 0,            vjust = 2.3)+  scale_color_manual(values = c("#CB5640","#65B0C6"))+  scale_y_discrete(expand = c(0.1,0))+  labs(title = "Enrichment of genes",       y=c("Down                                     Up"))