3.2 Annotation¶

Motivation¶

为了更好地理解您的数据并利用现有知识，了解数据中每个细胞的“细胞身份”是很重要的。根据已知（或有时未知）的细胞表型对数据中的细胞群体进行标注的过程称为“细胞注释”。尽管有许多方法可以注释细胞（例如，基于批次、疾病、性别等），在本笔记本中我们将重点放在“细胞类型”的注释上。

那么，什么是细胞类型？生物学家使用“细胞类型”这个术语来表示一种在数据集之间具有鲁棒性、可以根据特定标记（即蛋白质或基因转录物）的表达识别的细胞表型，且通常与特定功能相关。例如，浆细胞B是一种分泌抗体以对抗病原体的白细胞类型，可以通过特定标记来识别。了解样本中有哪些细胞类型对于理解您的数据至关重要。例如，知道肿瘤中存在特定免疫细胞类型或骨髓样本中存在异常的造血干细胞可能会为您研究的疾病提供宝贵的见解。

然而，像任何分类一样，类别的大小和它们之间划定的边界部分是主观的，并且可能随时间变化。例如，由于新技术允许更高分辨率地观察细胞，或者发现以前被认为没有生物学意义的特定“亚表型”具有重要的生物学意义（参见例如{cite}anno:KadurLakshminarasimhaMurthy2022）。因此，细胞类型通常进一步分类为“亚型”或“细胞状态”（例如，活化状态与静止状态），一些研究人员使用“细胞身份”一词来避免细胞类型、细胞亚型和细胞状态之间有时具有随意性的区别。关于这一主题的详细讨论，我们推荐Wagner等人的综述{cite}anno:Wagner2016以及最近发布的Zeng的综述{cite}anno:ZENG20222739。

类似地，多个细胞类型可以是单一连续体的一部分，其中一个细胞类型可能过渡或分化为另一个。例如，在造血过程中，细胞从干细胞分化为特定的免疫细胞类型。尽管通常在这一分化的早期和晚期之间划定了硬边界，但这些细胞的状态更准确地可以描述为介于较少分化和更多分化的细胞表型之间的分化坐标。我们将在后续章节讨论分化和细胞轨迹。

那么，我们如何注释单细胞数据中的细胞呢？有多种方法可以做到这一点，我们将在下面概述不同的方法。由于我们处理的是转录组数据，每种方法最终都基于特定基因或基因集的表达，或细胞之间的一般转录组相似性。

Environment setup¶

我们将过滤掉一些不会影响代码的弃用和性能警告：

In [2]:

import warnings

warnings.filterwarnings("ignore", category=DeprecationWarning)

import numba
from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning

warnings.simplefilter("ignore", category=NumbaDeprecationWarning)

import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) import numba from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning warnings.simplefilter("ignore", category=NumbaDeprecationWarning)

加载所需的模块：

In [3]:

import scanpy as sc
import pandas as pd
import numpy as np
import os
from scipy.sparse import csr_matrix
import seaborn as sns
import matplotlib.pyplot as plt
import celltypist
from celltypist import models
import scarches as sca
import urllib.request

import scanpy as sc import pandas as pd import numpy as np import os from scipy.sparse import csr_matrix import seaborn as sns import matplotlib.pyplot as plt import celltypist from celltypist import models import scarches as sca import urllib.request

