拿到一个多样品的单细胞转录组数据后，我们通常是使用批量读取每个单细胞转录组样品数据，然后harmony进行多样品整合，然后走seurat的标准单细胞降维聚类分群流程，这样的话会得到几十个单细胞亚群，它们是顺序编号的，但是我们第一层次通常是根据自己背诵好的基因把几十个单细胞亚群合并成为十个左右的生物学功能亚群。
但是， 最开始我们得到几十个单细胞亚群的时候，就需要对每个亚群找一下各自的单细胞亚群特异性高表达基因，通常是使用Seurat包的FindAllMarkers函数，这个函数的帮助文档写的是：Finds markers (differentially expressed genes) for each of the identity classes in a dataset ，默认使用 Wilcoxon Rank Sum test (default) 方法。
虽然FindAllMarkers是目前大家使用频率最高的函数，也有其它 ‘bimod’, ‘MAST’, ‘negbinom’, ‘poisson’ or ‘roc’ 等参数，但是速度都很慢！尤其是细胞数量超十万的时候，运行大半天是很正常的事情。最近看到交流群有小伙伴 推荐了一个解决方案，就是 COSG (COSine similarity-based marker Gene identification), 目前发表在 2022. Accurate and fast cell marker gene identification with COSG. Brief. Bioinform. bbab579.
文章写的是哪怕是100万个单细胞的数据集，也只需要两分钟左右，就可以快速找到各个单细胞亚群特异性高表达基因。它既有Python代码也有R版本分发：
- https://github.com/genecell/COSG (Python)
- https://github.com/genecell/COSGR (R).
测试数据
这里 我们仍然是选择 seurat-data包的测试数据pbmc3k哦：
代码语言：javascript复制# install.packages('devtools')
# devtools::install_github('satijalab/seurat-data')
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
DimPlot(sce,label = T,repel = T)
可以看到它是如下所示降维聚类分群：
第一层次降维聚类分群
这些信息都是数据集提供者已经做好了的，你也可以看看我们前面的例子：人人都能学会的单细胞聚类分群注释 ， 你必须要理解这个例子里面的降维聚类分群和注释。
Seurat包的FindAllMarkers函数
为了加快速度，我这里采取了4个线程（因为数据量很小，所以线程多一点也不怕，如果你是大数据而计算机配置很差，前往不要加线程！！！）
代码语言：javascript复制library(future)
# check the current active plan
plan()
plan("multiprocess", workers = 4)
plan()
start = Sys.time()
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
end = Sys.time()
dur = end-start
dur
可以看到是 Time difference of 1.608828 mins ，可以看到针对全部的9个单细胞亚群，都找到了符合阈值的单细胞亚群特异性高表达基因列表：
FindAllMarkers函数结果
很常规的结果啦，相信大家都很熟悉了，这7列信息。每个单细胞亚群都有几百个左右的基因被挑选出来，9个单细胞亚群就是2885个基因，当然了，每个基因的特异性不一样哦。
接下来是COSG
代码也超级简单的，因为就是cosg 函数替代了前面的 Seurat包的FindAllMarkers函数而已：
代码语言：javascript复制#remotes::install_github(repo = 'genecell/COSGR')
start = Sys.time()
library(COSG)
marker_cosg <- cosg(
sce,
groups='all',
assay='RNA',
slot='data',
mu=1,
n_genes_user=100)
end = Sys.time()
dur = end-start
dur
可以看到， 每个单细胞亚群默认就找到了100个单细胞亚群特异性高表达基因：
默认就找到了100个单细胞亚群特异性
可以看到， 这个结果跟前面的Seurat包的FindAllMarkers函数结果的形式就不一样的，前面的是长数据，现在这个 cosg 函数的结果是宽数据，有点类似于前面我们提到的：基于非负矩阵分解的单细胞降维聚类分群，找各个亚群的特征基因 。
其实 特异性的基因我们一般来说看十个就绰绰有余，要太多也没有必要。
简单的比较两个函数的结果：
首先看看 cosg 函数的结果是否都在前面的 Seurat包的FindAllMarkers函数结果里面：
代码语言：javascript复制findG = split(sce.markers$gene,sce.markers$cluster)
names(findG)
tmp = do.call(rbind ,
lapply(names(findG), function(x){
table(marker_cosg$names[,x] %in% findG[[x]])
}))
rownames(tmp) = names(findG)
tmp
如下所示：
代码语言：javascript复制> tmp
FALSE TRUE
Naive CD4 T 19 81
Memory CD4 T 37 63
CD14 Mono 15 85
B 60 40
CD8 T 71 29
FCGR3A Mono 32 68
NK 56 44
DC 80 20
Platelet 37 63
可以看到 DC和CD8 T 的表现比较差。因为前面提到了 特异性的基因我们一般来说看十个就绰绰有余，要太多也没有必要。所以我们试试看仅仅是看top10的基因，它们交集如何：
代码语言：javascript复制Naive CD4 T 7
Memory CD4 T 8
CD14 Mono 10
B 10
CD8 T 9
FCGR3A Mono 10
NK 10
DC 5
Platelet 10
这个时候可以看到仍然是DC表现很差，那么为什么 cosg 函数针对DC细胞亚群找到的特异性基因结果在前面的 Seurat包的FindAllMarkers函数结果里面超过半数都没有呢？
那，我们检测一下两个函数的各自的top10基因看看，代码是：
代码语言：javascript复制library(stringr)
library(ggplot2)
g1 =head( marker_cosg$names$DC ,10)
g2 = head(findG$DC,10)
th= theme(axis.text.x = element_text(angle = 90))
p1 = DotPlot(sce, features = g1,assay='RNA' ) th NoLegend()
p2 = DotPlot(sce, features = g2,assay='RNA' ) th
library(patchwork)
p1 p2
蛮有意思的，Seurat包的FindAllMarkers函数效果要好一点， cosg 函数虽然在速度方面有飞跃提升，也可以得到绝大部分细胞亚群的重要的的基因，但是在部分细胞数量比较少的单细胞亚群里面会出现一些不完善的结果，如下所示：
两个函数的效果对比
左边是 cosg 函数针对DC细胞亚群 的top10基因，右边是Seurat包的FindAllMarkers函数。
