9. 细胞注释Annotation

为了更好地理解数据并利用现有知识，弄清楚数据中每个细胞的“细胞身份”非常重要。根据已知（或有时未知）的细胞表型标记数据中的细胞群的过程称为“细胞注释”。
那么什么是细胞类型呢？生物学家使用术语“细胞类型”来表示一种细胞表型，该表型在不同数据集中具有稳定性，可以根据特定标记（即蛋白质或基因转录本）的表达进行识别，并且通常与特定功能相关。例如，血浆B细胞是一种白细胞，可以分泌用于对抗病原体的抗体，可以使用特定标记进行识别。了解样品中的细胞类型对于理解数据至关重要。例如，了解肿瘤中存在特定的免疫细胞类型或骨髓样本中存在不寻常的造血干细胞可以为正在研究的疾病提供有价值的见解。
然而，与任何分类一样，类别的大小以及它们之间的边界在一定程度上是主观的，并且可能随着时间的推移而变化，利用新技术可以更高分辨率地观察细胞，或者将本来被认为不具有生物学意义的特定“亚表型”赋予重要的生物学意义。因此，细胞类型通常被进一步分为“亚型”或“细胞状态”（例如激活与静息），一些研究人员使用术语“细胞身份”来指代细胞类型、细胞亚型和细胞状态。
同样，多种细胞类型可以是单个连续体的一部分，其中一种细胞类型可能转变或分化为另一种细胞类型。例如，在造血细胞中，干细胞分化为特定的免疫细胞类型。尽管通常会在这种分化的早期和晚期阶段之间划出界限，但这些细胞的状态可以通过分化程度较低和分化程度较高的细胞表型之间的分化坐标来更准确地描述。我们将在后续章节中讨论分化和细胞轨迹。
那么我们如何注释单细胞数据中的细胞呢？有多种方法可以做到这一点，我们将在下面概述不同的方法。当我们处理转录组数据时，每种方法最终都基于特定基因或基因集的表达，或细胞之间的转录组相似性。

9.1 环境配置

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

import warnings

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

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

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

加载需要的模块：

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警告也过滤：

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

我们将继续对前期预处理的scRNA-seq数据集进行工作，并对其进行注释。
设置图形参数：

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

9.2 加载数据

使用数据为该教程前一部分也使用的数据中的单个样本('site4-donor8')。此外，未通过QC的细胞已经被去除。

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

9.3 手动注释

执行细胞类型注释的经典或最古老的方法是基于已知与特定细胞类型相关的单个或一小组标记基因。 这种方法可以追溯到“scRNA-seq时代之前”，当时单细胞数据是低维的（例如，FACS数据的基因组由不超过30-40个基因组成）。这是一种快速且透明的数据注释方式 然而，当特定细胞类型不存在独特的标记时（通常是这种情况），这种方法可能会变得更加复杂，甚至不太客观，需要结合正确注释所需的标记或表达阈值。一组强大的标记基因和先验知识或注释经验可以在此提供帮助，但该方法存在不清楚和主观决策的风险。
在这种设置中，数据通常在注释之前进行聚类，以便我们可以注释细胞群，而不是进行每个细胞的调用。这不仅不那么费力，而且可以更好的排除噪声干扰：单个细胞可能没有特定标记的计数，即使它在该细胞中表达，这仅仅是由于单细胞数据固有的稀疏性。聚类能够检测整体基因表达高度相似的细胞。
最后，有两个方法可以处理基于标记基因的注释。一种选择是根据数据中预期的所有细胞类型的标记基因表进行工作，并检查这些簇在哪些细胞中表达。另一种选择是检查哪些基因在您定义的簇中高度表达，然后检查它们是否与已知的细胞类型或状态相关。如有必要，可以在这些方法之间来回切换。 我们将在下面展示两者的示例。

9.3.1 从标记marker注释细胞簇

首先使用基于已知标记的方法。我们将首先列出一组基于文献的骨髓细胞类型标记：之前研究特定细胞类型和亚型并报告这些细胞类型的标记基因的论文。请注意，蛋白质水平的标记（例如用于 FACS）有时在转录组数据中效果不佳，因此使用基于RNA的论文中的标记通常更有可能发挥作用。此外，有时一个数据集中的标记在其他数据集中效果不佳。因此，理想情况下，标记集可以跨多个数据集进行验证。最后，寻求领域内专家合作通常很有用：作为生物信息学家，尝试与对组织、生物学、预期细胞类型和标记等有更广泛知识的生物学家合作。

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对象中找到的基因作为该细胞类型的标记。

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-normalized counts数据构建UMAP，仅使用高变的基因。请注意，在生成UMAP之前，我们首先对标准化计数执行PCA以降低数据的维数。
首先，我们将原始计数存储在.layers[‘counts’]中，以便以后在需要时仍然可以访问它们。然后，我们将adata.X设置为SCRAN-normalized、对数转换的counts。

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

我们进一步将adata.var.highly_variable设置为高变的基因。Scanpy在下游计算中使用这个var列，例如下面的PCA，

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

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

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

根据PC计算领域图：

sc.pp.neighbors(adata)

并使用该领域图来计算数据的二维UMAP嵌入：

sc.tl.umap(adata)

现在使用计算的UMAP显示标记的表达。在本例中，我们将限制为B/浆细胞亚型。请注意上面的标记字典，我们的列表中有3个阴性标记：B1 B细胞的IGHD和IGHM，以及浆母细胞的PAX5，这些marker预计不表达或低度表达。
让我们列出我们想要显示标记的B细胞亚型：

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

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

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算法将数据分组定义为相似的细胞子集：

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

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

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

或者将簇号显示在UMAP图中：

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

这种聚类更加精细，将帮助我们更详细地注释数据。可以使用分辨率参数来找到最能捕获观察到的标记表达模式的设置。
往会看，你会发现簇4和5是一致表达Naive CD20+ B细胞标记的簇。我们还可以使用点图来可视化这一点：

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
}

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的可视化结果和上面的点图，我们现在可以开始注释簇：

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标记，而有几个簇表达部分标记。我们经常看到适用于一个数据集的标记不适用于其他数据集。这可能是由于测序深度的差异，也可能是由于数据集或样本之间的其他变异来源造成的。
让我们可视化目前的注释结果：

adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)
sc.pl.umap(adata, color=["manual_celltype_annotation"])

输出结果：

... storing 'manual_celltype_annotation' as categorical

9.3.2 从簇差异表达基因注释细胞

相反，我们可以计算每个簇的标记基因，然后查找是否可以将这些标记基因与任何已知的生物学（例如细胞类型和/或状态）联系起来。对于簇的标记基因计算，简单的方法（例如 Wilcoxon 秩和检验）被认为效果最好。
让我们计算每个簇与数据中其他细胞相比的差异表达基因：

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

我们可以使用标准的scanpy点图可视化每个簇中最高差异表达基因的表达：

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.

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

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.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是造血干细胞标记物。在UMAP中，我们可以看到这些基因在整个簇12中都有表达：

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的部分似乎属于该细胞类型：

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

这体现了基于标记的注释是多么复杂：它对于选择的簇分辨率、所拥有的标记集的稳健性和唯一性以及对数据中预期的细胞类型的了解很敏感。
出于这个原因，该领域正在部分尝试摆脱手动注释，而转向细胞自动注释。
下一部分，我将重点介绍自动注释细胞类型算法。
但现在，我们需要先将最后的注释信息存储在我们的adata中：

cl_annotation["12"] = "HSCs + MK/E prog (?)"
adata.obs["manual_celltype_annotation"] = adata.obs.leiden_2.map(cl_annotation)