/home/icb/lisa.sikkema/miniconda3/envs/best_practices_annotation/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
In order to use the mouse gastrulation seqFISH datsets, please install squidpy (see https://github.com/scverse/squidpy).
Created a temporary directory at /tmp/tmpihngzax_
Writing /tmp/tmpihngzax_/_remote_module_non_scriptable.py
In order to use sagenet models, please install pytorch geometric (see https://pytorch-geometric.readthedocs.io) and 
 captum (see https://github.com/pytorch/captum).
mvTCR is not installed. To use mvTCR models, please install it first using "pip install mvtcr"
multigrate is not installed. To use multigrate models, please install it first using "pip install multigrate".

再过滤一个 pandas 警告：

In [4]:

warnings.filterwarnings("ignore", category=pd.errors.PerformanceWarning)

warnings.filterwarnings("ignore", category=pd.errors.PerformanceWarning)

我们将继续处理之前预处理过的scRNA-seq数据集，现在将对其进行注释。

Set figure parameters:

In [5]:

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

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

Load data¶

让我们读取我们将在本教程中使用的玩具数据集。该数据集包括一本书其他部分中使用的同一数据的单个样本 ("site4-donor8")。此外，未通过质控的细胞已被移除。

In [6]:

adata = sc.read(
    filename="s4d8_clustered.h5ad",
    backup_url="https://figshare.com/ndownloader/files/41436666",
)

adata = sc.read( filename="s4d8_clustered.h5ad", backup_url="https://figshare.com/ndownloader/files/41436666", )

Manual annotation¶

进行细胞类型注释的经典或最古老的方法是基于已知与特定细胞类型相关的单个或少量标记基因。这种方法可以追溯到“scRNA-seq 时代之前”，那时的单细胞数据是低维的（例如包含不超过30-40个基因面板的FACS数据）。这是一种快速且透明的方法来注释您的数据。然而，当某一特定细胞类型不存在独特的标记基因（这种情况常见）时，这种方法可能会变得更加复杂，甚至更少客观，需要使用标记基因的组合或表达阈值进行适当的注释。一套稳健的标记基因集和先前的知识或注释经验可以在这里提供帮助，但这种方法存在决策不清晰和主观的风险。

在这种情况下，数据通常在注释之前进行聚类，以便我们可以对细胞群体进行注释，而不是逐个细胞进行判断。这不仅减少了工作量，还能更有效地应对噪音问题：即使特定标记在某个细胞中表达，单个细胞也可能没有相应的计数，这只是由于单细胞数据固有的稀疏性。聚类使得可以检测到在总体基因表达上高度相似的细胞，从而可以解决单个细胞水平的掉落问题。

最后，可以从两个角度来进行基于标记基因的注释。一种选择是从一个包含所有预期细胞类型的标记基因的表格入手，检查这些标记基因在哪些聚类中表达。另一种选择是检查您定义的聚类中哪些基因表达水平较高，然后检查这些基因是否与已知的细胞类型或状态相关。如果需要，可以在这两种方法之间来回切换。我们将在下面展示这两种方法的示例。

From markers to cluster annotation¶

让我们从已知标记基因的方法开始。首先，我们将在此列出一组基于文献的骨髓细胞类型标记。这些标记基因来自于先前研究特定细胞类型和亚型的论文，并报告了这些细胞类型的标记基因。请注意，蛋白质水平的标记（例如用于FACS的标记）有时在转录组数据中效果不佳，因此使用基于RNA的论文中的标记往往更有效。此外，有时在一个数据集中表现良好的标记在其他数据集中可能效果不佳。因此，理想的标记集应该在多个数据集中进行验证。最后，与专家合作通常是有益的：作为一名生物信息学家，尝试与对组织、生物学、预期的细胞类型和标记等有更广泛知识的生物学家合作。

In [7]:

marker_genes = {
    "CD14+ Mono": ["FCN1", "CD14"],
    "CD16+ Mono": ["TCF7L2", "FCGR3A", "LYN"],
    "ID2-hi myeloid prog": [
        "CD14",
        "ID2",
        "VCAN",
        "S100A9",
        "CLEC12A",
        "KLF4",
        "PLAUR",
    ],
    "cDC1": ["CLEC9A", "CADM1"],
    "cDC2": [
        "CST3",
        "COTL1",
        "LYZ",
        "DMXL2",
        "CLEC10A",
        "FCER1A",
    ],  # Note: DMXL2 should be negative
    "Normoblast": ["SLC4A1", "SLC25A37", "HBB", "HBA2", "HBA1", "TFRC"],
    "Erythroblast": ["MKI67", "HBA1", "HBB"],
    "Proerythroblast": [
        "CDK6",
        "SYNGR1",
        "HBM",
        "GYPA",
    ],  # Note HBM and GYPA are negative markers
    "NK": ["GNLY", "NKG7", "CD247", "GRIK4", "FCER1G", "TYROBP", "KLRG1", "FCGR3A"],
    "ILC": ["ID2", "PLCG2", "GNLY", "SYNE1"],
    "Lymph prog": [
        "VPREB1",
        "MME",
        "EBF1",
        "SSBP2",
        "BACH2",
        "CD79B",
        "IGHM",
        "PAX5",
        "PRKCE",
        "DNTT",
        "IGLL1",
    ],
    "Naive CD20+ B": ["MS4A1", "IL4R", "IGHD", "FCRL1", "IGHM"],
    "B1 B": [
        "MS4A1",
        "SSPN",
        "ITGB1",
        "EPHA4",
        "COL4A4",
        "PRDM1",
        "IRF4",
        "CD38",
        "XBP1",
        "PAX5",
        "BCL11A",
        "BLK",
        "IGHD",
        "IGHM",
        "ZNF215",
    ],  # Note IGHD and IGHM are negative markers
    "Transitional B": ["MME", "CD38", "CD24", "ACSM3", "MSI2"],
    "Plasma cells": ["MZB1", "HSP90B1", "FNDC3B", "PRDM1", "IGKC", "JCHAIN"],
    "Plasmablast": ["XBP1", "RF4", "PRDM1", "PAX5"],  # Note PAX5 is a negative marker
    "CD4+ T activated": ["CD4", "IL7R", "TRBC2", "ITGB1"],
    "CD4+ T naive": ["CD4", "IL7R", "TRBC2", "CCR7"],
    "CD8+ T": ["CD8A", "CD8B", "GZMK", "GZMA", "CCL5", "GZMB", "GZMH", "GZMA"],
    "T activation": ["CD69", "CD38"],  # CD69 much better marker!
    "T naive": ["LEF1", "CCR7", "TCF7"],
    "pDC": ["GZMB", "IL3RA", "COBLL1", "TCF4"],
    "G/M prog": ["MPO", "BCL2", "KCNQ5", "CSF3R"],
    "HSC": ["NRIP1", "MECOM", "PROM1", "NKAIN2", "CD34"],
    "MK/E prog": [
        "ZNF385D",
        "ITGA2B",
        "RYR3",
        "PLCB1",
    ],  # Note PLCB1 is a negative marker
}

marker_genes = { "CD14+ Mono": ["FCN1", "CD14"], "CD16+ Mono": ["TCF7L2", "FCGR3A", "LYN"], "ID2-hi myeloid prog": [ "CD14", "ID2", "VCAN", "S100A9", "CLEC12A", "KLF4", "PLAUR", ], "cDC1": ["CLEC9A", "CADM1"], "cDC2": [ "CST3", "COTL1", "LYZ", "DMXL2", "CLEC10A", "FCER1A", ], # Note: DMXL2 should be negative "Normoblast": ["SLC4A1", "SLC25A37", "HBB", "HBA2", "HBA1", "TFRC"], "Erythroblast": ["MKI67", "HBA1", "HBB"], "Proerythroblast": [ "CDK6", "SYNGR1", "HBM", "GYPA", ], # Note HBM and GYPA are negative markers "NK": ["GNLY", "NKG7", "CD247", "GRIK4", "FCER1G", "TYROBP", "KLRG1", "FCGR3A"], "ILC": ["ID2", "PLCG2", "GNLY", "SYNE1"], "Lymph prog": [ "VPREB1", "MME", "EBF1", "SSBP2", "BACH2", "CD79B", "IGHM", "PAX5", "PRKCE", "DNTT", "IGLL1", ], "Naive CD20+ B": ["MS4A1", "IL4R", "IGHD", "FCRL1", "IGHM"], "B1 B": [ "MS4A1", "SSPN", "ITGB1", "EPHA4", "COL4A4", "PRDM1", "IRF4", "CD38", "XBP1", "PAX5", "BCL11A", "BLK", "IGHD", "IGHM", "ZNF215", ], # Note IGHD and IGHM are negative markers "Transitional B": ["MME", "CD38", "CD24", "ACSM3", "MSI2"], "Plasma cells": ["MZB1", "HSP90B1", "FNDC3B", "PRDM1", "IGKC", "JCHAIN"], "Plasmablast": ["XBP1", "RF4", "PRDM1", "PAX5"], # Note PAX5 is a negative marker "CD4+ T activated": ["CD4", "IL7R", "TRBC2", "ITGB1"], "CD4+ T naive": ["CD4", "IL7R", "TRBC2", "CCR7"], "CD8+ T": ["CD8A", "CD8B", "GZMK", "GZMA", "CCL5", "GZMB", "GZMH", "GZMA"], "T activation": ["CD69", "CD38"], # CD69 much better marker! "T naive": ["LEF1", "CCR7", "TCF7"], "pDC": ["GZMB", "IL3RA", "COBLL1", "TCF4"], "G/M prog": ["MPO", "BCL2", "KCNQ5", "CSF3R"], "HSC": ["NRIP1", "MECOM", "PROM1", "NKAIN2", "CD34"], "MK/E prog": [ "ZNF385D", "ITGA2B", "RYR3", "PLCB1", ], # Note PLCB1 is a negative marker }

将数据子集化，只保留在我们的数据中检测到的标记基因。我们将遍历所有细胞类型，只保留在我们的adata对象中发现的基因作为该细胞类型的标记。这将防止我们开始绘图时出现错误。

In [8]:

marker_genes_in_data = dict()
for ct, markers in marker_genes.items():
    markers_found = list()
    for marker in markers:
        if marker in adata.var.index:
            markers_found.append(marker)
    marker_genes_in_data[ct] = markers_found

marker_genes_in_data = dict() for ct, markers in marker_genes.items(): markers_found = list() for marker in markers: if marker in adata.var.index: markers_found.append(marker) marker_genes_in_data[ct] = markers_found

为了查看这些标记基因的表达位置，我们可以使用数据的二维可视化，如UMAP。我们将在此基于经过scran标准化的计数数据构建UMAP，仅使用高偏离基因。请注意，我们首先对标准化计数进行PCA以降低数据维度，然后生成UMAP。

首先，我们将原始计数存储在 .layers['counts'] 中，这样如果以后需要，我们仍然可以访问它们。然后，我们将 adata.X 设置为经过 scran 标准化和对数转换的计数。

In [9]:

adata.layers["counts"] = adata.X
adata.X = adata.layers["scran_normalization"]

adata.layers["counts"] = adata.X adata.X = adata.layers["scran_normalization"]

此外，我们将 adata.var.highly_variable 设置为高度偏离的基因。Scanpy 在后续计算中使用此 var 列，例如下面的 PCA。

In [10]:

adata.var["highly_variable"] = adata.var["highly_deviant"]

adata.var["highly_variable"] = adata.var["highly_deviant"]

现在执行PCA。我们使用高度偏离的基因（上述设置为“highly variable”）来减少噪音并增强数据中的信号，并将组件数量设置为默认值n=50。对于单个样本的数据来说，50是一个较高的数值，但它将确保我们不会忽略数据中的重要变化。

In [11]:

sc.tl.pca(adata, n_comps=50, use_highly_variable=True)

sc.tl.pca(adata, n_comps=50, use_highly_variable=True)

基于主成分（PCs）计算邻居图：

In [12]:

sc.pp.neighbors(adata)

sc.pp.neighbors(adata)

并使用该邻居图计算数据的二维 UMAP 嵌入：

In [13]:

sc.tl.umap(adata)

sc.tl.umap(adata)

现在使用计算得出的 UMAP 显示标记基因的表达。在此示例中，我们将自己限制在 B/浆细胞亚型上。请注意，上述标记字典中有三个负标记：B1 B 中的 IGHD 和 IGHM 以及浆母细胞中的 PAX5，这意味着预期这些细胞类型不表达或低水平表达这些标记。

让我们列出我们想要显示标记的 B 细胞亚型：

In [14]:

B_plasma_cts = [
    "Naive CD20+ B",
    "B1 B",
    "Transitional B",
    "Plasma cells",
    "Plasmablast",
]

B_plasma_cts = [ "Naive CD20+ B", "B1 B", "Transitional B", "Plasma cells", "Plasmablast", ]

现在为每个 B 细胞亚型的每个标记绘制一个 UMAP。请注意，我们只能绘制在我们的数据中存在的标记。

In [15]:

for ct in B_plasma_cts:
    print(f"{ct.upper()}:")  # print cell subtype name
    sc.pl.umap(
        adata,
        color=marker_genes_in_data[ct],
        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
    )
    print("\n\n\n")  # print white space for legibility

for ct in B_plasma_cts: print(f"{ct.upper()}:") # print cell subtype name sc.pl.umap( adata, color=marker_genes_in_data[ct], 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 ) print("\n\n\n") # print white space for legibility

NAIVE CD20+ B:

B1 B:

TRANSITIONAL B:

PLASMA CELLS:

PLASMABLAST:

正如您所看到的，即使是单个细胞类型的标记也经常在数据的不同子集中表达，即单个标记通常不会唯一地在一个细胞类型中表达。相反，正是这些子集的交集将告诉您感兴趣的细胞类型在哪里。

您可能还注意到，标记基因通常是稀疏表达的，即通常只有某个细胞类型的一部分细胞检测到某个标记。这是由于scRNA-seq数据的特性：我们只测序细胞中总RNA分子的一个小子集，由于这种抽样，有时即使某些基因在细胞中表达，我们也不会在该细胞中检测到这些基因的转录本。因此，我们不会基于一组标记的最小表达阈值来注释单个细胞。相反，我们首先通过聚类将数据细分为相似细胞的组（即“分区”数据），从而考虑单个基因的“缺失转录本”，而是基于总体转录组相似性进行分组。然后，我们可以基于这些组的整体标记表达模式来注释这些聚类。

现在让我们对数据进行聚类。我们将使用 Leiden 算法 {cite}anno:Traag2019，如聚类章节中所讨论的那样，将数据分组为相似的细胞子集：

In [16]:

sc.tl.leiden(adata, resolution=1, key_added="leiden_1")

sc.tl.leiden(adata, resolution=1, key_added="leiden_1")

In [17]:

sc.pl.umap(adata, color="leiden_1")

sc.pl.umap(adata, color="leiden_1")

您可能会注意到，这种数据分区相当粗略，一些我们之前看到的标记表达模式并未被此聚类捕获。因此，我们可以通过更改聚类的分辨率参数来尝试更高分辨率的聚类：

In [18]:

sc.tl.leiden(adata, resolution=2, key_added="leiden_2")

sc.tl.leiden(adata, resolution=2, key_added="leiden_2")

In [19]:

sc.pl.umap(adata, color="leiden_2")

sc.pl.umap(adata, color="leiden_2")

或者在 UMAP 中显示聚类编号：

In [20]:

sc.pl.umap(adata, color="leiden_2", legend_loc="on data")

sc.pl.umap(adata, color="leiden_2", legend_loc="on data")

这种聚类更加精细，有助于我们更详细地注释数据。您可以调整分辨率参数，以找到最佳捕捉到您观察到的标记表达模式的设置。

向上滚动，您会看到簇4和簇6是持续表达 Naive CD20+ B 细胞标记的簇。我们还可以使用点图来可视化这一点：

In [21]:

B_plasma_markers = {
    ct: [m for m in ct_markers if m in adata.var.index]
    for ct, ct_markers in marker_genes.items()
    if ct in B_plasma_cts
}

B_plasma_markers = { ct: [m for m in ct_markers if m in adata.var.index] for ct, ct_markers in marker_genes.items() if ct in B_plasma_cts }

In [22]:

sc.pl.dotplot(
    adata,
    groupby="leiden_2",
    var_names=B_plasma_markers,
    standard_scale="var",  # standard scale: normalize each gene to range from 0 to 1
)

sc.pl.dotplot( adata, groupby="leiden_2", var_names=B_plasma_markers, standard_scale="var", # standard scale: normalize each gene to range from 0 to 1 )

通过结合对 UMAP 和上述点图的视觉检查，我们现在可以开始对这些簇进行注释。

In [23]:

cl_annotation = {
    "4": "Naive CD20+ B",
    "6": "Naive CD20+ B",
    "8": "Transitional B",
    "18": "B1 B",  # note that IGHD and IGHM are negative markers, in this case more lowly expressed than in the other B cell clusters
}

cl_annotation = { "4": "Naive CD20+ B", "6": "Naive CD20+ B", "8": "Transitional B", "18": "B1 B", # note that IGHD and IGHM are negative markers, in this case more lowly expressed than in the other B cell clusters }

您可能会注意到，B1 B 细胞的注释比较困难，没有一个簇表达所有的 B1 B 标记，但有几个簇表达了一些标记。我们经常发现，适用于一个数据集的标记在其他数据集中效果不佳。这可能是由于测序深度的差异，也可能是由于数据集或样本之间的其他变异来源。

让我们可视化目前的注释结果：

In [24]:

adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)

adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)

In [25]:

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

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

... storing 'manual_celltype_annotation' as categorical

从簇差异表达基因到簇注释¶

相反，我们可以计算每个簇的标记基因，然后查找是否可以将这些标记基因与已知的生物学信息（如细胞类型和/或状态）关联起来。对于簇的标记基因计算，简单的方法如Wilcoxon秩和检验被认为表现最好 {cite}anno:Pullin2022.05.09.490241。重要的是，由于簇的定义基于用于这些统计检验的相同数据，这些检验的p值将被夸大，如这里所述 {cite}anno:ZHANG2019383。

让我们计算每个簇相对于我们adata中其他细胞的差异表达基因：

In [26]:

sc.tl.rank_genes_groups(
    adata, groupby="leiden_2", method="wilcoxon", key_added="dea_leiden_2"
)

sc.tl.rank_genes_groups( adata, groupby="leiden_2", method="wilcoxon", key_added="dea_leiden_2" )

我们可以使用标准的scanpy点图可视化每个簇的顶级差异表达基因的表达情况：

In [27]:

sc.pl.rank_genes_groups_dotplot(
    adata, groupby="leiden_2", standard_scale="var", n_genes=5, key="dea_leiden_2"
)

sc.pl.rank_genes_groups_dotplot( adata, groupby="leiden_2", standard_scale="var", n_genes=5, key="dea_leiden_2" )

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

正如您在上面看到的，许多差异表达基因在多个簇中都高度表达。我们可以过滤差异表达基因，以选择更多簇特异性的差异表达基因：

In [28]:

sc.tl.filter_rank_genes_groups(
    adata,
    min_in_group_fraction=0.2,
    max_out_group_fraction=0.2,
    key="dea_leiden_2",
    key_added="dea_leiden_2_filtered",
)

sc.tl.filter_rank_genes_groups( adata, min_in_group_fraction=0.2, max_out_group_fraction=0.2, key="dea_leiden_2", key_added="dea_leiden_2_filtered", )

可视化过滤后的基因：

In [29]:

sc.pl.rank_genes_groups_dotplot(
    adata,
    groupby="leiden_2",
    standard_scale="var",
    n_genes=5,
    key="dea_leiden_2_filtered",
)

sc.pl.rank_genes_groups_dotplot( adata, groupby="leiden_2", standard_scale="var", n_genes=5, key="dea_leiden_2_filtered", )

让我们看看簇12，它似乎有一组相对独特的标记基因，包括CDK6、ETV6、NKAIN2和GNAQ。通过搜索我们了解到，NKAIN2和ETV6是造血干细胞的标记基因 {cite}anno:SHI20222234 {cite}anno:Wang1998-rx（NKAIN2也出现在我们上面的列表中）。在UMAP中，我们可以看到这些基因在簇12中广泛表达：

In [30]:

sc.pl.umap(
    adata,
    color=["CDK6", "ETV6", "NKAIN2", "GNAQ", "leiden_2"],
    vmax="p99",
    legend_loc="on data",
    frameon=False,
    cmap="Reds",
)

sc.pl.umap( adata, color=["CDK6", "ETV6", "NKAIN2", "GNAQ", "leiden_2"], vmax="p99", legend_loc="on data", frameon=False, cmap="Reds", )

然而，查看已知的巨核细胞/红细胞祖细胞（"MK/E prog"）标记时，我们发现簇12的一部分似乎属于这种细胞类型：

In [31]:

sc.pl.umap(
    adata,
    color=[
        "ZNF385D",
        "ITGA2B",
        "RYR3",
        "PLCB1",
    ],
    vmax="p99",
    legend_loc="on data",
    frameon=False,
    cmap="Reds",
)

sc.pl.umap( adata, color=[ "ZNF385D", "ITGA2B", "RYR3", "PLCB1", ], vmax="p99", legend_loc="on data", frameon=False, cmap="Reds", )

这凸显了基于标记的注释的复杂性：它对您选择的聚类分辨率、标记集的稳健性和独特性以及您对数据中预期细胞类型的了解都非常敏感。

因此，该领域部分地在尝试摆脱手动聚类注释，而转向自动化注释算法。本教程的其余部分将重点介绍这些选项。

在继续之前，将最后的注释信息存储在我们的adata中：

In [32]:

cl_annotation["12"] = "HSCs + MK/E prog (?)"

cl_annotation["12"] = "HSCs + MK/E prog (?)"

In [33]:

adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)

adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)

Automated annotation¶

General remarks¶

接下来讨论的方法将是用于自动化数据注释的方法，而不是手动注释。与上述展示的方法不同，每种方法都可以以自动化方式注释数据。它们基于不同的原理，有时需要预定义的标记集，有时则在预先存在的完整scRNA-seq数据集上进行训练。正如下面所讨论的，生成的注释质量可能会有所不同。因此，重要的是将这些方法视为注释过程的起点而不是终点。有关自动化注释方法的更详细讨论，请参阅几篇评论 {cite}anno:PASQUINI2021961，{cite}anno:Abdelaal2019。

正如所说，自动生成的注释质量可能会有所不同。更具体地说，注释的质量取决于：

1. 选择的分类器类型：先前的基准研究表明，不同类型的分类器通常表现相当，基于神经网络的方法通常并不优于通用模型如支持向量机或线性回归模型 {cite}anno:Abdelaal2019, {cite}anno:PASQUINI2021961, {cite}anno:Huang2021。

2. 分类器训练数据的质量：如果训练数据没有得到很好的注释或注释分辨率较低，分类器也会有同样的问题。同样，如果训练数据和/或其注释是噪声较大的，分类器的性能可能会较差。

3. 您自己的数据与分类器训练数据的相似性：例如，如果分类器是在 drop-seq 单细胞数据集上训练的，而您的数据是 10X 单核而非单细胞 drop-seq，这可能会降低注释质量。在包括多样数据集的跨数据集图谱上训练的分类器可能比在单一数据集上训练的分类器提供更稳健和高质量的注释。一个例子是 CellTypist（一个将更详细讨论的自动注释方法）分类器，该分类器在包括 14 个不同肺数据集的人类肺细胞图谱上训练 {cite}anno:Sikkema2023。这个模型在新的肺数据上可能表现得比在单一肺数据集上训练的模型更好。

上述几点强调了使用分类器可能存在的缺点，具体取决于训练数据和模型类型。然而，使用预训练分类器来注释数据也有几个重要的优势。首先，这是一个快速且简单的方法来注释您的数据。注释不需要下载或预处理训练数据，有时只需将您的数据上传到一个在线网页。其次，这些方法不依赖于将数据划分为簇，像手动注释那样。第三，预训练分类器使您能够直接利用先前研究中的知识和信息，例如高质量的注释。最后，使用这样的分类器有助于在整个领域内统一细胞类型定义，从而为这些定义的领域共识铺平道路。

最后，由于这些分类器通常不像手动标记基因注释那样透明，因此一个好的不确定性度量来量化注释不确定性将提高方法的质量和可用性。我们将在下面更详细地讨论这一点。

Marker gene-based classifiers¶

一类自动细胞类型注释方法依赖于预定义的标记基因集。细胞根据这些标记基因的表达水平被分类为不同的细胞类型。这类方法的例子包括 Garnett {cite}anno:Pliner2019 和 CellAssign {cite}anno:Zhang2019。这些模型所基于的标记基因集越稳健和通用，模型的性能就会越好。然而，与其他模型一样，它们可能会受到模型训练数据与需要标注数据之间的批次效应相关差异的影响。与基于较大基因集的模型（见下文）相比，这些方法的一个优势是它们更加透明：我们知道分类是基于哪些基因进行的。
我们不会在这个笔记本中展示基于标记的分类器的示例，但如果您感兴趣，我们鼓励您自己探索这些方法。

Classifiers based on a wider set of genes¶

值得注意的是，目前讨论的方法仅使用了数据中检测到的基因的小子集：通常每种细胞类型使用的标记基因集只有1到大约10个。另一种方法是使用输入更大基因集（数千个或更多基因）的分类器，从而更充分地利用scRNA-seq数据的广度。这类分类器在先前注释的数据集或图谱上进行训练。这些方法的例子包括CellTypist {cite}anno:Conde2022（参见 https://www.celltypist.org，可以将数据上传到门户网站以获取自动细胞注释）和Clustifyr {cite}anno:Fu2020。

让我们在数据上尝试使用 CellTypist。根据 CellTypist 教程（https://www.celltypist.org/tutorials），我们需要准备数据，将计数标准化为每个细胞10,000计数，然后进行log1p转换：

In [34]:

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

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 模型：

In [35]:

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

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

In [36]:

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

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

对于每个模型，我们可以查看它包含哪些细胞类型，以确认是否包括骨髓细胞类型：

In [37]:

model_high.cell_types

model_high.cell_types

Out[37]:

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)

In [38]:

model_low.cell_types

model_low.cell_types

Out[38]:

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)

