細(xì)胞類型注釋不同類別方法及常見的一些軟件或工具介紹

細(xì)胞類型注釋是scRNA-seq數(shù)據(jù)分析中的關(guān)鍵一環(huán)。傳統(tǒng)的手動(dòng)注釋方法通常包含兩個(gè)步驟：首先使用無(wú)監(jiān)督學(xué)習(xí)對(duì)細(xì)胞進(jìn)行聚類，然后依據(jù)聚類結(jié)果中典型的差異表達(dá)基因來(lái)將細(xì)胞類群注釋為不同的類型。這一過(guò)程十分繁瑣耗時(shí)，并且依賴先驗(yàn)知識(shí)，帶有一定的主觀傾向；相比之下，自動(dòng)注釋工具和方法具有快速簡(jiǎn)潔的優(yōu)點(diǎn)。scRNA-seq細(xì)胞的自動(dòng)注釋一般可概括為三種類型[1]，下圖展示了不同類別方法中常見的一些軟件或工具（圖1）。

圖1：scRNA-seq細(xì)胞注釋的典型策略及代表性方法

第一種方法是基于Marker基因的細(xì)胞自動(dòng)注釋，依據(jù)一些特定基因來(lái)對(duì)細(xì)胞進(jìn)行標(biāo)記，判定細(xì)胞類型，選定的Marker基因或基因集應(yīng)當(dāng)具有細(xì)胞類型特異性并且能夠穩(wěn)定表達(dá)的。因此，該方法非常依賴Marker的準(zhǔn)確性。獲取Marker基因的方法有很多種，通�？梢詮囊恍┕�共數(shù)據(jù)庫(kù)或者文獻(xiàn)中得到。CellMarker[2]和PanglaoDB[3]是較為常見的存儲(chǔ)人類和小鼠組織中多種細(xì)胞類型標(biāo)記物的在線資源。第二種方法是基于參考轉(zhuǎn)錄組數(shù)據(jù)集的細(xì)胞注釋。使用具有細(xì)胞類型標(biāo)記的數(shù)據(jù)集作為參考，通過(guò)關(guān)聯(lián)輸入數(shù)據(jù)的基因表達(dá)值，從而將參考數(shù)據(jù)中的細(xì)胞或簇標(biāo)簽轉(zhuǎn)移到具有類似基因表達(dá)譜的輸入數(shù)據(jù)中的未標(biāo)記細(xì)胞或簇。該方法較為依賴參考數(shù)據(jù)的來(lái)源和可靠性。最后，第三種方法利用了有監(jiān)督機(jī)器學(xué)習(xí)，將帶有細(xì)胞標(biāo)記的參考數(shù)據(jù)作為訓(xùn)練集，通過(guò)機(jī)器學(xué)習(xí)方法得到一個(gè)模型（分類器），然后用于預(yù)測(cè)待注釋細(xì)胞的細(xì)胞類型。同樣，該方法依賴參考數(shù)據(jù)的質(zhì)量。下面，我們從上述方法中各挑選一個(gè)工具進(jìn)行闡述。

 

01 CellAssign

CellAssign[4]是一種基于概率圖模型的方法，于2019年發(fā)表在Nature Methods上。它利用已有的細(xì)胞類型Marker基因的先驗(yàn)知識(shí)，通過(guò)推斷細(xì)胞類型的概率來(lái)將未知數(shù)據(jù)分配到不同的細(xì)胞類型（圖2）。具體來(lái)說(shuō)，CellAssign將每個(gè)細(xì)胞表示為一個(gè)隨機(jī)變量，該變量表示細(xì)胞屬于各細(xì)胞類型的概率。同時(shí)，它還考慮了批次和樣本效應(yīng)的影響，通過(guò)引入?yún)f(xié)變量來(lái)調(diào)整這些效應(yīng)。CellAssign還考慮了基因的表達(dá)變異性，并使用原始計(jì)數(shù)矩陣和負(fù)二項(xiàng)分布來(lái)建模。通過(guò)最大化后驗(yàn)概率估計(jì)，CellAssign可以推斷每個(gè)細(xì)胞屬于各個(gè)細(xì)胞類型的概率，并將細(xì)胞分配到具體的細(xì)胞類型中。與其他方法相比，CellAssign在準(zhǔn)確性和F1得分方面表現(xiàn)更好，并且能夠控制批次效應(yīng)和樣本效應(yīng)。

