前面给大家介绍过神兵利器——单细胞细胞类群基因marker鉴定新方法：COSG，现在给大家分享一个新的方法：genesorteR。并且会对这两个方法进行一个准确性和时间成本两个方面的简单比较。
首先在学习一个工具之前，我个人觉得还是得多多少少对它的原理进行一些了解，所以首先从我的角度来理解一下这个工具的原理，当然这部分是选择性的去进行。
genesorteR原理概述（optional reading）
首先我们要明白：一个好的基因marker无非具有的特征是其表达具有细胞类型特异性。正是基于这个标准，genesorteR对一个基因细胞类型表达特异性的水平进行了如下刻画：
首先，genesorteR和其它工具一样都是假定你已经聚好类，分好群。细胞共被分为k类（type），分别记为：
共有j个基因（gene），分别记为：
那么首先，当我们观测到某个基因表达时，这些表达的细胞中细胞类型
出现的概率就是：
使用条件概率的定义做如下变换：
对这里几个部分的含义进行分析：首先表示的是在细胞类型
中基因
表达的比例；
就是细胞类型为
的细胞的比例；
表示在所有细胞中检测到基因
表达的比例。所以综上，
这个参数就刻画了每个基因表达细胞类型特异性程度，显然其表达的细胞类型特异性程度越高，
就会在相应的那个细胞类型中越高。
进一步的，作者将这个参数进行了如下的变换，最终得到了：
最终，**这个参数将是一个
之间的值，越接近于1，就表明基因
越特异性的表达于细胞类群（类型）
。
更进一步，为了综合考量每个基因在不同细胞类型中的表达特异性程度，作者提出了对于每个基因基于熵的参数。（注意不是
哦）
所以一个标准的marker gene应该具备两个特质：（1）有一个比较高的
（在某一种细胞类型里面表达较普遍）；（2）有比较高的
（其余的细胞类型表达较少）。
实例代码分析
我们还是使用pbmc公共数据集来进行测试：
#install packages
devtools::install_github("mahmoudibrahim/genesorteR")
remotes::install_github(repo = 'genecell/COSGR')
library(Seurat)
library(SeuratData)
library(genesorteR)
library(COSG)
data("pbmc3k.final")
#identify gene marker with COSG
starttime <- Sys.time()
COSG_markers <- cosg(
pbmc3k.final,
groups = 'all',
assay = 'RNA',
slot = 'data',
mu = 1)
endtime <- Sys.time()
timecost <- endtime - starttime
print(timecost)
#Time difference of 0.3216481 secs
#identify gene marker with genesorteR
starttime <- Sys.time()
gs <- sortGenes(x = pbmc3k.final@assays$RNA@data,
classLabels = Idents(pbmc3k.final))
head(gs$specScore)
genesorteR_markers <- getMarkers(gs = gs, quant = 0.99)
pp <- plotMarkerHeat(gs$inputMat,
gs$inputClass,
genesorteR_markers$markers,
clusterGenes = T,
outs = T)
pp$gene_class_info
endtime <- Sys.time()
timecost <- endtime - starttime
print(timecost)
#Time difference of 1.358346 secs
几个函数和参数做个说明：
-
sortGenes函数就是在计算每个基因在不同的细胞类型中的 -
getMarkers函数选取了位于前1%（quant = 0.99）的值进行分类，最终鉴定出真正的marker基因。
最终genesorteR鉴定出来的marker：
S100A9 S100A8 FCER1A FCGR3A GNLY SDPR TMEM40 GP9 SPON2 FGFBP2
3 3 1 9 2 6 8 6 2 2
GZMK GZMA LTB PTCRA GNG11 FCN1 PRF1 IFITM3 MS4A7 MS4A1
4 4 7 8 6 3 2 9 9 5
AP001189.4 CD3D LDHB GZMB IL32 CCL5 ITGA2B CD79B CST3 CST7
6 7 7 2 7 4 6 5 3 4
CD79A NKG7 SEPT5 LGALS2
5 4 6 3
检查一下效果：
library(patchwork)
library(ggplot2)
#for cosg
p1 <- DotPlot(pbmc3k.final,
assay = 'RNA',
features = COSG_markers$names$B[1:10]) +
ggtitle(label = 'B cell Markers with COSG') +
theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- DotPlot(pbmc3k.final,
assay = 'RNA',
features = which(pp$gene_class_info == 5) %>% names()) +
ggtitle(label = 'B cell Markers with genesorteR') +
theme(axis.text.x = element_text(hjust = 1, angle = 45))
p1+p2
使用体验
