9.4 自動(dòng)注釋

9.4.1 總論

這一篇討論的方法將是自動(dòng)化的方法频丘，而不是手動(dòng)注釋數(shù)據(jù)。與上一篇展示的方法不同，這些方法中的每一種都可以使您以自動(dòng)化的方式對(duì)數(shù)據(jù)進(jìn)行注釋励烦。它們基于不同的原則贝或，有時(shí)需要預(yù)先定義的標(biāo)記集合吼过，或者在預(yù)先存在的完整scRNA-seq數(shù)據(jù)集上訓(xùn)練。
如前所述咪奖，自動(dòng)生成的注釋的質(zhì)量可能會(huì)有所不同盗忱。更具體地說(shuō)，注釋的質(zhì)量取決于：

• 選擇的分類(lèi)器類(lèi)型：之前的基準(zhǔn)研究表明羊赵，不同類(lèi)型的分類(lèi)器通常表現(xiàn)相當(dāng)趟佃，基于神經(jīng)網(wǎng)絡(luò)的方法通常不會(huì)優(yōu)于通用模型如支持向量機(jī)或線性回歸模型。
• 訓(xùn)練分類(lèi)器的數(shù)據(jù)的質(zhì)量昧捷。如果訓(xùn)練數(shù)據(jù)沒(méi)有很好地注釋或以低分辨率注釋闲昭，分類(lèi)器也會(huì)做同樣的事情。同樣靡挥，如果訓(xùn)練數(shù)據(jù)和/或其注釋有噪聲序矩，分類(lèi)器可能表現(xiàn)不佳。
• 自己的數(shù)據(jù)與訓(xùn)練分類(lèi)器的數(shù)據(jù)的相似度跋破。例如簸淀，如果分類(lèi)器是在drop-seq單細(xì)胞數(shù)據(jù)集上進(jìn)行訓(xùn)練的，并且數(shù)據(jù)是10X單核而不是單細(xì)胞drop-seq毒返，則這可能會(huì)降低注釋的質(zhì)量租幕。在包括多種數(shù)據(jù)集的跨數(shù)據(jù)集圖集上訓(xùn)練的分類(lèi)器可能比在單個(gè)數(shù)據(jù)集上訓(xùn)練的分類(lèi)器提供更穩(wěn)健和更好質(zhì)量的注釋。

上述幾點(diǎn)強(qiáng)調(diào)了使用分類(lèi)器可能存在的缺點(diǎn)饿悬，具體取決于訓(xùn)練數(shù)據(jù)和模型類(lèi)型令蛉。盡管如此，使用預(yù)先訓(xùn)練的分類(lèi)器來(lái)注釋數(shù)據(jù)有幾個(gè)重要的優(yōu)點(diǎn)。首先珠叔，它是一種快速蝎宇、簡(jiǎn)單的數(shù)據(jù)注釋方法。注釋不需要下載或預(yù)處理訓(xùn)練數(shù)據(jù)祷安，有時(shí)只需將數(shù)據(jù)上傳到在線網(wǎng)頁(yè)姥芥。其次，這些方法不像手動(dòng)注釋那樣依賴于將數(shù)據(jù)劃分為集群汇鞭。第三凉唐，預(yù)先訓(xùn)練的分類(lèi)器使您能夠直接利用以前研究中的知識(shí)和信息，例如高質(zhì)量的注釋霍骄。
最后台囱，由于這些分類(lèi)器的透明性往往不如基于人工標(biāo)記的標(biāo)注，一個(gè)好的不確定性度量量化標(biāo)注將提高方法的質(zhì)量和可用性读整。我們將在下面更廣泛地討論這個(gè)問(wèn)題簿训。

9.4.2 基于marker基因的分類(lèi)器

一類(lèi)自動(dòng)化細(xì)胞類(lèi)型注釋方法依賴于一組預(yù)定義的標(biāo)記基因。根據(jù)這些標(biāo)記基因的表達(dá)水平米间，細(xì)胞被分為細(xì)胞類(lèi)型强品。這些模型所基于的標(biāo)記基因集越穩(wěn)健、越通用屈糊，模型的性能就越好的榛。然而，與其他模型一樣逻锐，它們可能會(huì)受到模型訓(xùn)練數(shù)據(jù)和需要標(biāo)記的數(shù)據(jù)之間與批次效應(yīng)相關(guān)的差異的影響夫晌。與基于較大基因集的模型（見(jiàn)下文）相比，這些方法的優(yōu)點(diǎn)之一是它們更加透明：我們知道根據(jù)哪些基因進(jìn)行分類(lèi)谦去。

9.4.3 基于更廣泛的基因集的分類(lèi)器