看起来这些模型包括了许多不同的免疫细胞类型前体！

现在让我们运行这些模型。首先是粗粒度模型：

In [39]:

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

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以获取完整输出：

In [40]:

predictions_high_adata = predictions_high.to_adata()

predictions_high_adata = predictions_high.to_adata()

...并将结果复制到我们原始的 AnnData 对象中：

In [41]:

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

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

现在对更细的注释进行相同的操作：

In [42]:

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

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!

In [43]:

predictions_low_adata = predictions_low.to_adata()

predictions_low_adata = predictions_low.to_adata()

In [44]:

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

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

Now plot:

In [45]:

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

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

In [46]:

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

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

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

In [47]:

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

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.

Out[47]:

<Axes: >

该树状图部分反映了关于细胞类型关系的先验知识（例如，B 细胞大部分聚集在一起），但我们也观察到一些意外的模式：Tcm/Naive 辅助性 T 细胞与红系细胞和巨噬细胞聚集在一起，而不是与其他 T 细胞聚集在一起。这是一个警示信号！可能 Tcm/Naive 辅助性 T 细胞的注释是错误的。

现在让我们看看之前的手动注释：

In [48]:

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

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

可以看到，我们的初始 B 细胞注释与部分自动化的初始 B 细胞注释很好地对应。类似地，我们称为过渡 B 细胞的一部分在他们的注释中被称为“小前 B 细胞”，而我们的 B1 B 细胞与他们的记忆 B 细胞对应，这令人鼓舞！