圖 2：CellAssign概覽。(a) 使用CellAssign進(jìn)行細(xì)胞注釋，每個(gè)細(xì)胞被概率地分配到給定的細(xì)胞類型，而不需要任何手動(dòng)注釋或干預(yù)，并且考慮到批次效應(yīng)的影響。(b) 圖像展示了包含隨機(jī)變量、數(shù)據(jù)、分布假設(shè)等在內(nèi)的細(xì)胞類型分配概率模型。(c) 在CellAssign概率模型中使用的隨機(jī)變量的描述，以及它們的先驗(yàn)分布。

◆ 軟件的下載可參考如下方式（任選其一）：

# 從GitHub下載

1. install.packages("tensorflow")

2. tensorflow::install_tensorflow(extra_packages='tensorflow-probability')

3. install.packages("devtools")

4. devtools::install_github("Irrationone/cellassign")

# 通過(guò)Conda下載

conda install -c conda-forge -c bioconda r-cellassig

 

◆ 下面以軟件自帶的示例數(shù)據(jù)進(jìn)行基礎(chǔ)用法演示：

# 加載R包

1. library(SingleCellExperiment)

2.library(cellassign)

3. library(scran)

# 加載內(nèi)置數(shù)據(jù)集

4. data(example_sce)

5. data(example_marker_mat)

# 計(jì)算細(xì)胞標(biāo)準(zhǔn)化因子（示例數(shù)據(jù)已經(jīng)包含該部分信息，如果是新數(shù)據(jù)，需要重新計(jì)算。推薦使用scran包中的computeSumFactors 函數(shù)。注意：在計(jì)算標(biāo)準(zhǔn)化因子之前最好不要對(duì)表達(dá)矩陣的基因進(jìn)行篩選

6. # example_sce <- computeSumFactors(sce)

7. s <- sizeFactors(example_sce)

# 開始運(yùn)行cellassign。這一步關(guān)鍵的地方在于，必須要保證輸入的表達(dá)矩陣只包含marker基因的部分

8. fit <- cellassign(exprs_obj =

example_sce[rownames(example_marker_mat),],

marker_gene_info = example_marker_mat,

s = s,

learning_rate = 1e-2,

shrinkage = TRUE,

verbose = FALSE)

 

02 SingleR

在一篇研究肺纖維化相關(guān)巨噬細(xì)胞亞型的文章中[5]，作者同時(shí)開發(fā)了一種單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的細(xì)胞類型注釋工具——SingleR。它通過(guò)將未知類型的單細(xì)胞的基因表達(dá)與參考數(shù)據(jù)集中已知細(xì)胞類型的基因表達(dá)進(jìn)行相關(guān)性分析，來(lái)確定每個(gè)單細(xì)胞的細(xì)胞類型。它的注釋過(guò)程包括以下幾個(gè)步驟：
1、針對(duì)每個(gè)單細(xì)胞獨(dú)立進(jìn)行注釋。首先，計(jì)算單細(xì)胞表達(dá)與參考數(shù)據(jù)集中每個(gè)樣本的Spearman相關(guān)系數(shù)。這個(gè)相關(guān)性分析僅針對(duì)參考數(shù)據(jù)集中的變異基因進(jìn)行；
2、根據(jù)參考數(shù)據(jù)集中的命名注釋，將每個(gè)未知細(xì)胞的多個(gè)相關(guān)系數(shù)整合，為了防止由于參考樣本的異質(zhì)性而導(dǎo)致錯(cuò)誤分類，SingleR使用相關(guān)系數(shù)的80%作為閾值；
3、SingleR重新運(yùn)行相關(guān)性分析，但僅針對(duì)上一步中排名靠前的細(xì)胞類型。并且這個(gè)分析僅針對(duì)囊括的細(xì)胞類型之間的變異基因進(jìn)行。最低相關(guān)性值的細(xì)胞類型被移除（或者比最高值低0.05以上），然后重復(fù)這個(gè)步驟，直到只剩下兩個(gè)細(xì)胞類型為止。
最后一次運(yùn)行后，對(duì)應(yīng)于最高值的細(xì)胞類型被分配給單個(gè)細(xì)胞（圖3）。

 

