黄思源

V1

2022/07/25阅读:41主题:全栈蓝

【单细胞高级绘图】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

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

分类:

其他

标签:

其他

作者介绍

黄思源
V1