9.4 自动注释

9.4.1 总论

这一篇讨论的方法将是自动化的方法，而不是手动注释数据。与上一篇展示的方法不同，这些方法中的每一种都可以使您以自动化的方式对数据进行注释。它们基于不同的原则，有时需要预先定义的标记集合，或者在预先存在的完整scRNA-seq数据集上训练。
如前所述，自动生成的注释的质量可能会有所不同。更具体地说，注释的质量取决于：

• 选择的分类器类型：之前的基准研究表明，不同类型的分类器通常表现相当，基于神经网络的方法通常不会优于通用模型如支持向量机或线性回归模型。
• 训练分类器的数据的质量。如果训练数据没有很好地注释或以低分辨率注释，分类器也会做同样的事情。同样，如果训练数据和/或其注释有噪声，分类器可能表现不佳。
• 自己的数据与训练分类器的数据的相似度。例如，如果分类器是在drop-seq单细胞数据集上进行训练的，并且数据是10X单核而不是单细胞drop-seq，则这可能会降低注释的质量。在包括多种数据集的跨数据集图集上训练的分类器可能比在单个数据集上训练的分类器提供更稳健和更好质量的注释。

上述几点强调了使用分类器可能存在的缺点，具体取决于训练数据和模型类型。尽管如此，使用预先训练的分类器来注释数据有几个重要的优点。首先，它是一种快速、简单的数据注释方法。注释不需要下载或预处理训练数据，有时只需将数据上传到在线网页。其次，这些方法不像手动注释那样依赖于将数据划分为集群。第三，预先训练的分类器使您能够直接利用以前研究中的知识和信息，例如高质量的注释。
最后，由于这些分类器的透明性往往不如基于人工标记的标注，一个好的不确定性度量量化标注将提高方法的质量和可用性。我们将在下面更广泛地讨论这个问题。

9.4.2 基于marker基因的分类器

一类自动化细胞类型注释方法依赖于一组预定义的标记基因。根据这些标记基因的表达水平，细胞被分为细胞类型。这些模型所基于的标记基因集越稳健、越通用，模型的性能就越好。然而，与其他模型一样，它们可能会受到模型训练数据和需要标记的数据之间与批次效应相关的差异的影响。与基于较大基因集的模型（见下文）相比，这些方法的优点之一是它们更加透明：我们知道根据哪些基因进行分类。

9.4.3 基于更广泛的基因集的分类器

值得注意的是，到目前为止讨论的方法仅使用数据中检测到的基因的一小部分：通常每种细胞类型仅使用一组1至10个标记基因。另一种方法是使用一个分类器，该分类器将更大的基因集（数千个或更多）作为输入，从而更多地利用scRNA-seq数据的广度。此类分类器是在先前注释的数据集或图谱上进行训练的。其中的示例包括 CellTypist（另请参阅 https://www.celltypist.org，其中数据可以上传到门户以获取自动细胞注释）和 Clustifyr。
让我们在数据上尝试使用CellTypist。根据CellTypist教程(https://www.celltypist.org/tutorials)，我们需要准备数据，以便将计数标准化为每个细胞10,000个计数，然后进行log1p转换：

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()

我们现在将下载免疫细胞的celltypist模型：

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

输出结果：

📜 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模型（这些模型对免疫细胞类型进行更精细的注释级别（低）和更粗略的注释级别（高））：

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

对于其中的每一个，我们都可以查看它包含哪些细胞类型，以查看是否包含骨髓细胞类型：

model_high.cell_types

输出结果：

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

输出结果：

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)

看起来这些模型包括许多不同的免疫细胞类型祖细胞。
现在让我们运行模型。首先是更粗略的注释级别：

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

输出结果：

🔬 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!

将预测转换为adata以获得完整的输出......

predictions_high_adata = predictions_high.to_adata()

然后将结果复制到原始的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"
]

现在对于更精细的注释级别也是如此：

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

输出结果：