然而，您也会注意到我们的 HSC + MK/E prog 簇在他们的细注释中被注释为 T 细胞和 HSCs/多能祖细胞的混合体，因此这些注释部分是矛盾的。查看两个注释的置信分数，我们发现大部分细胞的注释都是在相对低置信度下完成的，这表明这些注释不能在没有仔细验证和手动审核的情况下直接采用！

以下是簇12在细化的CellTypist标签方面的细分：

In [49]:

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

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

Out[49]:

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 细胞或巨核细胞前体，这再次仅部分对应于我们的注释。

In [50]:

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

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

Out[50]:

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”（根据上面的树状图已标记为可能错误注释）的细胞实际上来自红系细胞系！

In [51]:

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

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

Annotation by mapping to a reference¶

一种注释数据的最终方法是基于将数据映射到现有的、已注释的单细胞参考上，然后使用生成的联合嵌入来执行标签转移。这个参考可以是你之前手动注释的单个样本，然后你希望将这些注释转移到其余的数据集中。或者，它可以是已发布的、理想情况下经过精心整理的现有参考。在这种情况下，我们将“新数据”，即需要映射和注释的数据，称为“查询（query）”。
目前有多种现有方法可以执行这种“查询到参考的映射”，包括scArches {cite}anno:Lotfollahi2022、Symphony {cite}anno:Kang2021和Azimuth (Seurat) {cite}anno:HAO20213573。所有这些方法都使您能够将新数据集映射到现有参考上，而无需重新整合参考数据，也无需访问完整的参考数据。