圖 3：SingleR原理示意圖

◆ SingleR包自帶一些參考數(shù)據(jù)集，詳細(xì)信息可通過(guò)官方網(wǎng)站查詢，下面對(duì)它的基礎(chǔ)用法進(jìn)行演示：

# 從Bioconductor下載

1. BiocManager::install("SingleR")

# 加載R包及測(cè)試數(shù)據(jù)

2. library(celldex)

3. library(SingleR)

4. library(scRNAseq)

5. hpca.se <- HumanPrimaryCellAtlasData()

6. hESCs <- LaMannoBrainData('human-es')

7. hESCs <- hESCs[,1:100]

# 注釋

8. pred.hesc <- SingleR(test = hESCs, ref = hpca.se, assay.type.test=1, labels = hpca.se$label.main)

 

 

03 CellTypist

CellTypist[6]是一個(gè)用于自動(dòng)注釋人類組織中免疫細(xì)胞的工具。它的原理是通過(guò)整合來(lái)自不同組織的細(xì)胞數(shù)據(jù)，并使用機(jī)器學(xué)習(xí)方法訓(xùn)練模型，實(shí)現(xiàn)對(duì)細(xì)胞類型的準(zhǔn)確分類（圖4）。研究人員收集了來(lái)自20個(gè)不同組織的細(xì)胞數(shù)據(jù)，并對(duì)這些數(shù)據(jù)進(jìn)行深度篩選和整合，以獲得兩個(gè)層次的細(xì)胞類型信息。然后，他們使用邏輯回歸和隨機(jī)梯度下降學(xué)習(xí)的方法訓(xùn)練了模型。模型的性能通過(guò)精確度、召回率和全局F1分?jǐn)?shù)進(jìn)行評(píng)估，結(jié)果顯示在高層次和低層次的細(xì)胞類型分類上都達(dá)到了約0.9的準(zhǔn)確性。CellTypist能夠識(shí)別出不同的細(xì)胞亞群，包括T細(xì)胞、B細(xì)胞、單核巨噬細(xì)胞等，此外，它還能夠自動(dòng)注釋細(xì)胞的詳細(xì)亞型，提供了對(duì)免疫細(xì)胞群體的深入了解。它的預(yù)測(cè)結(jié)果對(duì)不同數(shù)據(jù)集之間的差異具有魯棒性，包括基因表達(dá)稀疏性和批次效應(yīng)等。

圖 4：CellTypist原理示意圖。CellTypist的工作流程包括數(shù)據(jù)收集、處理、模型訓(xùn)練和細(xì)胞類型預(yù)測(cè)

CellTypist既可以作為一個(gè)注釋工具使用，也可以作為一個(gè)免疫細(xì)胞相關(guān)的數(shù)據(jù)庫(kù)（圖5）。在其官網(wǎng)上，我們可以通過(guò)Encyclopedia查詢特定細(xì)胞類型的相關(guān)Marker信息，也可以輸入特定基因來(lái)查看其在不同組織不同細(xì)胞類型中的表達(dá)情況（圖6a-c）。

圖 5：CellTypist官方網(wǎng)站示意圖

而作為工具使用時(shí)，CellTypist提供了多種運(yùn)行方法。你可以通過(guò)命令行工具進(jìn)行細(xì)胞注釋，可以在Python環(huán)境下使用，可以通過(guò)Docker運(yùn)行。此外，你也同樣可以在其官網(wǎng)進(jìn)行在線分析（圖6d）。下面我們介紹幾種簡(jiǎn)單的用法。

