### 读入测试数据 ###
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)

这是一个已经注释过的数据，这时通常我们会利用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计算结果"))

 可以看出，面对相同的数据，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.

欢迎在发文/毕业时向我们分享你的喜悦~