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

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

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

我们的复现效果如下：

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

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转化为ID
Gene_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 symbol
diff_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"))

这个图还是非常nice的，可以用到自己文章中，希望分享对你有用！