🔬 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"
]

现在画图：

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

输出结果：

... 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,
)

了解这些注释质量的一种方法是查看观察到的细胞类型相似性是否符合我们的预期：

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

输出结果：

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.

该树状图部分反映了关于细胞类型关系的先验知识（例如B细胞主要聚集在一起），但我们也观察到一些意想不到的模式：Tcm/Naive helper T细胞与红细胞和巨噬细胞聚集，而不是与其他T细胞聚集。这是一个危险信号！可能Tcm/Naive helper T细胞注释是错误的。
现在让我们看一下之前的手动注释：

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

可以看到幼稚B细胞注释与自动细胞注释的幼稚B细胞一部分非常对应。同样，我们所说的过渡性B细胞的一部分在注释中被称为“small pre-B cells”，而我们的B1 B细胞对应于它们的记忆B细胞。
然而，您还会注意到，我们的HSC + MK/E prog簇在精细注释中被注释为T细胞和 HSC/多能祖细胞的混合物，因此这些注释部分是矛盾的。查看两个注释的置信度分数，我们发现大部分细胞的注释是以相对较低的置信度完成的，这是一个有用的指示，表明如果没有仔细验证和手动审查，这些注释就无法复制！
请参阅此处根据精细celltypist标签对簇12进行的细分：

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

输出结果：

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标签中，我们观察到不同的模式：我们的簇12主要注释为B细胞或巨核细胞前体，同样仅部分对应于我们的注释。

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

输出结果：

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

因此，我们看到这种自动注释仅部分对应于我们的手动注释，甚至其本身的粗标签和细标签之间是矛盾的。上面讨论了此失败的可能原因。
这强调了自动注释算法应谨慎使用，并应被视为注释数据的起点，而不是最终注释。 最终，已知标记基因的表达仍然是细胞类型注释最被接受的支持。
为了强调这一点，让我们看一下红细胞谱系的标记：血球蛋白B。最有可能注释为“Tcm/Naive helper T”的细胞（已根据上面的树状图标记为可能错误注释）来自红细胞谱系。

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

9.4.4 通过映射到参考进行注释

最后一种标注数据的方式是基于将数据映射到已有的、已标注的单细胞参考，然后使用产生的联合嵌入进行标签迁移。例如，此参考可以是您之前手动注释的单个样本，之后您希望将这些注释传输到其他数据集。或者，它也可以是已发布且最好精心策划的现有参考资料。在这种情况下，我们将“新数据”，即要映射和注释的数据称为“query”。
有多种现有方法可以执行此类“query-to-reference映射”，包括scArches、Symphony和Azimuth (Seurat) 。 所有这些方法都使您能够将新数据集映射到现有参考，而无需重新集成参考中的数据，也无需访问完整的参考数据。
我们将使用scArches作为基于参考映射的标签传输的示例，它以现有（variational autoencoder-based）模型为基础，该模型将参考数据嵌入低维、批量校正的空间中。然后，它稍微扩展该模型，以将看不见的数据集映射到相同的“潜在空间”（即低维嵌入）。该模型扩展还可以学习和消除映射数据集中存在的批次效应。
现在，我们将展示如何使用scArches将数据映射到引用，并使用此映射执行从引用到新数据的标签传输（“query”）。
让我们首先准备数据以映射到参考。scArches是一种能够使现有参考模型适应新数据的方法，需要原始的、非标准化的计数。因此，我们将保留计数layer并从要映射的数据中删除所有其他layers。我们也将.X设置为这些原始计数。

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"]

此外，重要的是我们使用与训练参考模型相同的输入特征（即基因），并且将这些特征按相同的顺序排列。参考模型的特征信息与模型一起存储。让我们加载特征表。

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

该表包含基因名称和基因ID。由于基因ID通常不太容易因基因组注释版本而发生变化，因此我们将使用它们来对我们的数据进行子集化。因此，我们将数据和参考模型特征的行名称设置为gene_ids。重要的是，我们必须确保还存储基因名称以供以后使用：这些比基因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)

