神兵利器——单细胞细胞类群基因marker鉴定新方法：COSG

在单细胞数据分析当中，当我们聚类分群完成之后，紧接着就是细胞类型注释，细胞类型的注释就离不开 基因marker ，即在目标细胞类群和其它细胞类群之间呈现出不同表达模式特征的基因，这样我们就能根据这些marker对照现存的marker list或相关数据库来做我们自己的细胞类型注释。

目前关于细胞类群基因marker鉴定的方法最常用的当属 Seurat 官方的 FindMarkers与FindAllMarkers函数，在这里给大家分享的是最新online的工具：COSG。

工具来源

COSG工具于2022年3月在线发表于 Brief Bioinformatics 上，题为：Accurate and fast cell marker gene identification with COSG，该工具同时有Python和R版本，这意味着无论你是 Seurat 用户还是 Scanpy 用户都可以使用这个工具。

原理概述（仅供参考）

传统的基因marker鉴定方法基于欧式距离，即目标细胞类群与其它细胞之间的基因表达差异来进行marker的鉴定。但这样的鉴定方法可能存在潜在的问题：当基因marker在目标细胞类群以外的其它细胞中的一小部分高表达时（如某个小的细胞亚型），传统的鉴定方法仍然会将其鉴定为基因marker，因为小部分的高表达对于整体的细胞群来说并不具备决定性意义，仍然会认为这个基因marker在其余的细胞群中是低表达的。

那么为了解决这个潜在的问题，中科院遗传所的Xiu-Jie Wang老师们开发出了基于 余弦值 的基因marker鉴定新方法：COSG（COSine similarity-based marker Gene identification）。为什么余弦值就能解决上面的这个问题？答案就在于余弦值不依赖于向量的模，在单细胞分析的背景下就是不依赖于基因的表达量，而依赖于基因的表达模式。举个简单的例子：基因A在细胞1和2中的表达量分别是2和4，基因B在细胞1和2中的表达量分别是1和2，如果根据欧式距离的算法，显然基因A和基因B会被定义成为不同的表达模式，但是如果我们基于余弦值的话，两个基因在该细胞集合中的表达向量之间的夹角就为：

在这个背景下，两个基因表达向量之间的夹角为0，也就是这两个基因在这个细胞集合中的表达模式实际上是一样的。

基于上面的原理，COSG 是怎么鉴定细胞类群的基因marker的呢？总结为以下几步：

• 第一步：基于现有的分群情况，COSG首先对每个细胞类群鉴定出一个 artificial gene，这个基因的表达特征是：只在目标细胞类群中表达，且不在其它任何一个细胞类群中有表达，这个基因就是每个细胞类群最理想的基因marker了；
• 第二步：假设一共有k个细胞，那么每个基因的表达情况就是一个 k维的向量 （在每个细胞中的表达量作为一个维度），那么对于每个基因和每个细胞类群，COSG会做如下操作：首先，计算该基因在目标细胞类群中与该目标类群artificial gene的表达向量之间的夹角；再计算该基因在其它细胞类群中与其它细胞类群的artificial gene的表达向量之间的夹角。最终鉴定出来的目标细胞类群的基因marker应该有如下特征：与目标细胞类群的artificial gene表达向量之间的夹角越小越好（即有相似的表达模式）而与其它细胞类群的artificial gene表达向量之间的夹角越大越好（即有相反的表达模式）。

这样，COSG就完成了细胞类群基因marker的鉴定了。

代码

• 第一步：安装COSG R包。

# install.packages('remotes')
remotes::install_github(repo = 'genecell/COSGR')

• 第二步：基因marker鉴定

library(Seurat)
library(SeuratData)
library(DT)

#加载内置数据集
data("pbmc3k")
pbmc3k <- pbmc3k.final

invisible(gc())

#查看细胞类群
celltype <- table(pbmc3k@meta.data$seurat_annotations)
datatable(as.data.frame(celltype))

#经典的FindAllMarkers()函数
starttime <- Sys.time()
Seurat_markers <- FindAllMarkers(pbmc3k, only.pos = TRUE)
endtime <- Sys.time()
timecost <- endtime - starttime
print(timecost)

最终耗时：35.88439s。
再使用COSG来试试：为了验证COSG的效率，这里每个cluster统一返回900个marker（相比较于Seurat返回的数量最大值800多）。

library(COSG)
library(dplyr)
Seurat_markers %>% count(cluster)
starttime <- Sys.time()
COSG_markers <- cosg(
  pbmc3k,
  groups='all',
  assay='RNA',
  slot='data',
  mu=1,
  n_genes_user=900)
endtime <- Sys.time()
timecost <- endtime - starttime
print(timecost)

最终耗时：1.838692s。

再来看看效果（以B细胞为例）：

library(patchwork)
library(ggplot2)
p1 <- DotPlot(pbmc3k, 
              assay = 'RNA',
              features = subset(Seurat_markers, cluster == 'B')$gene[1:10]) + 
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) + 
  ggtitle(label = 'B cell Markers with Seurat') +
  NoLegend()
p2 <- DotPlot(pbmc3k, 
              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))
p1 + p2

感觉结果真的还不错哦~

最后感谢曾健明老师及时与我们分享这个新的工具！