不知道大家有没有在使用Seurat的FindAllMarkers时常常觉得计算时间过长，人生苦短，等待的时间真的很煎熬。genesorteR是一个用于单细胞数据分析的R软件包。它计算一个特异性分数来对每个细胞簇中的所有基因进行排序。然后，它可以使用这个排序来找到一组标记基因，或者找到高度可变或差异表达的基因。也就是说，genesorteR可以很好的替代FindAllMarkers，并且秒出结果！genesorteR既适用于scRNA-Seq数据，也适用于scATAC-Seq等稀疏单细胞数据，我们简单演示一下。测试数据大家自己扫码获取吧：
1.png
通过百度网盘分享的文件：genesorteR链接：https://pan.baidu.com/s/1h7CZ2MqXiGWc-nIYbnCcBQ?pwd=cn1m 提取码：cn1m
2.jpg
### 读入测试数据 ###library(Seurat)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,## which was just loaded, were retired in October 2023.## Please refer to R-spatial evolution reports for details, especially## https://r-spatial.org/r/2023/05/15/evolution4.html.## It may be desirable to make the sf package available;## package maintainers should consider adding sf to Suggests:.
## Attaching SeuratObject
scRNA <- readRDS('pbmcrenamed.rds')DimPlot(scRNA)
1.png
这是一个已经注释过的数据，这时通常我们会利用Seurat内置的FindAllMarkers()来计算各细胞类型的marker基因，我们来看一下计算时间：
system.time({seu_marker <- FindAllMarkers(scRNA)})
## Calculating cluster VSMC
## Calculating cluster T cell
## Calculating cluster Macro
## Calculating cluster B cell
## Calculating cluster Fibro
## Calculating cluster Myeloid cells
## Calculating cluster Mono
## Calculating cluster Neut
## Calculating cluster EC
## user system elapsed ## 163.181 0.967 164.148
可以看出计算marker需要花费大量的时间，这里给大家介绍一个加速marker计算的新方法：
# 装包：if(!require(devtools))install.packages("devtools")
## Loading required package: devtools
## Loading required package: usethis
if(!require(genesorteR))devtools::install_github("mahmoudibrahim/genesorteR")
## Loading required package: genesorteR
## Loading required package: Matrix
system.time( { # 计算： gs <- sortGenes(scRNA@assays$RNA@data, Idents(scRNA)) # 获取marker列表： gs_marker <- getMarkers(gs, quant = 0.99) })
## Warning in sortGenes(scRNA@assays$RNA@data, Idents(scRNA)): A Friendly Warning:## Some genes were removed because they were zeros in all cells after## binarization. You probably don't need to do anything but you might want to look## into this. Maybe you forgot to pre-filter the genes? You can also use a## different binarization method. Excluded genes are available in the output under## '$removed'.
## 'as(<dgeMatrix>, "dgCMatrix")' is deprecated.## Use 'as(., "CsparseMatrix")' instead.## See help("Deprecated") and help("Matrix-deprecated").
## user system elapsed ## 1.357 0.108 1.466
可以说计算速度有质的飞跃。
我们对比一下两种方式筛选出的top5marker：
library(dplyr)
## ## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':## ## filter, lag
## The following objects are masked from 'package:base':## ## intersect, setdiff, setequal, union
# 查看Seurat筛选出的top5 marker基因seu_marker_top5 <- seu_marker %>% group_by(cluster) %>% slice_max(n = 5, order_by = avg_log2FC)# 将genesorteR计算的top5基因整理为数据框spec_scores <- gs$specScore# 获取基因名genes <- rownames(spec_scores)# 创建一个空的数据框来存储结果top5_genes_df <- data.frame()# 遍历每个 cluster (specScore 的列)for (cluster in colnames(spec_scores)) { # 获取当前 cluster 的 specScore 列并绑定基因名 cluster_scores <- data.frame( gene = genes, specScore = spec_scores[, cluster], cluster = cluster ) # 按 specScore 降序排列并取前5个基因 top5_genes <- cluster_scores %>% arrange(desc(specScore)) %>% head(5) # 将结果合并到总的数据框中 top5_genes_df <- rbind(top5_genes_df, top5_genes)}# 打印最终的 top 5 基因数据框top5_genes_df
## gene specScore cluster## Itga8 Itga8 0.6400709 VSMC## Npnt Npnt 0.6387903 VSMC## Ppp1r14a Ppp1r14a 0.5918881 VSMC## Ccn3 Ccn3 0.5916800 VSMC## Lmod1 Lmod1 0.5912409 VSMC## Cd3g Cd3g 0.7788462 T cell## Ms4a4b Ms4a4b 0.7617241 T cell## Cd3d Cd3d 0.7570093 T cell## Trbc2 Trbc2 0.7180583 T cell## Skap1 Skap1 0.6912000 T cell## Clec4a3 Clec4a3 0.5825532 Macro## Ms4a6c Ms4a6c 0.5259740 Macro## Clec4a1 Clec4a1 0.5159770 Macro## Sirpb1c Sirpb1c 0.4932967 Macro## Ccr2 Ccr2 0.4490722 Macro## Cd79a Cd79a 0.8695652 B cell## Ms4a1 Ms4a1 0.8549495 B cell## Ly6d Ly6d 0.8335849 B cell## Iglc3 Iglc3 0.8227174 B cell## Iglc2 Iglc2 0.8151579 B cell## Serpinf1 Serpinf1 0.7032967 Fibro## Lum Lum 0.6982759 Fibro## Dpt Dpt 0.6812162 Fibro## Dcn Dcn 0.6097561 Fibro## Pcolce Pcolce 0.6073494 Fibro## Alox12 Alox12 0.6377953 Myeloid cells## Itga2b Itga2b 0.6267361 Myeloid cells## Clec1b Clec1b 0.6203165 Myeloid cells## Treml1 Treml1 0.6150893 Myeloid cells## Gp5 Gp5 0.5632184 Myeloid cells## C1qc C1qc 0.7455682 Mono## C1qa C1qa 0.7342353 Mono## C3ar1 C3ar1 0.7101370 Mono## C1qb C1qb 0.7077895 Mono## Tnf Tnf 0.7004167 Mono## Mmp9 Mmp9 0.7680000 Neut## Retnlg Retnlg 0.7458678 Neut## Slfn4 Slfn4 0.7348544 Neut## Ly6g Ly6g 0.7272727 Neut## Cxcr2 Cxcr2 0.7054839 Neut## Egfl7 Egfl7 0.9209184 EC## Ptprb Ptprb 0.8709677 EC## Cdh5 Cdh5 0.8501149 EC## Ecscr Ecscr 0.7810714 EC## Mmrn2 Mmrn2 0.7401316 EC
library(ggplot2)library(patchwork)(VlnPlot(scRNA,features = seu_marker_top5$gene,stack = T)+ ggtitle("Seurat计算结果"))/ (VlnPlot(scRNA,features = top5_genes_df$gene,stack = T) + ggtitle("genesorteR计算结果"))
1.png
可以看出，面对相同的数据，genesorteR不仅计算速度更快，marker的特异性也更高一些。
参考：https://doi.org/10.1101/676379
最后、如果以上内容对你有帮助，欢迎在文章的Acknowledgement中加上这一段：
Since Biomamba and his wechat public account team produce bioinformatics tutorials and share code with annotation, we thank Biomamba for their guidance in bioinformatics and data analysis for the current study.
欢迎在发文/毕业时向我们分享你的喜悦~