值得注意的是慷丽，到目前為止討論的方法僅使用數(shù)據(jù)中檢測(cè)到的基因的一小部分：通常每種細(xì)胞類(lèi)型僅使用一組1至10個(gè)標(biāo)記基因。另一種方法是使用一個(gè)分類(lèi)器鳄哭，該分類(lèi)器將更大的基因集（數(shù)千個(gè)或更多）作為輸入要糊，從而更多地利用scRNA-seq數(shù)據(jù)的廣度。此類(lèi)分類(lèi)器是在先前注釋的數(shù)據(jù)集或圖譜上進(jìn)行訓(xùn)練的妆丘。其中的示例包括 CellTypist（另請(qǐng)參閱 https://www.celltypist.org锄俄，其中數(shù)據(jù)可以上傳到門(mén)戶以獲取自動(dòng)細(xì)胞注釋）和 Clustifyr。
讓我們?cè)跀?shù)據(jù)上嘗試使用CellTypist勺拣。根據(jù)CellTypist教程(https://www.celltypist.org/tutorials)奶赠，我們需要準(zhǔn)備數(shù)據(jù)，以便將計(jì)數(shù)標(biāo)準(zhǔn)化為每個(gè)細(xì)胞10,000個(gè)計(jì)數(shù)药有，然后進(jìn)行l(wèi)og1p轉(zhuǎn)換：

adata_celltypist = adata.copy()  # make a copy of our adata
adata_celltypist.X = adata.layers["counts"]  # set adata.X to raw counts
sc.pp.normalize_per_cell(
    adata_celltypist, counts_per_cell_after=10**4
)  # normalize to 10,000 counts per cell
sc.pp.log1p(adata_celltypist)  # log-transform
# make .X dense instead of sparse, for compatibility with celltypist:
adata_celltypist.X = adata_celltypist.X.toarray()

我們現(xiàn)在將下載免疫細(xì)胞的celltypist模型：

models.download_models(
    force_update=True, model=["Immune_All_Low.pkl", "Immune_All_High.pkl"]
)

輸出結(jié)果：

?? Retrieving model list from server https://celltypist.cog.sanger.ac.uk/models/models.json
?? Total models in list: 19
?? Storing models in /home/icb/lisa.sikkema/.celltypist/data/models
?? Total models to download: 2
?? Downloading model [1/2]: Immune_All_Low.pkl
?? Downloading model [2/2]: Immune_All_High.pkl

讓我們嘗試一下Immune_All_Low和Immune_All_High模型（這些模型對(duì)免疫細(xì)胞類(lèi)型進(jìn)行更精細(xì)的注釋級(jí)別（低）和更粗略的注釋級(jí)別（高））：

model_low = models.Model.load(model="Immune_All_Low.pkl")
model_high = models.Model.load(model="Immune_All_High.pkl")

對(duì)于其中的每一個(gè)毅戈，我們都可以查看它包含哪些細(xì)胞類(lèi)型苹丸，以查看是否包含骨髓細(xì)胞類(lèi)型：

model_high.cell_types

輸出結(jié)果：

array(['B cells', 'B-cell lineage', 'Cycling cells', 'DC', 'DC precursor',
       'Double-negative thymocytes', 'Double-positive thymocytes', 'ETP',
       'Early MK', 'Endothelial cells', 'Epithelial cells',
       'Erythrocytes', 'Erythroid', 'Fibroblasts', 'Granulocytes',
       'HSC/MPP', 'ILC', 'ILC precursor', 'MNP', 'Macrophages',
       'Mast cells', 'Megakaryocyte precursor',
       'Megakaryocytes/platelets', 'Mono-mac', 'Monocyte precursor',
       'Monocytes', 'Myelocytes', 'Plasma cells', 'Promyelocytes',
       'T cells', 'pDC', 'pDC precursor'], dtype=object)

model_low.cell_types

輸出結(jié)果：

array(['Age-associated B cells', 'Alveolar macrophages', 'B cells',
       'CD16+ NK cells', 'CD16- NK cells', 'CD8a/a', 'CD8a/b(entry)',
       'CMP', 'CRTAM+ gamma-delta T cells', 'Classical monocytes',
       'Cycling B cells', 'Cycling DCs', 'Cycling NK cells',
       'Cycling T cells', 'Cycling gamma-delta T cells',
       'Cycling monocytes', 'DC', 'DC precursor', 'DC1', 'DC2', 'DC3',
       'Double-negative thymocytes', 'Double-positive thymocytes', 'ELP',
       'ETP', 'Early MK', 'Early erythroid', 'Early lymphoid/T lymphoid',
       'Endothelial cells', 'Epithelial cells', 'Erythrocytes',
       'Erythrophagocytic macrophages', 'Fibroblasts',
       'Follicular B cells', 'Follicular helper T cells', 'GMP',
       'Germinal center B cells', 'Granulocytes', 'HSC/MPP',
       'Hofbauer cells', 'ILC', 'ILC precursor', 'ILC1', 'ILC2', 'ILC3',
       'Intermediate macrophages', 'Intestinal macrophages',
       'Kidney-resident macrophages', 'Kupffer cells',
       'Large pre-B cells', 'Late erythroid', 'MAIT cells', 'MEMP', 'MNP',
       'Macrophages', 'Mast cells', 'Megakaryocyte precursor',
       'Megakaryocyte-erythroid-mast cell progenitor',
       'Megakaryocytes/platelets', 'Memory B cells',
       'Memory CD4+ cytotoxic T cells', 'Mid erythroid', 'Migratory DCs',
       'Mono-mac', 'Monocyte precursor', 'Monocytes', 'Myelocytes',
       'NK cells', 'NKT cells', 'Naive B cells',
       'Neutrophil-myeloid progenitor', 'Neutrophils',
       'Non-classical monocytes', 'Plasma cells', 'Plasmablasts',
       'Pre-pro-B cells', 'Pro-B cells',
       'Proliferative germinal center B cells', 'Promyelocytes',
       'Regulatory T cells', 'Small pre-B cells', 'T(agonist)',
       'Tcm/Naive cytotoxic T cells', 'Tcm/Naive helper T cells',
       'Tem/Effector helper T cells', 'Tem/Effector helper T cells PD1+',
       'Tem/Temra cytotoxic T cells', 'Tem/Trm cytotoxic T cells',
       'Transitional B cells', 'Transitional DC', 'Transitional NK',
       'Treg(diff)', 'Trm cytotoxic T cells', 'Type 1 helper T cells',
       'Type 17 helper T cells', 'gamma-delta T cells', 'pDC',
       'pDC precursor'], dtype=object)

看起來(lái)這些模型包括許多不同的免疫細(xì)胞類(lèi)型祖細(xì)胞。
現(xiàn)在讓我們運(yùn)行模型苇经。首先是更粗略的注釋級(jí)別：

predictions_high = celltypist.annotate(
    adata_celltypist, model=model_high, majority_voting=True
)

輸出結(jié)果：

?? Input data has 9370 cells and 31208 genes
?? Matching reference genes in the model
?? 6047 features used for prediction
?? Scaling input data
??? Predicting labels
? Prediction done!
?? Detected a neighborhood graph in the input object, will run over-clustering on the basis of it
?? Over-clustering input data with resolution set to 10
??? Majority voting the predictions
? Majority voting done!

將預(yù)測(cè)轉(zhuǎn)換為adata以獲得完整的輸出......

predictions_high_adata = predictions_high.to_adata()

然后將結(jié)果復(fù)制到原始的AnnData object：

adata.obs["celltypist_cell_label_coarse"] = predictions_high_adata.obs.loc[
    adata.obs.index, "majority_voting"
]
adata.obs["celltypist_conf_score_coarse"] = predictions_high_adata.obs.loc[
    adata.obs.index, "conf_score"
]

現(xiàn)在對(duì)于更精細(xì)的注釋級(jí)別也是如此：

predictions_low = celltypist.annotate(
    adata_celltypist, model=model_low, majority_voting=True
)

輸出結(jié)果：

?? Input data has 9370 cells and 31208 genes
?? Matching reference genes in the model
?? 6047 features used for prediction
?? Scaling input data
??? Predicting labels
? Prediction done!
?? Detected a neighborhood graph in the input object, will run over-clustering on the basis of it
?? Over-clustering input data with resolution set to 10
??? Majority voting the predictions
? Majority voting done!

predictions_low_adata = predictions_low.to_adata()

adata.obs["celltypist_cell_label_fine"] = predictions_low_adata.obs.loc[
    adata.obs.index, "majority_voting"
]
adata.obs["celltypist_conf_score_fine"] = predictions_low_adata.obs.loc[
    adata.obs.index, "conf_score"
]

現(xiàn)在畫(huà)圖：

sc.pl.umap(
    adata,
    color=["celltypist_cell_label_coarse", "celltypist_conf_score_coarse"],
    frameon=False,
    sort_order=False,
    wspace=1,
)

輸出結(jié)果：

... storing 'manual_celltype_annotation' as categorical

sc.pl.umap(
    adata,
    color=["celltypist_cell_label_fine", "celltypist_conf_score_fine"],
    frameon=False,
    sort_order=False,
    wspace=1,
)

了解這些注釋質(zhì)量的一種方法是查看觀察到的細(xì)胞類(lèi)型相似性是否符合我們的預(yù)期：

sc.pl.dendrogram(adata, groupby="celltypist_cell_label_fine")

輸出結(jié)果：

WARNING: dendrogram data not found (using key=dendrogram_celltypist_cell_label_fine). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.

該樹(shù)狀圖部分反映了關(guān)于細(xì)胞類(lèi)型關(guān)系的先驗(yàn)知識(shí)（例如B細(xì)胞主要聚集在一起）赘理，但我們也觀察到一些意想不到的模式：Tcm/Naive helper T細(xì)胞與紅細(xì)胞和巨噬細(xì)胞聚集，而不是與其他T細(xì)胞聚集扇单。這是一個(gè)危險(xiǎn)信號(hào)商模！可能Tcm/Naive helper T細(xì)胞注釋是錯(cuò)誤的。
現(xiàn)在讓我們看一下之前的手動(dòng)注釋：

sc.pl.umap(
    adata,
    color=["manual_celltype_annotation"],
    frameon=False,
)

可以看到幼稚B細(xì)胞注釋與自動(dòng)細(xì)胞注釋的幼稚B細(xì)胞一部分非常對(duì)應(yīng)蜘澜。同樣施流，我們所說(shuō)的過(guò)渡性B細(xì)胞的一部分在注釋中被稱為“small pre-B cells”，而我們的B1 B細(xì)胞對(duì)應(yīng)于它們的記憶B細(xì)胞鄙信。
然而瞪醋，您還會(huì)注意到，我們的HSC + MK/E prog簇在精細(xì)注釋中被注釋為T(mén)細(xì)胞和 HSC/多能祖細(xì)胞的混合物扮碧，因此這些注釋部分是矛盾的趟章。查看兩個(gè)注釋的置信度分?jǐn)?shù)，我們發(fā)現(xiàn)大部分細(xì)胞的注釋是以相對(duì)較低的置信度完成的慎王，這是一個(gè)有用的指示，表明如果沒(méi)有仔細(xì)驗(yàn)證和手動(dòng)審查宏侍，這些注釋就無(wú)法復(fù)制赖淤！
請(qǐng)參閱此處根據(jù)精細(xì)celltypist標(biāo)簽對(duì)簇12進(jìn)行的細(xì)分：