现在，让我们检查一下查询数据中是否有我们需要的所有基因：

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

输出结果：

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(),
)

输出结果：

Number of genes found in query dataset: 3998

我们缺少一些基因。我们将手动添加这些基因并将其计数设置为 0，因为我们的数据中似乎未检测到这些基因。让我们为那些只有零值的缺失基因创建一个AnnData对象（包括我们的原始计数层，它将用于映射）。之后我们会将其连接到我们自己的AnnData对象。

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

将我们的原始数据连接到缺失的基因数据。为了确保我们可以毫无错误地执行此连接，我们将从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",
)

现在对模型中使用的基因进行子集并正确排序：

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

检查adata基因名称是否与所需的基因顺序完全对应：

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

True

现在可以将基因索引设置回基因名称以便于解释：

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’]作为我们的批处理变量。因此，我们将检查整个样本是否已将其设置为一个值：

adata_to_map_augmented.obs.batch.unique()

输出结果：

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

现在我们来谈谈我们的参考模型。我们的参考模型越好，我们的标签传输性能就越好。使用注释良好的参考文献集成了许多不同的数据集并且与您的数据很好地匹配（相同的器官，相同的单细胞技术等）是理想的选择：这样的模型在各种数据集和批次上训练，因此预期对批次效应更加稳健。然而，并非所有组织都存在这样的参考。在本教程中，我们将使用骨髓样本上训练的参考模型，不包括我们将要绘制的样本。参考模型是一个scvi模型（用于数据集成），它生成输入数据的低维集成嵌入。
我们首先加载模型并向其传递我们想要映射的数据：

# 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,
)

输出结果：

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.)

我们现在将更新此参考模型，以便我们可以将我们自己的数据（“查询”）嵌入到与参考相同的潜在空间中。这需要使用scArches对我们的查询数据进行训练：

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

输出结果：

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]

现在我们已经更新了模型，我们可以计算query的（理想情况下批量校正的）潜在表示：

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

我们现在可以使用这个新计算的低维嵌入作为可视化和聚类的基础。让我们使用基于scVI的数据表示来计算新的UMAP。

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

为了了解基于映射的UMAP是否具有普遍意义，让我们看一下一些标记以及它们的表达是否本地化到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
)

现在重要的一步是我们可以将query数据的推断潜在空间嵌入与现有的参考嵌入结合起来。使用这种联合嵌入，我们不仅能够可视化并将两者聚类在一起，但我们也可以进行从query到参考的标签转移。
让我们加载参考嵌入：这通常与现有图谱一起公开提供。

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

我们将存储一个变量，指定这些细胞来自引用。

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

让我们看看这个参照对象是什么：

ref_emb

输出结果：

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个维度( in .X )共同表示参考细胞的潜在空间嵌入。我们针对自己的数据计算的查询嵌入也有10个维度。参考和查询的10个维度是相同的，可以合并。
此外，它在.obs['cell_type']中具有细胞类型标签。我们将使用这些标签来标注我们自己的数据。
为了执行标签迁移，我们首先将参考数据和查询数据使用10维嵌入进行拼接。为了达到这个目的，我们将从我们的查询数据中创建与我们从参考(与. X下的嵌入)中获得的相同类型的AnnData对象，并将两者进行连接。有了这一点，我们可以联合分析引用和查询，包括进行从一个到另一个的转换。

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可视化联合嵌入。

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

在这个UMAP中，query和引用的(部分)混合是一个很好的信号，当映射完全失败时，在UMAP中经常会看到query和引用的完全分离。
现在让我们从参考文献中看一下细胞类型的注释。所有来自查询的细胞都被设置为NA，因为它们还没有注释，并以黑色显示。
我们将把这个数字做得更大一点，这样我们就可以很好地读懂这个图例：

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所示，我们可以通过查看周围参考的细胞类型来猜测我们每个细胞的细胞类型(黑色)。这正是基于最近邻图的标签传递方法所做的：对于每个查询细胞，它检查其相邻参考细胞中最常见的细胞类型。参考细胞来自单个细胞类型的比例越高，对标签传递越有信心。
让我们执行基于knn的标签传递。
首先我们建立了标签迁移模型：

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 ... 

