单细胞数据分析-R语言对分群结果的top基因循环做富集分析

原创

小胡子刺猬的生信学习123

发布于 2022-05-20 18:16:41

2.1K0

发布于 2022-05-20 18:16:41

在单细胞的数据分析当中，每个亚群的top基因是十分重要的，因为这一部分的基因主要是代表了这一亚群的高表达基因，为了后面的分群鉴定，主要是通过seurat的findallmarkers这个函数进行计算。可以参考这个博主的文章，对源码解析的很细https://www.jianshu.com/p/f5c8f9ea84af，同时对应着这个函数的解析http://www.idata8.com/rpackage/Seurat/FindAllMarkers.html。

top50markergene的提取

为了加大后面的鉴定结果，我主要提取了每个亚群与其他亚群相比差别较大的前50的基因，来做后面的富集分析。

library(Seurat)
markers <- FindAllMarkers(object = obj, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.36)
top = 50 # can be adjusted as needed
TopMarkers <- markers %>% filter(p_val_adj < 0.01) %>% group_by(cluster) %>%
top_n(n = top, wt = avg_log2FC)
write.csv(TopMarkers,'CellType/TopMarkers.csv', row.names=F)

得到了每个亚群的marker基因，准备来循环做GO和KEGG的富集分析。

目前主要的问题是我们是在一张表里面有每个亚群的基因，所以需要首先将每个亚群的基因循环读到一个文件，然后在将基因的ID进行转换，然后进行富集分析。

每个亚群的基因循环读进一个文件

首先是需要对每个基因的id进行转换，我主要是用的clusterProfiler进行转换，但是我做的这个物种的id信息我是主要在phytozome上下载的，然后用的OrgDb的加载文件是在ncbi上下的，所以两个数据库的id号不同，我需要先在biodbnet进行全部的转换，读到一个新的表格里面，然后在进行转换，我这里主要是用的最近新学的dplyr包里面的函数，大家可以看一下这个博主的文章:https://blog.csdn.net/weixin_42437924/article/details/108701670，里面写的很详细，在前面进行两个数据集取交集的时候，大大的节省了时间，不需要在excel表上进行操作。

先尝试进行单个亚群的富集分析

首先是需要对每个基因的id进行转换，我主要是用的clusterProfiler进行转换，但是我做的这个物种的id信息我是主要在phytozome上下载的，然后用的OrgDb的加载文件是在ncbi上下的，所以两个数据库的id号不同，我需要先在biodbnet进行全部的转换，读到一个新的表格里面，然后在进行转换，我这里主要是用的最近新学的dplyr包里面的函数，大家可以看一下这个博主的文章:https://blog.csdn.net/weixin_42437924/article/details/108701670，里面写的很详细，在前面进行两个数据集取交集的时候，大大的节省了时间，不需要在excel表上进行操作。

同时为了提高后面在富集分析后的工作效率，提前将每个基因的注释结果也读到新的那个表格里面，然后进行整理，有利于后面在翻看注释文件的时候，节省时间。

将所需要的两个表格进行整合后。

##TopMarkers为前面获得的每个亚群的top50的高表达的基因，ann为自己手动整理的注释及基因转换id的文件，将TopMarkers的geneid为标准，进行取交集，获得TopMarkers里面基因的注释结果和geneid号
library(dplyr)
locus <- left_join(TopMarkers, ann, by = "geneID")
##尝试提取每个cluster的基因号，分别输出一个数据框结果，然后准备前面的单个成功后，尝试做循环将所有的亚群进行富集分析
cluster0 <- (which(pro[,3] == 0))
cluster_0 <- pro[cluster0,]
cluster_0<-na.omit(cluster_0)
ego = enrichGO(cluster_0$ENTREZID, OrgDb = ara, keyType="ENTREZID", pvalueCutoff=1, qvalueCutoff=1, ont='all')
dotplot(ego, split="ONTOLOGY",showCategory = 20,label_format=100) + facet_grid(ONTOLOGY~., scale="free")
compKEGG = enrichKEGG(cluster_0$ENTREZID, organism="ath", pvalueCutoff=1, keyType='kegg')
dotplot(compKEGG, showCategory = 15,label_format=100, title = "KEGG Pathway Enrichment Analysis")

通过目前的尝试，以上的代码没有发生报错的现象，因此我目前开始准备写循环，进行亚群的批量富集分析。主要也是参考我前几次肺癌文章里面的批量读取cellranger的gz文件的语句，然后进行更改。

循环读入每个亚群的结果

##首先写一个xsl的文件，将cluster读进去，这里如果亚群数目少，可以选择第2种方法，这里可以参照以前教程里面的excle的表格的模板
library(readxl)
cluster <- read_excel("name.xls", range = cell_cols("A:A")) %>% .$cluster_id
##首先进行go富集，保存csv文件和dot的图
##第2种方法：for (i in c(1, 2, 3,, 4, 5, 6,7))
for (i in seq_along(cluster)){
 ego = enrichGO(gene = eval(parse(text = paste0("cluster", i)))$ENTREZID, OrgDb = gly, keyType="ENTREZID", pvalueCutoff=1, qvalueCutoff=1, ont='all')
 print(dotplot(ego, split="ONTOLOGY",showCategory = 20,label_format=100) + facet_grid(ONTOLOGY~., scale="free")) 
 ggsave(paste0("cluster_", i, "_go.png"), path = "F:/", width = 30, height = 45, units = "cm")
 write.csv(ego, paste0("F:/cluster_", i, "_go.csv"),  quote = FALSE, row.names = FALSE)
}
##然后进行kegg富集分析
for (i in seq_along(cluster)){
 compKEGG = enrichKEGG(eval(parse(text = paste0("cluster", i)))$ENTREZID, organism="gmx", pvalueCutoff=1, keyType='kegg')
 print(dotplot(compKEGG, showCategory = 20,label_format=100, title = "KEGG Pathway Enrichment Analysis"))
 ggsave(paste0("cluster_", i, "_kegg.png"), path = "F:/", width = 15, height = 20, units = "cm")
 write.csv(compKEGG, paste0("F:/cluster_", i, "_kegg.csv"),  quote = FALSE, row.names = FALSE)
}

如果自己需要其他的图片可以参照print那里的方法进行其他的修改。