pd.crosstab(adata.obs.leiden_2, adata.obs.celltypist_cell_label_fine).loc[
    "12", :
].sort_values(ascending=False)

輸出結(jié)果：

celltypist_cell_label_fine
Naive B cells                  98
HSC/MPP                        97
Classical monocytes            56
Pro-B cells                    28
Tcm/Naive helper T cells       13
Tem/Temra cytotoxic T cells     1
Alveolar macrophages            0
CD16+ NK cells                  0
MAIT cells                      0
Memory B cells                  0
Mid erythroid                   0
NK cells                        0
Non-classical monocytes         0
Small pre-B cells               0
Tem/Trm cytotoxic T cells       0
Name: 12, dtype: int64

在較粗略的celltypist標(biāo)簽中，我們觀察到不同的模式：我們的簇12主要注釋為B細(xì)胞或巨核細(xì)胞前體谅河，同樣僅部分對(duì)應(yīng)于我們的注釋咱旱。

pd.crosstab(adata.obs.leiden_2, adata.obs.celltypist_cell_label_coarse).loc[
    "12", :
].sort_values(ascending=False)

輸出結(jié)果：

celltypist_cell_label_coarse
B cells                    98
Megakaryocyte precursor    97
ILC                        53
B-cell lineage             28
Erythroid                  13
Monocytes                   3
T cells                     1
DC                          0
Name: 12, dtype: int64

