在單細(xì)胞數(shù)據(jù)分析當(dāng)中，當(dāng)我們聚類分群完成之后，緊接著就是細(xì)胞類型注釋，細(xì)胞類型的注釋就離不開 基因marker ，即在目標(biāo)細(xì)胞類群和其它細(xì)胞類群之間呈現(xiàn)出不同表達(dá)模式特征的基因，這樣我們就能根據(jù)這些marker對(duì)照現(xiàn)存的marker list或相關(guān)數(shù)據(jù)庫(kù)來(lái)做我們自己的細(xì)胞類型注釋。
目前關(guān)于細(xì)胞類群基因marker鑒定的方法最常用的當(dāng)屬 Seurat 官方的 FindMarkers與FindAllMarkers函數(shù)，在這里給大家分享的是最新online的工具：COSG。
工具來(lái)源
COSG工具于2022年3月在線發(fā)表于 Brief Bioinformatics 上，題為：Accurate and fast cell marker gene identification with COSG，該工具同時(shí)有Python和R版本，這意味著無(wú)論你是 Seurat 用戶還是 Scanpy 用戶都可以使用這個(gè)工具。
原理概述（僅供參考）
傳統(tǒng)的基因marker鑒定方法基于歐式距離，即目標(biāo)細(xì)胞類群與其它細(xì)胞之間的基因表達(dá)差異來(lái)進(jìn)行marker的鑒定。但這樣的鑒定方法可能存在潛在的問(wèn)題：當(dāng)基因marker在目標(biāo)細(xì)胞類群以外的其它細(xì)胞中的一小部分高表達(dá)時(shí)（如某個(gè)小的細(xì)胞亞型），傳統(tǒng)的鑒定方法仍然會(huì)將其鑒定為基因marker，因?yàn)樾〔糠值母弑磉_(dá)對(duì)于整體的細(xì)胞群來(lái)說(shuō)并不具備決定性意義，仍然會(huì)認(rèn)為這個(gè)基因marker在其余的細(xì)胞群中是低表達(dá)的。
那么為了解決這個(gè)潛在的問(wèn)題，中科院遺傳所的Xiu-Jie Wang老師們開發(fā)出了基于 余弦值 的基因marker鑒定新方法：COSG（COSine similarity-based marker Gene identification）。為什么余弦值就能解決上面的這個(gè)問(wèn)題？答案就在于余弦值不依賴于向量的模，在單細(xì)胞分析的背景下就是不依賴于基因的表達(dá)量，而依賴于基因的表達(dá)模式。舉個(gè)簡(jiǎn)單的例子：基因A在細(xì)胞1和2中的表達(dá)量分別是2和4，基因B在細(xì)胞1和2中的表達(dá)量分別是1和2，如果根據(jù)歐式距離的算法，顯然基因A和基因B會(huì)被定義成為不同的表達(dá)模式，但是如果我們基于余弦值的話，兩個(gè)基因在該細(xì)胞集合中的表達(dá)向量之間的夾角就為：
在這個(gè)背景下，兩個(gè)基因表達(dá)向量之間的夾角為0，也就是這兩個(gè)基因在這個(gè)細(xì)胞集合中的表達(dá)模式實(shí)際上是一樣的。
基于上面的原理，COSG 是怎么鑒定細(xì)胞類群的基因marker的呢？總結(jié)為以下幾步：
- 第一步：基于現(xiàn)有的分群情況，COSG首先對(duì)每個(gè)細(xì)胞類群鑒定出一個(gè) artificial gene，這個(gè)基因的表達(dá)特征是：只在目標(biāo)細(xì)胞類群中表達(dá)，且不在其它任何一個(gè)細(xì)胞類群中有表達(dá)，這個(gè)基因就是每個(gè)細(xì)胞類群最理想的基因marker了；
- 第二步：假設(shè)一共有k個(gè)細(xì)胞，那么每個(gè)基因的表達(dá)情況就是一個(gè) k維的向量 （在每個(gè)細(xì)胞中的表達(dá)量作為一個(gè)維度），那么對(duì)于每個(gè)基因和每個(gè)細(xì)胞類群，COSG會(huì)做如下操作：首先，計(jì)算該基因在目標(biāo)細(xì)胞類群中與該目標(biāo)類群artificial gene的表達(dá)向量之間的夾角；再計(jì)算該基因在其它細(xì)胞類群中與其它細(xì)胞類群的artificial gene的表達(dá)向量之間的夾角。最終鑒定出來(lái)的目標(biāo)細(xì)胞類群的基因marker應(yīng)該有如下特征：與目標(biāo)細(xì)胞類群的artificial gene表達(dá)向量之間的夾角越小越好（即有相似的表達(dá)模式）而與其它細(xì)胞類群的artificial gene表達(dá)向量之間的夾角越大越好（即有相反的表達(dá)模式）。
這樣，COSG就完成了細(xì)胞類群基因marker的鑒定了。
代碼
-
第一步：安裝COSG R包。
# install.packages('remotes')
remotes::install_github(repo = 'genecell/COSGR')
-
第二步：基因marker鑒定
library(Seurat)
library(SeuratData)
library(DT)
#加載內(nèi)置數(shù)據(jù)集
data("pbmc3k")
pbmc3k <- pbmc3k.final
invisible(gc())
#查看細(xì)胞類群
celltype <- table(pbmc3k@meta.data$seurat_annotations)
datatable(as.data.frame(celltype))
#經(jīng)典的FindAllMarkers()函數(shù)
starttime <- Sys.time()
Seurat_markers <- FindAllMarkers(pbmc3k, only.pos = TRUE)
endtime <- Sys.time()
timecost <- endtime - starttime
print(timecost)
最終耗時(shí)：35.88439s。
再使用COSG來(lái)試試：為了驗(yàn)證COSG的效率，這里每個(gè)cluster統(tǒng)一返回900個(gè)marker（相比較于Seurat返回的數(shù)量最大值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)
最終耗時(shí)：1.838692s。
再來(lái)看看效果（以B細(xì)胞為例）：
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
感覺結(jié)果真的還不錯(cuò)哦~
最后感謝曾健明老師及時(shí)與我們分享這個(gè)新的工具！