由于查询到参考的映射涉及将新数据嵌入到现有的低维参考数据表示中，这种低维表示的维度和轴在从查询中学习之前基本上是预定义的。因此，学习和整合查询中可能存在的未见过的变异（包括新的生物变异，如未见过的细胞类型或状态，以及新的技术变异，即需要移除的未见过的批次效应）对这些模型来说可能是一个挑战。结果，查询数据与参考数据的整合可能并不总是最优的，批次效应也可能无法完全从联合的查询-参考嵌入中移除。然而，由于细胞类型标签转移不一定需要完美的整合，只需要在嵌入中相同细胞类型的接近距离，即使是一个不完美的映射在注释你的数据时也仍然非常有帮助。

以scArches为例，它是一种基于参考映射的标签转移方法，其基础是一个现有的（基于变分自编码器的）模型，该模型将参考数据嵌入到一个低维的、经过批次校正的空间中。然后，它略微扩展了该模型，使其能够将未见过的数据集映射到相同的“潜在空间”（即低维嵌入）。这种模型扩展还能够学习和移除映射数据集中存在的批次效应。

现在我们将展示如何使用scArches将数据映射到参考中，并使用此映射将标签从参考转移到新数据（“查询”）。

{admonition}
请注意，如果您没有访问GPU，scArches可能无法运行，或者运行非常缓慢。因此，您可能需要在计算集群/服务器上运行此部分笔记本。