因此，我們看到這種自動(dòng)注釋僅部分對(duì)應(yīng)于我們的手動(dòng)注釋绷耍，甚至其本身的粗標(biāo)簽和細(xì)標(biāo)簽之間是矛盾的吐限。上面討論了此失敗的可能原因。
這強(qiáng)調(diào)了自動(dòng)注釋算法應(yīng)謹(jǐn)慎使用褂始，并應(yīng)被視為注釋數(shù)據(jù)的起點(diǎn)诸典，而不是最終注釋。 最終崎苗，已知標(biāo)記基因的表達(dá)仍然是細(xì)胞類(lèi)型注釋最被接受的支持狐粱。
為了強(qiáng)調(diào)這一點(diǎn)，讓我們看一下紅細(xì)胞譜系的標(biāo)記：血球蛋白B胆数。最有可能注釋為“Tcm/Naive helper T”的細(xì)胞（已根據(jù)上面的樹(shù)狀圖標(biāo)記為可能錯(cuò)誤注釋）來(lái)自紅細(xì)胞譜系肌蜻。

sc.pl.umap(adata, color="HBB", cmap="Reds", frameon=False, sort_order=False)

9.4.4 通過(guò)映射到參考進(jìn)行注釋

最后一種標(biāo)注數(shù)據(jù)的方式是基于將數(shù)據(jù)映射到已有的、已標(biāo)注的單細(xì)胞參考必尼，然后使用產(chǎn)生的聯(lián)合嵌入進(jìn)行標(biāo)簽遷移蒋搜。例如，此參考可以是您之前手動(dòng)注釋的單個(gè)樣本，之后您希望將這些注釋傳輸?shù)狡渌麛?shù)據(jù)集豆挽∷嵝荩或者，它也可以是已發(fā)布且最好精心策劃的現(xiàn)有參考資料祷杈。在這種情況下斑司，我們將“新數(shù)據(jù)”，即要映射和注釋的數(shù)據(jù)稱為“query”但汞。
有多種現(xiàn)有方法可以執(zhí)行此類(lèi)“query-to-reference映射”宿刮，包括scArches、Symphony和Azimuth (Seurat) 私蕾。 所有這些方法都使您能夠?qū)⑿聰?shù)據(jù)集映射到現(xiàn)有參考僵缺，而無(wú)需重新集成參考中的數(shù)據(jù)，也無(wú)需訪問(wèn)完整的參考數(shù)據(jù)踩叭。
我們將使用scArches作為基于參考映射的標(biāo)簽傳輸?shù)氖纠某保袁F(xiàn)有（variational autoencoder-based）模型為基礎(chǔ)，該模型將參考數(shù)據(jù)嵌入低維容贝、批量校正的空間中自脯。然后，它稍微擴(kuò)展該模型斤富，以將看不見(jiàn)的數(shù)據(jù)集映射到相同的“潛在空間”（即低維嵌入）膏潮。該模型擴(kuò)展還可以學(xué)習(xí)和消除映射數(shù)據(jù)集中存在的批次效應(yīng)。
現(xiàn)在满力，我們將展示如何使用scArches將數(shù)據(jù)映射到引用焕参，并使用此映射執(zhí)行從引用到新數(shù)據(jù)的標(biāo)簽傳輸（“query”）。
讓我們首先準(zhǔn)備數(shù)據(jù)以映射到參考油额。scArches是一種能夠使現(xiàn)有參考模型適應(yīng)新數(shù)據(jù)的方法叠纷，需要原始的、非標(biāo)準(zhǔn)化的計(jì)數(shù)潦嘶。因此涩嚣，我們將保留計(jì)數(shù)layer并從要映射的數(shù)據(jù)中刪除所有其他layers。我們也將.X設(shè)置為這些原始計(jì)數(shù)衬以。

adata_to_map = adata.copy()
for layer in list(adata_to_map.layers.keys()):
    if layer != "counts":
        del adata_to_map.layers[layer]
adata_to_map.X = adata_to_map.layers["counts"]

此外缓艳，重要的是我們使用與訓(xùn)練參考模型相同的輸入特征（即基因），并且將這些特征按相同的順序排列看峻。參考模型的特征信息與模型一起存儲(chǔ)阶淘。讓我們加載特征表。

reference_model_features = pd.read_csv(
    "https://figshare.com/ndownloader/files/41436645", index_col=0
)

該表包含基因名稱和基因ID互妓。由于基因ID通常不太容易因基因組注釋版本而發(fā)生變化溪窒，因此我們將使用它們來(lái)對(duì)我們的數(shù)據(jù)進(jìn)行子集化坤塞。因此，我們將數(shù)據(jù)和參考模型特征的行名稱設(shè)置為gene_ids澈蚌。重要的是摹芙，我們必須確保還存儲(chǔ)基因名稱以供以后使用：這些比基因ID更容易理解。