◆ 軟件的下載：

# pip方式下載

1. pip install celltypist

# conda方式下載

2. conda install -c bioconda -c conda-forge celltypist
◆ Python中使用：

# 模塊導(dǎo)入

1. import celltypist

2. from celltypist import models

# 下載（初次使用時(shí)）并導(dǎo)入模型

3. models.download_models()

4. model = models.Model.load(model = 'Immune_All_Low.pkl')

# 細(xì)胞注釋

# 注意1：輸入數(shù)據(jù)的類型比較寬松，可以是環(huán)境內(nèi)部變量，可以是外部文件；

# 可以是表達(dá)矩陣，也可以是AnnData對(duì)象；

# 注意2：當(dāng)輸入表達(dá)矩陣時(shí)，需要保持行為細(xì)胞，列為基因；可以通過(guò)添加參數(shù)”--transpose-input”轉(zhuǎn)置

5. predictions = celltypist.annotate(input_file, model = model)

◆ 命令行使用：

# 下載模型（初次使用時(shí)）

1. celltypist –update-models

# 細(xì)胞注釋

2. celltypist --indata /path/to/input/file --model Immune_All_Low.pkl --outdir /path/to/outdir

圖 6：CellTypist在線工具。(a) 細(xì)胞或基因檢索界面。(b)查詢特定細(xì)胞類型后相關(guān)信息展示。(c)查詢特定基因后熱圖展示。(d) CellTypist在線注釋工具界面展示。

盡管自動(dòng)注釋工具可以自動(dòng)預(yù)測(cè)scRNA-seq數(shù)據(jù)的細(xì)胞類型，但與傳統(tǒng)的基于Marker的方法一樣，它們?nèi)匀恍枰�關(guān)于細(xì)胞類型的一些先驗(yàn)知識(shí)。前文我們提到了自動(dòng)注釋工具的三種分類方法，但究其本質(zhì)，仍然都是從一個(gè)reference開始的，無(wú)論是Marker還是數(shù)據(jù)集。因此，盡管使用同一種方法，但如果使用不同的參考數(shù)據(jù)（訓(xùn)練模型、細(xì)胞特異性Marker等），仍然可能導(dǎo)致截然不同的注釋結(jié)果。所以，在初期的準(zhǔn)備階段，如何選擇參考數(shù)據(jù)需要格外謹(jǐn)慎。

另外，依據(jù)注釋分辨率水平的不同，這些工具還可以被進(jìn)一步分成cell-based和cluster-based，即注釋到單個(gè)細(xì)胞或注釋到簇。研究發(fā)現(xiàn)[7]，參考數(shù)據(jù)的組成也會(huì)顯著影響自動(dòng)注釋工具的性能。測(cè)試中，注釋結(jié)果的準(zhǔn)確率會(huì)隨著參考數(shù)據(jù)中每種細(xì)胞類型的數(shù)量的增加而逐漸提高。當(dāng)參考數(shù)據(jù)中包含的某一細(xì)胞類型數(shù)量較小時(shí)，使用注釋單個(gè)細(xì)胞的方法可能不再合適。

在一些基準(zhǔn)測(cè)試中，不同的研究給出了不同的結(jié)論。比如在Tamim等人的研究中[8]，作者使用了27套測(cè)試數(shù)據(jù)，基于準(zhǔn)確性、未分類細(xì)胞的百分比和計(jì)算時(shí)間等系統(tǒng)性評(píng)估了22種注釋工具或方法。結(jié)果表明，支持向量機(jī)（SVM）分類器的整體表現(xiàn)最好。在另幾項(xiàng)研究中[7,9]，作者都沒(méi)有明確表明哪種方法在所有評(píng)估指標(biāo)上表現(xiàn)最好，而是在不同的指標(biāo)下，不同方法具有不同的表現(xiàn)性能。因此，對(duì)于自動(dòng)注釋來(lái)說(shuō)，更重要的可能不是工具的選擇，而是參考數(shù)據(jù)的質(zhì)量以及與待注釋樣本的匹配性。