首先，让我们准备数据以进行映射到参考中。scArches方法允许我们将现有的参考模型适应到新数据上，需要原始的、未标准化的计数。因此，我们将保留计数层，并从待映射的adata中移除所有其他层。我们还将 .X 设置为这些原始计数。

In [52]:

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

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

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

In [55]:

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

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

该表包含基因名称和基因ID。由于基因ID通常在基因组注释版本中更少发生变化，我们将使用这些基因ID对数据进行子集化。因此，我们将为adata和参考模型特征都设置基因ID作为行名称。重要的是，我们必须确保也存储基因名称以供以后使用：基因名称比基因ID更容易理解。

In [56]:

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

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

In [57]:

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

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

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

In [58]:

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

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

Total number of genes needed for mapping: 4000

In [59]:

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

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对象连接起来。

In [60]:

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

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

In [61]:

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

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

将原始adata与缺失基因的adata连接起来。为了确保我们可以无错误地进行此连接，我们将从varm中删除PCA矩阵。

In [62]:

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

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

In [63]:

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 = sc.concat( [adata_to_map, missing_gene_adata], axis=1, join="outer", index_unique=None, merge="unique", )

现在，将数据子集化为模型中使用的基因，并按正确顺序排列：

In [64]:

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

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

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

In [65]:

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

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