adata_to_map.var["gene_names"] = adata_to_map.var.index
adata_to_map.var.set_index("gene_ids", inplace=True)

reference_model_features["gene_names"] = reference_model_features.index
reference_model_features.set_index("gene_ids", inplace=True)

現(xiàn)在宛瞄，讓我們檢查一下查詢數(shù)據(jù)中是否有我們需要的所有基因：

print("Total number of genes needed for mapping:", reference_model_features.shape[0])

輸出結(jié)果：

Total number of genes needed for mapping: 4000

print(
    "Number of genes found in query dataset:",
    adata_to_map.var.index.isin(reference_model_features.index).sum(),
)

輸出結(jié)果：

Number of genes found in query dataset: 3998

我們?nèi)鄙僖恍┗蚋『獭Ｎ覀儗⑹謩?dòng)添加這些基因并將其計(jì)數(shù)設(shè)置為 0，因?yàn)槲覀兊臄?shù)據(jù)中似乎未檢測(cè)到這些基因份汗。讓我們?yōu)槟切┲挥辛阒档娜笔Щ騽?chuàng)建一個(gè)AnnData對(duì)象（包括我們的原始計(jì)數(shù)層盈电，它將用于映射）。之后我們會(huì)將其連接到我們自己的AnnData對(duì)象杯活。

missing_genes = [
    gene_id
    for gene_id in reference_model_features.index
    if gene_id not in adata_to_map.var.index
]

missing_gene_adata = sc.AnnData(
    X=csr_matrix(np.zeros(shape=(adata.n_obs, len(missing_genes))), dtype="float32"),
    obs=adata.obs.iloc[:, :1],
    var=reference_model_features.loc[missing_genes, :],
)
missing_gene_adata.layers["counts"] = missing_gene_adata.X

將我們的原始數(shù)據(jù)連接到缺失的基因數(shù)據(jù)匆帚。為了確保我們可以毫無(wú)錯(cuò)誤地執(zhí)行此連接，我們將從varm中刪除PCA矩陣旁钧。

if "PCs" in adata_to_map.varm.keys():
    del adata_to_map.varm["PCs"]

adata_to_map_augmented = sc.concat(
    [adata_to_map, missing_gene_adata],
    axis=1,
    join="outer",
    index_unique=None,
    merge="unique",
)

現(xiàn)在對(duì)模型中使用的基因進(jìn)行子集并正確排序：

adata_to_map_augmented = adata_to_map_augmented[
    :, reference_model_features.index
].copy()

檢查adata基因名稱是否與所需的基因順序完全對(duì)應(yīng)：

(adata_to_map_augmented.var.index == reference_model_features.index).all()

True

現(xiàn)在可以將基因索引設(shè)置回基因名稱以便于解釋：

adata_to_map_augmented.var["gene_ids"] = adata_to_map_augmented.var.index
adata_to_map_augmented.var.set_index("gene_names", inplace=True)

最后吸重，該參考模型使用adata.obs[‘batch’]作為我們的批處理變量。因此歪今，我們將檢查整個(gè)樣本是否已將其設(shè)置為一個(gè)值：

adata_to_map_augmented.obs.batch.unique()

輸出結(jié)果：

['12']
Categories (1, object): ['12']

現(xiàn)在我們來(lái)談?wù)勎覀兊膮⒖寄Ｐ秃啃摇Ｎ覀兊膮⒖寄Ｐ驮胶茫覀兊臉?biāo)簽傳輸性能就越好彤委。使用注釋良好的參考文獻(xiàn)集成了許多不同的數(shù)據(jù)集并且與您的數(shù)據(jù)很好地匹配（相同的器官鞭铆，相同的單細(xì)胞技術(shù)等）是理想的選擇：這樣的模型在各種數(shù)據(jù)集和批次上訓(xùn)練，因此預(yù)期對(duì)批次效應(yīng)更加穩(wěn)健焦影。然而，并非所有組織都存在這樣的參考封断。在本教程中斯辰，我們將使用骨髓樣本上訓(xùn)練的參考模型，不包括我們將要繪制的樣本坡疼。參考模型是一個(gè)scvi模型（用于數(shù)據(jù)集成）彬呻，它生成輸入數(shù)據(jù)的低維集成嵌入。
我們首先加載模型并向其傳遞我們想要映射的數(shù)據(jù)：

# loading model.pt from figshare
if not os.path.exists("./reference_model"):
    os.mkdir("./reference_model")
elif not os.path.exists("./reference_model/model.pt"):
    urllib.request.urlretrieve(
        "https://figshare.com/ndownloader/files/41436648",
        filename="reference_model/model.pt",
    )

scarches_model = sca.models.SCVI.load_query_data(
    adata=adata_to_map_augmented,
    reference_model="./reference_model",
    freeze_dropout=True,
)

輸出結(jié)果：

INFO     File ./reference_model/model.pt already downloaded   

Unable to initialize backend 'cuda': module 'jaxlib.xla_extension' has no attribute 'GpuAllocatorConfig'
Unable to initialize backend 'rocm': module 'jaxlib.xla_extension' has no attribute 'GpuAllocatorConfig'
Unable to initialize backend 'tpu': INVALID_ARGUMENT: TpuPlatform is not available.
Unable to initialize backend 'plugin': xla_extension has no attributes named get_plugin_device_client. Compile TensorFlow with //tensorflow/compiler/xla/python:enable_plugin_device set to true (defaults to false) to enable this.
No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

