【单细胞高级绘图】07.KEGG富集结果展示

这一节画的图是比较新的，图中我用红色箭头标出的是pathway一级注释信息(big annotation，自己想的，非专有名词)，纵轴花花绿绿的标注是pathway的二级注释(small annotation)。「如何获取注释」算一个难点，我上一讲也已经讲过：KEGG通路的从属/注释信息如何获取。

整个图反映的是有多少基因落到了对应的分类里面。

「辩证地看」，整张图都是pathway注释，没有具体的pathway名称，跟平常做的富集分析很不一样。把图里面的二级注释换成具体的pathway会更好。另外，这个图中横坐标是基因数，但我觉得在富集分析中这个数值并不重要（基因ratio比单个数值重要），我们可以换成p值。

---

在开始画图之前，需要整理一个表格，表格中至少包含：「pathway ID、pathway description、pathway注释/big annotation、几个富集指标、被富集到的基因」。

• 「前三个信息」，上一讲的表格已经整理好了，可以拿来用；
• 「几个富集指标」在 clusterprofiler的输出结果中也有；
• 「被富集到的基因」貌似clusterprofiler做kegg富集不能直接显示基因symbol，所以返回的结果需要稍微加工一下。

---

整理表格的代码名称为run.R，如下所示：

library(Seurat)
library(tidyverse)
library(xlsx)

testseu=readRDS("testseu.rds")
Idents(testseu)="anno_new"

### 找差异基因 #########################################################################
marker_celltype=FindAllMarkers(testseu,logfc.threshold = 0.8,only.pos = T)
# 过滤
marker_celltype=marker_celltype%>%filter(p_val_adj < 0.01)
marker_celltype$d=marker_celltype$pct.1-marker_celltype$pct.2
marker_celltype=marker_celltype%>%filter(d > 0.2)
marker_celltype=marker_celltype%>%arrange(cluster,desc(avg_log2FC))
marker_celltype=as.data.frame(marker_celltype)
write.xlsx(marker_celltype,file = "markers_log2fc0.8_padj0.01_pctd0.2.xlsx",row.names = F,col.names = T)

### 富集分析 ###########################################################################
library(clusterProfiler)
library(org.Hs.eg.db)
R.utils::setOption("clusterProfiler.download.method","auto") #https://github.com/YuLab-SMU/clusterProfiler/issues/256

source("syEnrich.R")
syEnrich(marker_celltype,outpath = "markers_log2fc0.8_padj0.01_pctd0.2")

### 挑一类细胞来作为演示 #######################################################
kegg.res=read.xlsx("markers_log2fc0.8_padj0.01_pctd0.2.KEGG.xls",sheetIndex = 1,as.data.frame = T,header = T)
kegg.res=kegg.res%>%filter(p.adjust < 0.05)
kegg.res=kegg.res%>%filter(cluster == "Endothelial")

# 导入上一讲的文件
kegg_info=read.xlsx("kegg_info.xlsx",sheetIndex = 1,startRow = 3)
kegg_info=kegg_info[,c("ID","Pathway","big.annotion")]

# 合并两个表格
kegg.res$ID=str_replace(kegg.res$ID,"hsa","")
kegg.res=kegg.res%>%inner_join(kegg_info,by = "ID")
write.table(kegg.res,file = "kegg.res.txt",quote = F,sep = "\t",row.names = F,col.names = T)

我以单细胞分析中的kegg富集分析作为演示，只取其中一个cluster的富集结果来画图。

上述代码中间用到的富集代码叫syEnrich.R，这个文件只需要输入单细胞seurat对象运行FindAllMarkers得到的差异基因，就可以返回GO/KEGG富集结果，同时被富集到某个通路的基因symbol也会被列出。

运行run.R之后，最终的表格如下图所示：

---

然后开始画图，代码名称为3.plot.R，这里就不演示了，最终可以得到的图如下：

获取代码

包含这张图会用到的所有代码，数据整理以及画图，超贴心有没有！

这个系列都会采取「限时公开」的方式共享代码，24小时内是免费的。超过这个时间如何获取，后台回复2022A可知（可能需要你动动小手转发一下）。

代码和测试数据的网盘链接如下: 链接：https://pan.baidu.com/s/1hebbeQH4DgYA8JqFWXO4dg 提取码：abo5

觉得代码有用的话，可以给个三/二/一连（文末点赞、分享、点下小广告） 帮帮忙好嘛:(前几次涨粉少得可怜