Out[65]:

True

现在我们可以将基因索引重新设置为基因名称，以便于理解：

In [66]:

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_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']作为批次变量。因此，我们将检查是否在整个样本中将其设置为一个值：

In [67]:

adata_to_map_augmented.obs.batch.unique()

adata_to_map_augmented.obs.batch.unique()

Out[67]:

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

现在让我们讨论一下参考模型。参考模型越好，标签转移的效果就越好。使用集成了许多不同数据集并与您的数据匹配良好的（相同器官、相同单细胞技术等）精心注释的参考模型是理想的选择：这种模型在各种数据集和批次上训练，因此预计对批次效应更具鲁棒性。然而，并非所有组织都存在这样的参考模型。在本教程中，我们将使用一个参考模型，该模型是用我们在本书中使用的骨髓样本（排除将要映射的样本）训练的。该参考模型是一个scvi模型（用于数据集成），生成输入数据的低维集成嵌入，详情请参阅scvi的发表文章{cite}anno:Lopez2018-zc。请注意，这只是为本教程生成的示例模型，不应在其他情况下使用。
我们将从加载模型并将其传递给我们要映射的adata开始：

In [68]:

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

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

In [69]:

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

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对查询数据进行训练：

In [70]:

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

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]

现在我们已经更新了模型，可以计算查询数据的（理想情况下经过批次校正的）潜在表示：

In [71]:

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

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

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

In [72]:

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

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

为了查看基于映射的UMAP是否具有一般意义，让我们看看几个标记物，并查看它们的表达是否集中在UMAP的特定部分：

In [73]:

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
)

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 )

现在的关键步骤是我们可以将查询数据的推断潜在空间嵌入与现有的参考嵌入结合起来。使用这个联合嵌入，我们不仅可以一起可视化和聚类这两个数据集，还可以从查询到参考进行标签转移。
让我们加载参考嵌入：现有图谱通常会公开提供这些嵌入。

In [74]:

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

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