现在我们执行标签转移：

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!

然后存储结果至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"]

让我们将结果传递给我们查询的一个数据对象，这个数据对象也有我们的UMAP和基因计数，这样我们就可以将所有这些可视化在一起。

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"
]

我们现在可以在我们先前计算的自己数据的UMAP中可视化转移的标签：
让我们再把图形的大小设置得更小一些：

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

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

... storing 'transf_cell_type' as categorical

基于每个query细胞的邻居，我们不仅可以猜测这些细胞所属的细胞类型，而且还可以生成该标签的确定性度量：如果一个细胞有来自几个不同细胞类型的邻居，我们的猜测将是高度不确定的。这与评估我们在多大程度上可以"信任"转移的标签有关，我们将不确定度得分可视化：

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

让我们检查每个细胞类型标签的标签传递不确定度水平有多高。这给我们的第一印象是哪些注释更有争议性/需要更多的人工检查。

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")

你会注意到，例如祖细胞往往比其他细胞类型更难区分。在我们的注释中，"Other T 细胞"这一非特异的分类也是如此。我们看到了pDCs，这是一种已知的转录截然不同的细胞类型，因此更容易识别和标记。
为了将这种不确定性信息包含在我们转移的标签中，我们可以将不确定性得分高于0.2的细胞设置为"未知"：

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"

让我们看看我们的注释在经过这种过滤之后的样子。注意图例中的未知颜色和UMAP。

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

... storing 'transf_cell_type_certain' as categorical

为了便于识别，我们只对"未知"细胞进行着色。这将使我们更容易看到其中有多少。你可以用其他任何细胞类型的标签做同样的事情。

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

这些细胞有很多，将需要特别仔细的人工审查。然而，围绕"未知细胞"的低不确定性注释将使我们首先了解我们可以期望每个细胞属于什么细胞类型。
现在让我们来看看我们更为确定的注释。我们将检查几个细胞类型(在此随机选择)在多大程度上与上述已知的标记基因相匹配。在现实中，这应该对所有的注释进行系统的整理。

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

对于这些细胞类型，我们在字典中都有标记。让我们对所有新注释的细胞类型进行plot标记表达。你会注意到，标记的表达一般对应着自动化的注释。

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
)

正如您所看到的，标记组通常在用匹配标签注释的细胞中表达最高。这意味着这些标签可能（至少部分）正确！
让我们再次回到充满不确定性的UMAP：

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

这种不确定性不仅可以帮助我们识别算法不确定细胞属于哪种细胞类型的区域（例如，因为它位于两个注释的表型之间），而且还可以突出显示看不见的细胞类型或新的细胞状态。例如，您的参考可能由健康细胞组成，而您的查询可能来自患病样本。然后，不确定性得分可以突出显示特定于疾病的细胞状态，因为它们可能没有来自始终来自单一细胞类型的参考的邻居。特别是当您的参考基于大量数据集时，不确定性分数对于标记可能值得研究的查询数据部分非常有用。因此，基于参考的标签传输不仅可以帮助您注释数据，还可以加快数据的探索和解释。然而，与任何指标一样，这些不确定性分数通常并不完美，在某些情况下无法突出新的细胞类型或状态。
与本篇中讨论的任何方法一样，传输注释的质量取决于“训练数据”（在本例中为参考）及其注释的质量、模型的质量以及您自己的匹配数据与训练数据！
因此，应始终通过使用标记基因表达的手动检查来验证转移注释的质量，并且可能需要对初始注释进行细化。

最后，如果需要，存储adata对象：

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