最近很多博主捅了这篇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的，可以用到自己文章中，希望分享对你有用！
