ixxmu
commented
2 months ago genesorteR——比FindAllMarkers快一百倍！ by Biomamba 生信基地
### 读入测试数据 ###
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 SeuratObjectscRNA <- readRDS('pbmcrenamed.rds')
DimPlot(scRNA)这是一个已经注释过的数据，这时通常我们会利用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: usethisif(!require(genesorteR))devtools::install_github("mahmoudibrahim/genesorteR")## Loading required package: genesorteR## Loading required package: Matrixsystem.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 EClibrary(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计算结果")) 可以看出，面对相同的数据，
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.欢迎在发文/毕业时向我们分享你的喜悦~
认真学完下面的内容，3~10分的SCI理应是囊中之物(当然也不要真的去水文章)
单细胞系列教程目录
如何联系我们
https://mp.weixin.qq.com/s/U_F6dfyK4oK8V4AWF9RQkQ