我们将存储一个变量，以指定这些细胞来自参考数据。

In [75]:

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

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

让我们看看这个参考对象中有什么：

In [76]:

ref_emb

ref_emb

Out[76]:

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个维度（在 .X 中），这些维度一起代表了参考细胞的潜在空间嵌入。我们为自己的数据计算的查询嵌入也有10个维度。参考和查询的10个维度是相同的，可以合并！
此外，它在 .obs['cell_type'] 中有细胞类型标签。我们将使用这些标签来注释我们自己的数据。

为了进行标签转移，我们将首先使用10维嵌入将参考数据和查询数据连接起来。为此，我们将从查询数据中创建与参考数据相同类型的AnnData对象（嵌入在 .X 中），并将两者连接起来。这样，我们可以联合分析参考和查询数据，包括从一个数据集到另一个数据集的转移。

In [77]:

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

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

In [78]:

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

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

In [79]:

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

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

让我们使用UMAP可视化联合嵌入。

In [80]:

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

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

我们可以基于UMAP直观地获得一个关于参考数据和查询数据是否良好集成的初步印象：

In [81]:

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

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

... storing 'reference_or_query' as categorical

在这个UMAP中查询数据和参考数据的（部分）混合是一个好兆头！当映射完全失败时，通常会在UMAP中看到查询数据和参考数据的完全分离。

现在让我们看看参考数据的细胞类型注释。查询数据中的所有细胞在这里都设置为NA，因为它们还没有注释，并显示为黑色。

我们将把这个图放大一些，以便更好地阅读图例：

In [82]:

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

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

In [83]:

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

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的标签转移。

First we set up the label transfer model:

In [84]:

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

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

现在我们进行标签转移：

In [85]:

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

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中：

In [86]:

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

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

让我们将结果转移到查询的adata对象中，该对象还包含我们的UMAP和基因计数，这样我们就可以将所有这些一起可视化。

In [87]:

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

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

We can now visualize the transferred labels in our previously calculated UMAP of our own data:

Let's set the figure size smaller again:

In [88]:

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

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

In [89]:

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

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

... storing 'transf_cell_type' as categorical

基于每个查询细胞的邻居，我们不仅可以猜测这些细胞所属的细胞类型，还可以生成一个衡量标签确定性的指标：如果一个细胞的邻居来自几种不同的细胞类型，我们的猜测将非常不确定。这对于评估我们可以在多大程度上“信任”转移的标签是很重要的！让我们可视化不确定性评分：

In [90]:

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

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

让我们检查每个细胞类型标签的标签转移不确定性水平。这可以让我们初步了解哪些注释更有争议/需要更多的人工检查。

In [91]:

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

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

你会注意到，例如，前体细胞通常比其他细胞类型更难区分。同样，我们注释中的"其他T细胞"类别也相对不具体。在最右边，我们可以看到pDCs，这是一种已知在转录上非常不同的细胞类型，因此更容易识别和标注。

为了在我们的转移标签中纳入这种不确定性信息，我们可以将不确定性评分高于0.2的细胞设置为“未知”：

"我们可以将不确定性评分高于0.2的细胞设置为“未知”"

In [92]:

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"

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中的“未知”颜色。

In [93]:

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

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

... storing 'transf_cell_type_certain' as categorical

为了便于阅读，我们可以只为“未知”细胞着色。这将使我们更容易看到这些细胞的数量。您也可以对其他任何细胞类型标签进行相同操作。

In [94]:

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

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

这些细胞数量相当多！这些细胞需要特别仔细的人工审核。然而，围绕“未知细胞”的低不确定性注释已经可以让我们初步了解每个细胞可能属于的细胞类型。

现在让我们看看更确定的注释。我们将检查一些细胞类型（这里随机选择）以了解参考转移的注释与我们上面已知的标记基因的匹配程度。实际上，这应该系统地对所有注释进行检查！

In [95]:

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

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

方便的是，对于每种细胞类型，我们的字典中都有相应的标记基因。让我们为所有新注释的细胞类型绘制标记基因的表达情况。您会注意到，标记基因的表达通常与自动注释相对应，这是一个好兆头！

In [96]:

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
)

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：

In [97]:

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

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

不确定性不仅帮助我们识别出算法对某个细胞类型所属的细胞不确定的区域（例如，因为它位于两个注释的表型之间），还可以突出未见过的细胞类型或新的细胞状态。例如，您的参考数据可能由健康细胞组成，而您的查询可能来自患病样本。不确定性评分可以突出疾病特异性的细胞状态，因为这些细胞可能没有来自参考数据的、始终来自单一细胞类型的邻居。尤其是当您的参考数据基于大量数据集时，不确定性评分有助于标记查询数据中可能值得研究的部分。因此，基于参考的标签转移不仅有助于注释您的数据，还可以加快数据的探索和解释。然而，像任何指标一样，这些不确定性评分往往并不完美，在某些情况下无法突出新细胞类型或状态。关于不确定性指标的更广泛讨论，请参见例如 {cite}anno:Engelmann2019。

如同本笔记本中讨论的任何方法一样，转移注释的质量取决于“训练数据”（在本例中是参考数据）及其注释的质量、模型的质量以及您自己的数据与训练数据的匹配程度！

因此，总是需要通过手动检查来验证转移注释的质量，并利用标记基因表达进行检查，可能还需要对初始注释进行细化。

Finally, store your adata object if wanted:

In [98]:

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

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