我們現(xiàn)在將更新此參考模型柄瑰，以便我們可以將我們自己的數(shù)據(jù)（“查詢”）嵌入到與參考相同的潛在空間中闸氮。這需要使用scArches對(duì)我們的查詢數(shù)據(jù)進(jìn)行訓(xùn)練：

scarches_model.train(max_epochs=500, plan_kwargs=dict(weight_decay=0.0))

輸出結(jié)果：

INFO: GPU available: True (cuda), used: True
GPU available: True (cuda), used: True
INFO: TPU available: False, using: 0 TPU cores
TPU available: False, using: 0 TPU cores
INFO: IPU available: False, using: 0 IPUs
IPU available: False, using: 0 IPUs
INFO: HPU available: False, using: 0 HPUs
HPU available: False, using: 0 HPUs
INFO: LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Epoch 500/500: 100%|██████████| 500/500 [03:55<00:00,  2.24it/s, v_num=1, train_loss_step=673, train_loss_epoch=653]

INFO: `Trainer.fit` stopped: `max_epochs=500` reached.
`Trainer.fit` stopped: `max_epochs=500` reached.

Epoch 500/500: 100%|██████████| 500/500 [03:55<00:00,  2.12it/s, v_num=1, train_loss_step=673, train_loss_epoch=653]

現(xiàn)在我們已經(jīng)更新了模型，我們可以計(jì)算query的（理想情況下批量校正的）潛在表示：

adata.obsm["X_scVI"] = scarches_model.get_latent_representation()

我們現(xiàn)在可以使用這個(gè)新計(jì)算的低維嵌入作為可視化和聚類(lèi)的基礎(chǔ)教沾。讓我們使用基于scVI的數(shù)據(jù)表示來(lái)計(jì)算新的UMAP蒲跨。

sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.umap(adata)

為了了解基于映射的UMAP是否具有普遍意義，讓我們看一下一些標(biāo)記以及它們的表達(dá)是否本地化到UMAP的特定部分：

sc.pl.umap(
    adata,
    color=["IGHD", "IGHM", "PRDM1"],
    vmin=0,
    vmax="p99",  # set vmax to the 99th percentile of the gene count instead of the maximum, to prevent outliers from making expression in other cells invisible. Note that this can cause problems for extremely lowly expressed genes.
    sort_order=False,  # do not plot highest expression on top, to not get a biased view of the mean expression among cells
    frameon=False,
    cmap="Reds",  # or choose another color map e.g. from here: https://matplotlib.org/stable/tutorials/colors/colormaps.html
)

現(xiàn)在重要的一步是我們可以將query數(shù)據(jù)的推斷潛在空間嵌入與現(xiàn)有的參考嵌入結(jié)合起來(lái)授翻。使用這種聯(lián)合嵌入或悲，我們不僅能夠可視化并將兩者聚類(lèi)在一起孙咪，但我們也可以進(jìn)行從query到參考的標(biāo)簽轉(zhuǎn)移。
讓我們加載參考嵌入：這通常與現(xiàn)有圖譜一起公開(kāi)提供巡语。

ref_emb = sc.read(
    filename="reference_embedding.h5ad",
    backup_url="https://figshare.com/ndownloader/files/41376264",
)

我們將存儲(chǔ)一個(gè)變量翎蹈，指定這些細(xì)胞來(lái)自引用。

ref_emb.obs["reference_or_query"] = "reference"

讓我們看看這個(gè)參照對(duì)象是什么：

ref_emb

輸出結(jié)果：

AnnData object with n_obs × n_vars = 86332 × 10
    obs: 'donor', 'batch', 'site', 'cell_type', 'reference_or_query'
    uns: 'neighbors', 'umap'
    obsm: 'X_umap'
    obsp: 'connectivities', 'distances'

可以看到男公，它只有10個(gè)維度( in .X )共同表示參考細(xì)胞的潛在空間嵌入荤堪。我們針對(duì)自己的數(shù)據(jù)計(jì)算的查詢嵌入也有10個(gè)維度。參考和查詢的10個(gè)維度是相同的枢赔，可以合并澄阳。
此外，它在.obs['cell_type']中具有細(xì)胞類(lèi)型標(biāo)簽糠爬。我們將使用這些標(biāo)簽來(lái)標(biāo)注我們自己的數(shù)據(jù)寇荧。
為了執(zhí)行標(biāo)簽遷移，我們首先將參考數(shù)據(jù)和查詢數(shù)據(jù)使用10維嵌入進(jìn)行拼接执隧。為了達(dá)到這個(gè)目的揩抡，我們將從我們的查詢數(shù)據(jù)中創(chuàng)建與我們從參考(與. X下的嵌入)中獲得的相同類(lèi)型的AnnData對(duì)象，并將兩者進(jìn)行連接镀琉。有了這一點(diǎn)峦嗤，我們可以聯(lián)合分析引用和查詢，包括進(jìn)行從一個(gè)到另一個(gè)的轉(zhuǎn)換屋摔。

adata_emb = sc.AnnData(X=adata.obsm["X_scVI"], obs=adata.obs)

adata_emb.obs["reference_or_query"] = "query"

emb_ref_query = sc.concat(
    [ref_emb, adata_emb],
    axis=0,
    join="outer",
    index_unique=None,
    merge="unique",
)

讓我們用UMAP可視化聯(lián)合嵌入烁设。

sc.pp.neighbors(emb_ref_query)
sc.tl.umap(emb_ref_query)

