9. 細胞注釋Annotation

為了更好地理解數據并利用現有知識，弄清楚數據中每個細胞的“細胞身份”非常重要。根據已知（或有時未知）的細胞表型標記數據中的細胞群的過程稱為“細胞注釋”。
那么什么是細胞類型呢？生物學家使用術語“細胞類型”來表示一種細胞表型售担，該表型在不同數據集中具有穩(wěn)定性，可以根據特定標記（即蛋白質或基因轉錄本）的表達進行識別署辉，并且通常與特定功能相關族铆。例如，血漿B細胞是一種白細胞哭尝，可以分泌用于對抗病原體的抗體哥攘，可以使用特定標記進行識別。了解樣品中的細胞類型對于理解數據至關重要材鹦。例如逝淹，了解腫瘤中存在特定的免疫細胞類型或骨髓樣本中存在不尋常的造血干細胞可以為正在研究的疾病提供有價值的見解。
然而侠姑，與任何分類一樣创橄，類別的大小以及它們之間的邊界在一定程度上是主觀的，并且可能隨著時間的推移而變化莽红，利用新技術可以更高分辨率地觀察細胞，或者將本來被認為不具有生物學意義的特定“亞表型”賦予重要的生物學意義邦邦。因此安吁，細胞類型通常被進一步分為“亞型”或“細胞狀態(tài)”（例如激活與靜息），一些研究人員使用術語“細胞身份”來指代細胞類型燃辖、細胞亞型和細胞狀態(tài)鬼店。
同樣，多種細胞類型可以是單個連續(xù)體的一部分黔龟，其中一種細胞類型可能轉變或分化為另一種細胞類型妇智。例如滥玷，在造血細胞中，干細胞分化為特定的免疫細胞類型巍棱。盡管通常會在這種分化的早期和晚期階段之間劃出界限惑畴，但這些細胞的狀態(tài)可以通過分化程度較低和分化程度較高的細胞表型之間的分化坐標來更準確地描述。我們將在后續(xù)章節(jié)中討論分化和細胞軌跡航徙。
那么我們如何注釋單細胞數據中的細胞呢如贷？有多種方法可以做到這一點，我們將在下面概述不同的方法到踏。當我們處理轉錄組數據時杠袱，每種方法最終都基于特定基因或基因集的表達，或細胞之間的轉錄組相似性窝稿。

9.1 環(huán)境配置

我們將過濾掉一些不影響我們代碼的棄用和性能警告：

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)

我們將繼續(xù)對前期預處理的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 手動注釋

執(zhí)行細胞類型注釋的經典或最古老的方法是基于已知與特定細胞類型相關的單個或一小組標記基因潮梯。 這種方法可以追溯到“scRNA-seq時代之前”骗灶，當時單細胞數據是低維的（例如，FACS數據的基因組由不超過30-40個基因組成）秉馏。這是一種快速且透明的數據注釋方式 然而耙旦，當特定細胞類型不存在獨特的標記時（通常是這種情況），這種方法可能會變得更加復雜萝究，甚至不太客觀免都，需要結合正確注釋所需的標記或表達閾值。一組強大的標記基因和先驗知識或注釋經驗可以在此提供幫助帆竹，但該方法存在不清楚和主觀決策的風險绕娘。
在這種設置中，數據通常在注釋之前進行聚類栽连，以便我們可以注釋細胞群险领，而不是進行每個細胞的調用。這不僅不那么費力秒紧，而且可以更好的排除噪聲干擾：單個細胞可能沒有特定標記的計數绢陌，即使它在該細胞中表達，這僅僅是由于單細胞數據固有的稀疏性熔恢。聚類能夠檢測整體基因表達高度相似的細胞脐湾。
最后，有兩個方法可以處理基于標記基因的注釋叙淌。一種選擇是根據數據中預期的所有細胞類型的標記基因表進行工作秤掌，并檢查這些簇在哪些細胞中表達愁铺。另一種選擇是檢查哪些基因在您定義的簇中高度表達，然后檢查它們是否與已知的細胞類型或狀態(tài)相關闻鉴。如有必要茵乱，可以在這些方法之間來回切換。 我們將在下面展示兩者的示例椒拗。

9.3.1 從標記marker注釋細胞簇

首先使用基于已知標記的方法似将。我們將首先列出一組基于文獻的骨髓細胞類型標記：之前研究特定細胞類型和亞型并報告這些細胞類型的標記基因的論文。請注意蚀苛，蛋白質水平的標記（例如用于 FACS）有時在轉錄組數據中效果不佳在验，因此使用基于RNA的論文中的標記通常更有可能發(fā)揮作用。此外堵未，有時一個數據集中的標記在其他數據集中效果不佳腋舌。因此，理想情況下渗蟹，標記集可以跨多個數據集進行驗證块饺。最后，尋求領域內專家合作通常很有用：作為生物信息學家雌芽，嘗試與對組織授艰、生物學、預期細胞類型和標記等有更廣泛知識的生物學家合作世落。

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
}

僅保留在我們的數據中檢測到的標記的子集淮腾。我們將循環(huán)遍歷所有細胞類型，并僅保留在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之前专钉，我們首先對標準化計數執(zhí)行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"]

現在執(zhí)行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分子總量的一小部分進行測序，并且由于這種二次采樣审编，我們有時不會對細胞中特定基因的轉錄本進行采樣撼班，即使它們已表達在那個細胞里。因此垒酬，我們不會根據一組標記的最小表達閾值來注釋單個細胞砰嘁。相反，我們首先通過聚類將數據細分為相似細胞組（即“分區(qū)”數據）勘究，從而解釋單個基因的“缺失轉錄本”矮湘，而不是根據整體轉錄組相似性進行分組。然后我們可以根據它們的整體標記表達模式來注釋這些簇口糕。
現在讓我們對數據進行聚類缅阳。我們將使用“聚類”章節(jié)中討論的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")

這種聚類更加精細，將幫助我們更詳細地注釋數據向族『茄啵可以使用分辨率參數來找到最能捕獲觀察到的標記表達模式的設置。
往會看件相，你會發(fā)現簇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 從簇差異表達基因注釋細胞

相反，我們可以計算每個簇的標記基因辩稽，然后查找是否可以將這些標記基因與任何已知的生物學（例如細胞類型和/或狀態(tài)）聯系起來惧笛。對于簇的標記基因計算，簡單的方法（例如 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”）的已知標記了赵，我們發(fā)現簇12的部分似乎屬于該細胞類型：

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

這體現了基于標記的注釋是多么復雜：它對于選擇的簇分辨率潜支、所擁有的標記集的穩(wěn)健性和唯一性以及對數據中預期的細胞類型的了解很敏感。
出于這個原因柿汛，該領域正在部分嘗試擺脫手動注釋冗酿，而轉向細胞自動注釋。
下一部分，我將重點介紹自動注釋細胞類型算法已烤。
但現在鸠窗，我們需要先將最后的注釋信息存儲在我們的adata中：

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