在单细胞的数据分析当中，每个亚群的top基因是十分重要的，因为这一部分的基因主要是代表了这一亚群的高表达基因，为了后面的分群鉴定，主要是通过seurat的findallmarkers这个函数进行计算。可以参考这个博主的文章，对源码解析的很细https://www.jianshu.com/p/f5c8f9ea84af，同时对应着这个函数的解析http://www.idata8.com/rpackage/Seurat/FindAllMarkers.html。
top50markergene的提取
为了加大后面的鉴定结果，我主要提取了每个亚群与其他亚群相比差别较大的前50的基因，来做后面的富集分析。
代码语言：javascript复制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表上进行操作。
同时为了提高后面在富集分析后的工作效率，提前将每个基因的注释结果也读到新的那个表格里面，然后进行整理，有利于后面在翻看注释文件的时候，节省时间。
将所需要的两个表格进行整合后。
代码语言：javascript复制##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文件的语句，然后进行更改。
循环读入每个亚群的结果
代码语言：txt复制##首先写一个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那里的方法进行其他的修改。
循环后的文件夹结果总结
主要是需要先把自己要做富集分析的cluster读到R中，然后进行循环语句的读写，R中的循环语句主要注意的是自己用的是什么数据，需要怎么读入文件中。目前是批量完了，还没有报错，做完了，可以跟公司的结果进行对比，查看数据质量的重复性。