基于UMAP，我們可以直觀地獲得參考和query是否集成良好：

sc.pl.umap(
    emb_ref_query,
    color=["reference_or_query"],
    sort_order=False,
    frameon=False,
)

... storing 'reference_or_query' as categorical

在這個(gè)UMAP中钓试，query和引用的(部分)混合是一個(gè)很好的信號(hào)装黑，當(dāng)映射完全失敗時(shí)，在UMAP中經(jīng)常會(huì)看到query和引用的完全分離弓熏。
現(xiàn)在讓我們從參考文獻(xiàn)中看一下細(xì)胞類(lèi)型的注釋恋谭。所有來(lái)自查詢的細(xì)胞都被設(shè)置為NA，因?yàn)樗鼈冞€沒(méi)有注釋挽鞠，并以黑色顯示疚颊。
我們將把這個(gè)數(shù)字做得更大一點(diǎn)，這樣我們就可以很好地讀懂這個(gè)圖例：

sc.set_figure_params(figsize=(8, 8))

sc.pl.umap(
    emb_ref_query,
    color=["cell_type"],
    sort_order=False,
    frameon=False,
    legend_loc="on data",
    legend_fontsize=10,
    na_color="black",
)

如UMAP所示信认，我們可以通過(guò)查看周?chē)鷧⒖嫉募?xì)胞類(lèi)型來(lái)猜測(cè)我們每個(gè)細(xì)胞的細(xì)胞類(lèi)型(黑色)材义。這正是基于最近鄰圖的標(biāo)簽傳遞方法所做的：對(duì)于每個(gè)查詢細(xì)胞，它檢查其相鄰參考細(xì)胞中最常見(jiàn)的細(xì)胞類(lèi)型嫁赏。參考細(xì)胞來(lái)自單個(gè)細(xì)胞類(lèi)型的比例越高其掂，對(duì)標(biāo)簽傳遞越有信心。
讓我們執(zhí)行基于knn的標(biāo)簽傳遞橄教。
首先我們建立了標(biāo)簽遷移模型：

knn_transformer = sca.utils.knn.weighted_knn_trainer(
    train_adata=ref_emb,
    train_adata_emb="X",  # location of our joint embedding
    n_neighbors=15,
)

Weighted KNN with n_neighbors = 15 ... 

現(xiàn)在我們執(zhí)行標(biāo)簽轉(zhuǎn)移：

labels, uncert = sca.utils.knn.weighted_knn_transfer(
    query_adata=adata_emb,
    query_adata_emb="X",  # location of our embedding, query_adata.X in this case
    label_keys="cell_type",  # (start of) obs column name(s) for which to transfer labels
    knn_model=knn_transformer,
    ref_adata_obs=ref_emb.obs,
)

finished!

然后存儲(chǔ)結(jié)果至adata中：

adata_emb.obs["transf_cell_type"] = labels.loc[adata_emb.obs.index, "cell_type"]
adata_emb.obs["transf_cell_type_unc"] = uncert.loc[adata_emb.obs.index, "cell_type"]

讓我們將結(jié)果傳遞給我們查詢的一個(gè)數(shù)據(jù)對(duì)象，這個(gè)數(shù)據(jù)對(duì)象也有我們的UMAP和基因計(jì)數(shù)，這樣我們就可以將所有這些可視化在一起椒功。

adata.obs.loc[adata_emb.obs.index, "transf_cell_type"] = adata_emb.obs[
    "transf_cell_type"
]
adata.obs.loc[adata_emb.obs.index, "transf_cell_type_unc"] = adata_emb.obs[
    "transf_cell_type_unc"
]

我們現(xiàn)在可以在我們先前計(jì)算的自己數(shù)據(jù)的UMAP中可視化轉(zhuǎn)移的標(biāo)簽：
讓我們?cè)侔褕D形的大小設(shè)置得更小一些：

sc.set_figure_params(figsize=(5, 5))

sc.pl.umap(adata, color="transf_cell_type", frameon=False)

... storing 'transf_cell_type' as categorical

基于每個(gè)query細(xì)胞的鄰居，我們不僅可以猜測(cè)這些細(xì)胞所屬的細(xì)胞類(lèi)型翩迈，而且還可以生成該標(biāo)簽的確定性度量：如果一個(gè)細(xì)胞有來(lái)自幾個(gè)不同細(xì)胞類(lèi)型的鄰居，我們的猜測(cè)將是高度不確定的盔夜。這與評(píng)估我們?cè)诙啻蟪潭壬峡梢?信任"轉(zhuǎn)移的標(biāo)簽有關(guān)负饲，我們將不確定度得分可視化：

sc.pl.umap(adata, color="transf_cell_type_unc", frameon=False)

讓我們檢查每個(gè)細(xì)胞類(lèi)型標(biāo)簽的標(biāo)簽傳遞不確定度水平有多高。這給我們的第一印象是哪些注釋更有爭(zhēng)議性/需要更多的人工檢查喂链。

fig, ax = plt.subplots(figsize=(8, 3))
ct_order = (
    adata.obs.groupby("transf_cell_type")
    .agg({"transf_cell_type_unc": "median"})
    .sort_values(by="transf_cell_type_unc", ascending=False)
)
sns.boxplot(
    adata.obs,
    x="transf_cell_type",
    y="transf_cell_type_unc",
    color="grey",
    ax=ax,
    order=ct_order.index,
)
ax.tick_params(rotation=90, axis="x")