最后，Zoe A. Clarke等人在21年給出了一份關(guān)于scRNA-seq細(xì)胞注釋的建議[10]。如圖所示（圖7），對(duì)于一個(gè)新的數(shù)據(jù)來(lái)說(shuō)，細(xì)胞注釋應(yīng)當(dāng)按照如下步驟進(jìn)行。①自動(dòng)注釋：建議首先嘗試使用自動(dòng)注釋方法。自動(dòng)注釋方法可以快速進(jìn)行細(xì)胞注釋，并且在處理大規(guī)模數(shù)據(jù)集時(shí)效率高。自動(dòng)注釋的結(jié)果可能會(huì)受到所選擇的參考數(shù)據(jù)的影響，因此可以嘗試不同的Marker基因組合或數(shù)據(jù)集來(lái)優(yōu)化注釋結(jié)果；②手動(dòng)注釋：如果自動(dòng)注釋的結(jié)果不夠準(zhǔn)確或存在沖突/缺失的細(xì)胞標(biāo)簽，就需要進(jìn)行手動(dòng)注釋。手動(dòng)注釋需要人工檢查每個(gè)細(xì)胞的特征，并參考各種資源來(lái)確定其功能和細(xì)胞類型。通常，可以通過(guò)可視化每個(gè)已知Marker基因在二維數(shù)據(jù)圖上的表達(dá)情況來(lái)進(jìn)行手動(dòng)注釋。此外，還可以通過(guò)查閱文獻(xiàn)和挖掘已有的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)來(lái)尋找額外的Marker基因；③最后，細(xì)胞注釋的結(jié)果可以進(jìn)行獨(dú)立驗(yàn)證，例如新的驗(yàn)證實(shí)驗(yàn)，或與互補(bǔ)數(shù)據(jù)進(jìn)行比較，例如空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)。

圖 7：細(xì)胞注釋過(guò)程由三個(gè)主要步驟組成:自動(dòng)細(xì)胞注釋、手動(dòng)注釋和驗(yàn)證。

 

參考文獻(xiàn)：

[1] Su, M., Pan, T., Chen, QZ. et al. Data analysis guidelines for single-cell RNA-seq in biomedical studies and clinical applications. Military Med Res 9, 68 (2022).

[2] Hu C, Li T, Xu Y, Zhang X, Li F, Bai J, Chen J, Jiang W, Yang K, Ou Q, Li X, Wang P, Zhang Y. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Res. 2023 Jan 6;51(D1):D870-D876.

[3] Oscar Franzén, Li-Ming Gan, Johan L M Björkegren, PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data, Database, Volume 2019, 2019, baz046, doi:10.1093/database/baz046

[4] Zhang, A.W., O’Flanagan, C., Chavez, E.A. et al. Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling. Nat Methods 16, 1007–1015 (2019).

[5] Aran, D., Looney, A.P., Liu, L. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 20, 163–172 (2019).

[6] C. Domínguez Conde et al. ,Cross-tissue immune cell analysis reveals tissue-specific features in humans.Science 376,eabl5197(2022).

[7] Huang Q, Liu Y, Du Y, Garmire LX. Evaluation of Cell Type Annotation R Packages on Single-cell RNA-seq Data. Genomics Proteomics Bioinformatics. 2021 Apr;19(2):267-281.

[8] Abdelaal T, Michielsen L, Cats D, Hoogduin D, Mei H, Reinders MJT, Mahfouz A. A comparison of automatic cell identification methods for single-cell RNA sequencing data. Genome Biol. 2019 Sep 9;20(1):194.

[9] Xie B, Jiang Q, Mora A, Li X. Automatic cell type identification methods for single-cell RNA sequencing. Comput Struct Biotechnol J. 2021 Oct 20;19:5874-5887.

[10] Clarke, Z.A., Andrews, T.S., Atif, J. et al. Tutorial: guidelines for annotating single-cell transcriptomic maps using automated and manual methods. Nat Protoc 16, 2749–2764 (2021).