你會(huì)注意到返十，例如祖細(xì)胞往往比其他細(xì)胞類(lèi)型更難區(qū)分。在我們的注釋中椭微，"Other T 細(xì)胞"這一非特異的分類(lèi)也是如此洞坑。我們看到了pDCs，這是一種已知的轉(zhuǎn)錄截然不同的細(xì)胞類(lèi)型蝇率，因此更容易識(shí)別和標(biāo)記迟杂。
為了將這種不確定性信息包含在我們轉(zhuǎn)移的標(biāo)簽中，我們可以將不確定性得分高于0.2的細(xì)胞設(shè)置為"未知"：

adata.obs["transf_cell_type_certain"] = adata.obs.transf_cell_type.tolist()
adata.obs.loc[
    adata.obs.transf_cell_type_unc > 0.2, "transf_cell_type_certain"
] = "Unknown"

讓我們看看我們的注釋在經(jīng)過(guò)這種過(guò)濾之后的樣子本慕。注意圖例中的未知顏色和UMAP排拷。

sc.pl.umap(adata, color="transf_cell_type_certain", frameon=False)

... storing 'transf_cell_type_certain' as categorical

為了便于識(shí)別，我們只對(duì)"未知"細(xì)胞進(jìn)行著色锅尘。這將使我們更容易看到其中有多少监氢。你可以用其他任何細(xì)胞類(lèi)型的標(biāo)簽做同樣的事情。

sc.pl.umap(adata, color="transf_cell_type_certain", groups="Unknown")

這些細(xì)胞有很多藤违，將需要特別仔細(xì)的人工審查浪腐。然而，圍繞"未知細(xì)胞"的低不確定性注釋將使我們首先了解我們可以期望每個(gè)細(xì)胞屬于什么細(xì)胞類(lèi)型顿乒。
現(xiàn)在讓我們來(lái)看看我們更為確定的注釋牛欢。我們將檢查幾個(gè)細(xì)胞類(lèi)型(在此隨機(jī)選擇)在多大程度上與上述已知的標(biāo)記基因相匹配。在現(xiàn)實(shí)中淆游，這應(yīng)該對(duì)所有的注釋進(jìn)行系統(tǒng)的整理。

cell_types_to_check = [
    "CD14+ Mono",
    "cDC2",
    "NK",
    "B1 B",
    "CD4+ T activated",
    "T naive",
    "MK/E prog",
]

對(duì)于這些細(xì)胞類(lèi)型隔盛，我們?cè)谧值渲卸加袠?biāo)記犹菱。讓我們對(duì)所有新注釋的細(xì)胞類(lèi)型進(jìn)行plot標(biāo)記表達(dá)。你會(huì)注意到吮炕，標(biāo)記的表達(dá)一般對(duì)應(yīng)著自動(dòng)化的注釋腊脱。

sc.pl.dotplot(
    adata,
    var_names={
        ct: marker_genes_in_data[ct] for ct in cell_types_to_check
    },  # gene names grouped by cell type in a dictionary
    groupby="transf_cell_type_certain",
    standard_scale="var",  # normalize gene scores from 0 to 1
)

正如您所看到的，標(biāo)記組通常在用匹配標(biāo)簽注釋的細(xì)胞中表達(dá)最高龙亲。這意味著這些標(biāo)簽可能（至少部分）正確陕凹！
讓我們?cè)俅位氐匠錆M不確定性的UMAP：

sc.pl.umap(
    adata, color=["transf_cell_type_unc", "transf_cell_type_certain"], frameon=False
)

這種不確定性不僅可以幫助我們識(shí)別算法不確定細(xì)胞屬于哪種細(xì)胞類(lèi)型的區(qū)域（例如悍抑，因?yàn)樗挥趦蓚€(gè)注釋的表型之間），而且還可以突出顯示看不見(jiàn)的細(xì)胞類(lèi)型或新的細(xì)胞狀態(tài)杜耙。例如搜骡，您的參考可能由健康細(xì)胞組成，而您的查詢可能來(lái)自患病樣本佑女。然后记靡，不確定性得分可以突出顯示特定于疾病的細(xì)胞狀態(tài)，因?yàn)樗鼈兛赡軟](méi)有來(lái)自始終來(lái)自單一細(xì)胞類(lèi)型的參考的鄰居团驱。特別是當(dāng)您的參考基于大量數(shù)據(jù)集時(shí)摸吠，不確定性分?jǐn)?shù)對(duì)于標(biāo)記可能值得研究的查詢數(shù)據(jù)部分非常有用。因此嚎花，基于參考的標(biāo)簽傳輸不僅可以幫助您注釋數(shù)據(jù)寸痢，還可以加快數(shù)據(jù)的探索和解釋。然而紊选，與任何指標(biāo)一樣啼止，這些不確定性分?jǐn)?shù)通常并不完美，在某些情況下無(wú)法突出新的細(xì)胞類(lèi)型或狀態(tài)丛楚。
與本篇中討論的任何方法一樣族壳，傳輸注釋的質(zhì)量取決于“訓(xùn)練數(shù)據(jù)”（在本例中為參考）及其注釋的質(zhì)量、模型的質(zhì)量以及您自己的匹配數(shù)據(jù)與訓(xùn)練數(shù)據(jù)趣些！
因此仿荆，應(yīng)始終通過(guò)使用標(biāo)記基因表達(dá)的手動(dòng)檢查來(lái)驗(yàn)證轉(zhuǎn)移注釋的質(zhì)量，并且可能需要對(duì)初始注釋進(jìn)行細(xì)化坏平。

最后拢操，如果需要，存儲(chǔ)adata對(duì)象：

# adata.write("./annotation_adata_out.h5ad")