单细胞 RNA 测序分析的当前最佳实践：教程-文献精读80

Current best practices in single‐cell RNA‐seq analysis: a tutorial

单细胞 RNA 测序分析的当前最佳实践：教程

摘要

单细胞 RNA 测序使基因表达的研究达到了前所未有的分辨率。这项技术的前景吸引了越来越多的用户使用单细胞分析方法。随着更多分析工具的出现，如何在这个领域中找到合适的工具并建立一个最新的分析流程变得愈发困难。在本文中，我们详细介绍了一个典型的单细胞 RNA 测序分析的步骤，包括预处理（质量控制、归一化、数据校正、特征选择和降维）以及细胞和基因层面的下游分析。我们基于独立的比较研究，为这些步骤制定了当前的最佳实践建议。我们将这些最佳实践建议整合到一个工作流程中，并应用于一个公开数据集，以进一步展示这些步骤在实际中的应用。我们的案例研究文档可以在 https://www.github.com/theislab/single-cell-tutorial 找到。本综述将作为该领域新入门者的工作流程教程，并帮助已有经验的用户更新他们的分析流程。

引言

近年来，单细胞 RNA 测序（scRNA-seq）显著推动了我们对生物系统的理解。我们不仅能够研究斑马鱼、青蛙和扁虫的细胞异质性（Briggs 等, 2018; Plass 等, 2018; Wagner 等, 2018），还发现了此前被忽视的细胞群体（Montoro 等, 2018; Plasschaert 等, 2018）。这种技术的巨大潜力激励了计算生物学家开发出多种分析工具（Rostom 等, 2017）。尽管该领域已经投入了相当多的努力来提高单个工具的可用性，但由于该领域相对不成熟，缺乏标准化，这对新手来说是一个门槛。在本文中，我们为 scRNA-seq 分析提供了教程，并概述了当前的最佳实践，以为未来的分析标准化奠定基础。

标准化的挑战包括分析方法数量的增长（截至2019年3月7日已有385种工具）以及数据集规模的迅速膨胀（Angerer 等, 2017; Zappia 等, 2018）。我们不断发现新的数据应用方式。例如，最近已经可以预测分化过程中细胞的命运（La Manno 等, 2018）。尽管分析工具的不断改进有助于产生新的科学见解，但也增加了标准化的难度。

标准化的另一个挑战在于技术层面。用于 scRNA-seq 数据的分析工具是用多种编程语言编写的，主要是 R 和 Python（Zappia 等, 2018）。虽然跨环境支持正在增长（预印本: Scholz 等, 2018），但编程语言的选择往往也意味着分析工具的选择。诸如 Seurat（Butler 等, 2018）、Scater（McCarthy 等, 2017）或 Scanpy（Wolf 等, 2018）等流行平台提供了集成的环境来开发分析流程，并包含大量的分析工具箱。然而，出于必要，这些平台通常只限于使用其相应编程语言开发的工具。当前可用的 scRNA-seq 分析教程也受到语言限制，许多教程围绕上述平台展开（R 和 Bioconductor 工具: https://github.com/drisso/bioc2016singlecell 和 https://hemberg-lab.github.io/scRNA.seq.course/; Lun 等, 2016b; Seurat: Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurat; Scanpy: Tutorials — scanpy）。

鉴于上述挑战，我们不试图建立一个标准化的分析流程，而是概述了当前的最佳实践和常用工具，独立于编程语言。我们将引导读者逐步完成 scRNA-seq 分析流程的各个步骤（图1），介绍当前的最佳实践，并讨论分析中的陷阱和未解决的问题。在由于工具的新颖性和缺乏比较而无法确定最佳实践的地方，我们列出了一些常用的工具。我们概述的步骤从读取或计数矩阵开始，最终达到潜在的分析终点。较早的预处理步骤已在 Lun 等人（2016b）中介绍。一个整合了当前最佳实践的详细案例研究可在我们的 GitHub 上找到：https://github.com/theislab/single-cell-tutorial/。在这里，我们将当前的最佳实践应用于实际的工作流程示例，以分析一个公开数据集。该分析流程在一个包含 rpy2 的 Jupyter–Ipython notebook 中集成了 R 和 Python 工具。通过现有文档，该流程可作为模板轻松适应和使用。

图1. 典型单细胞 RNA 测序分析流程示意图

原始测序数据经过处理和比对，生成计数矩阵，作为工作流程的起点。计数数据随后进行预处理和下游分析。子图是基于 Haber 等人（2017）的肠上皮数据，使用最佳实践工作流程生成的。

预处理与可视化

测序设备生成的原始数据经过处理后，得到分子计数矩阵（计数矩阵）或读取计数矩阵（读取矩阵），具体取决于在单细胞文库构建协议中是否使用了唯一分子标识符（UMIs）（参见框1，了解分析之前的实验步骤概览）。诸如 Cell Ranger（Zheng 等, 2017）、indrops（Klein 等, 2015）、SEQC（Azizi 等, 2018）或 zUMIs（Parekh 等, 2018）等原始数据处理流程涵盖了读取质量控制（QC）、将读取分配给其细胞条形码和 mRNA 分子来源（也称为“去重”）、基因组比对和定量。生成的读取或计数矩阵的维度为条形码数量 x 转录本数量。这里使用“条形码”一词而非“细胞”，因为所有分配到同一条形码的读取不一定来自同一细胞。一个条形码可能会错误地标记多个细胞（双重标签）或没有标记任何细胞（空滴/孔）。

box1：实验性 scRNA-seq 工作流程的关键要素

从生物样本生成单细胞数据需要多个步骤。典型工作流程包括单细胞分离、单细胞隔离、文库构建和测序。我们在此简要概述这些阶段。关于不同协议的详细解释和比较，请参阅 Ziegenhain 等（2017）；Macosko 等（2015）；Svensson 等（2017）。

• 单细胞实验的输入材料通常来自生物组织样本。第一步是通过单细胞分离生成单细胞悬液，在此过程中，组织被消化分解。

• 为了单独分析每个细胞中的 mRNA，必须对细胞进行隔离。单细胞隔离方式因实验协议而异。板基技术将细胞隔离到板上的孔中，而液滴基方法依赖于将每个细胞捕获在自己的微流体液滴中。无论哪种方法，都可能出现错误，导致多个细胞被一起捕获（双重或多重标签）、不可存活的细胞被捕获，或没有捕获任何细胞（空滴/孔）。空滴在液滴基方法中特别常见，因为该方法依赖于低浓度细胞流来控制双重标签的发生率。

• 每个孔或液滴都包含必要的化学物质来破坏细胞膜并进行文库构建。文库构建是捕获细胞内 mRNA，将其逆转录为 cDNA 分子并扩增的过程。由于每个细胞在隔离的环境中完成此过程，每个细胞的 mRNA 可以被标记上特定于孔或液滴的条形码。此外，许多实验协议还会使用唯一分子标识符（UMI）来标记捕获的分子。为了增加检测概率，细胞的 cDNA 在测序前被扩增。UMI 使我们能够区分相同 mRNA 分子的扩增拷贝和来自同一基因的独立 mRNA 分子的读取。

• 在文库构建之后，细胞 cDNA 文库被标记上细胞条形码，并根据协议选择性地加入 UMI。这些文库被合并（多重标记）以进行测序。测序生成的读取数据会经历质量控制、根据分配的条形码分组（去重）和读取处理流程中的比对。对于基于 UMI 的协议，读取数据可以进一步去重以生成捕获 mRNA 分子的计数数据（计数数据）。

尽管读取数据和计数数据在测量噪声水平上有所不同，典型分析流程中的处理步骤是相同的。为简便起见，在本教程中我们将数据称为计数矩阵。当读取矩阵和计数矩阵的结果有所不同时，会特别提及读取矩阵。

质量控制

在分析单细胞基因表达数据之前，我们必须确保所有的细胞条形码数据都对应于活细胞。细胞质量控制（QC）通常基于三个 QC 协变量进行：每个条形码的计数数量（计数深度）、每个条形码的基因数量，以及来自线粒体基因的计数在每个条形码中的比例（Ilicic 等, 2016；Griffiths 等, 2018）。这些 QC 协变量的分布被检查以寻找异常峰值，通过设定阈值来过滤掉这些异常值（图2）。这些异常条形码可能对应于濒死细胞、细胞膜破裂的细胞或双重标签。例如，具有低计数深度、检测到的基因数量少且线粒体计数比例高的条形码通常表明细胞的胞质 mRNA 已经通过破裂的膜泄漏出去，因此只有位于线粒体中的 mRNA 仍然保留（图2）。相反，具有意外高计数和大量检测到的基因的细胞可能表示双重标签。因此，通常使用高计数深度阈值来过滤潜在的双重标签。

近来的三个双重标签检测工具提供了更优雅且潜在更好的解决方案（DoubletDecon: 预印本: DePasquale 等, 2018；Scrublet: Wolock 等, 2019；Doublet Finder: McGinnis 等, 2018）。

图2. 来自 Haber 等人（2017）的鼠肠上皮数据集的质量控制指标及筛选决策的图示

(A) 每个细胞的计数深度直方图。较小的直方图放大了计数深度低于 4000 的区域。在约 1200 计数处检测到峰值，根据此设置了 1500 的阈值。(B) 每个细胞检测到的基因数量的直方图。在大约 400 个基因处可见一个小的噪声峰值。这些细胞通过在 700 个基因处设定的阈值（红线）被过滤掉。(C) 从高到低计数深度的分布。此可视化类似于 Cell Ranger 输出中的对数–对数图，用于筛选空滴。它显示了在约 1500 计数处计数深度开始迅速下降的“拐点”。(D) 基因数量与计数深度的关系图，按线粒体读取的比例着色。线粒体读取比例仅在检测到少量基因的低计数细胞中特别高。这些细胞通过我们的计数和基因数量阈值被过滤掉。联合可视化计数和基因阈值显示出联合筛选效果，表明较低的基因阈值可能已足够。[2019年7月5日，首次在线发布后添加更正：面板 B 中的 x 轴标签从“计数深度”更正为“基因数量”。]

单独考虑任何一个 QC 协变量可能会导致对细胞信号的误解。例如，线粒体计数比例较高的细胞可能参与呼吸过程。同样，其他 QC 协变量也具有生物学意义。低计数和/或基因数量的细胞可能对应于静止的细胞群体，而高计数的细胞可能体积更大。确实，不同细胞的分子计数可能差异很大（参见项目 GitHub 上的案例研究）。因此，在做单变量阈值决策时，QC 协变量应联合考虑（图 2D），并尽可能设定宽松的阈值，以避免意外过滤掉活细胞群体。未来，考虑多变量 QC 依赖关系的过滤模型可能提供更灵敏的 QC 选项。

包含异质细胞类型混合的数据集可能会出现多个 QC 协变量峰。例如，图 2D 显示了具有不同 QC 分布的两个细胞群体。如果没有进行之前的筛选步骤（请注意，Cell Ranger 也执行细胞 QC），则应仅将最低计数深度和条形码每基因的峰值视为非活细胞。另一个筛选指南是所选阈值下被筛除的细胞比例。对于高计数筛选，此比例不应超过预期的双重标签率。

除了检查细胞的完整性外，QC 步骤还必须在转录本层面进行。原始计数矩阵通常包含超过 20,000 个基因。通过过滤掉在少数细胞中表达的非信息性基因，可以显著减少这个数量。设置该阈值的指南是使用感兴趣的最小细胞簇大小，并为掉失效应留有余地。例如，过滤掉在少于 20 个细胞中表达的基因可能会使检测少于 20 个细胞的细胞簇变得困难。对于掉失率较高的数据集，此阈值也可能使检测更大的簇变得复杂。阈值的选择应随数据集中的细胞数量和预期的下游分析而变化。

进一步的 QC 可以直接在计数数据上进行。环境基因表达是指不来源于条形码细胞的计数，而是来自文库构建前被其他裂解细胞污染的 mRNA。这些额外的环境计数可能会扭曲下游分析，例如标志基因识别或其他差异表达测试，尤其是在不同样本间水平变化时。由于大量空滴存在于液滴基 scRNA-seq 数据集中，可以利用这些空滴建模环境 RNA 表达谱以校正这些效应。最近开发的 SoupX（预印本: Young & Behjati, 2018）使用此方法直接校正计数数据。忽略下游分析中的强环境基因的实用方法也被用于解决该问题（Angelidis 等, 2019）。

质量控制旨在确保数据质量足以进行下游分析。由于“足够的数据质量”无法事先确定，其依据下游分析性能（如簇注释）来判断。因此，在分析数据时可能需要多次重新审视质量控制决策。通常，从宽松的 QC 阈值开始，并调查这些阈值的效果，然后再回去执行更严格的 QC 是有益的。这种方法对于包含异质细胞群的数据集尤其相关，其中细胞类型或状态可能被误解为低质量的异常细胞。在低质量数据集中，严格的 QC 阈值可能是必要的。数据集的质量可以通过实验 QC 指标确定（参见附录补充文本 S2）。在这种迭代的 QC 优化中，应注意避免数据窥视。QC 阈值不应根据统计测试的结果进行调整。相反，可以通过数据集中 QC 协变量的分布和聚类在可视化中评估 QC。

陷阱与建议：

• 通过在基因数量、计数深度和线粒体读取比例上找到异常峰值来进行质量控制。联合考虑这些协变量，而不是单独考虑。

• 尽可能宽松地设置质量控制阈值；如果下游聚类难以解释，则重新审视质量控制。

• 如果样本间的质量控制协变量分布不同，应针对每个样本单独确定质量控制阈值，以考虑样本质量差异，如在 Plasschaert 等人（2018）中的做法。

归一化

计数矩阵中的每个计数代表一个细胞 mRNA 分子的成功捕获、逆转录和测序（参见框1）。相同细胞的计数深度可能由于这些步骤中的固有变异性而不同。因此，在基于计数数据比较细胞基因表达时，任何差异可能仅由于抽样效应引起。归一化通过缩放计数数据来解决此问题，以获得细胞间正确的相对基因表达丰度。

有许多用于整体基因表达的归一化方法（预印本：Pachter, 2011; Dillies 等, 2013）。虽然其中一些方法已应用于 scRNA-seq 分析，但单细胞数据特有的变异源（如技术掉失，即抽样导致的零计数）促使了专门针对 scRNA-seq 的归一化方法的发展（Lun 等, 2016a；Vallejos 等, 2017）。

最常用的归一化协议是计数深度缩放，也称为“每百万计数”（CPM）归一化。这种方法来自整体表达分析，并使用与每个细胞的计数深度成比例的大小因子对计数数据进行归一化。该方法的变体会用不同的 10 倍因子或数据集中每个细胞的中位计数深度来缩放大小因子。CPM 归一化假设数据集中所有细胞最初含有相同数量的 mRNA 分子，计数深度差异仅由于抽样引起。该假设与降采样协议一致，即从数据中随机抽取读取或计数，使所有细胞保留预设的最大计数或以下。虽然降采样丢弃了数据，但它也增加了技术掉失率，这是 CPM 和其他全局缩放归一化方法所没有的。因此，降采样可以更真实地反映在类似计数深度下的细胞表达谱。

由于单细胞数据集通常由具有不同大小和分子计数的异质细胞群体组成，通常更适合使用更复杂的归一化方法。例如，Weinreb 等人（2018）使用了 CPM 的一种简单扩展，当计算其大小因子时排除了占总计数至少 5% 的基因。这种方法允许少数高表达基因的分子计数变化。Scran 的基于合并的大小因子估计方法则允许更多的细胞异质性（Lun 等, 2016a）。在这里，大小因子通过对基因进行线性回归估计，在合并细胞以避免技术掉失效应后进行。该方法将细胞间的差异限制在不到 50% 的基因上，并在独立比较中始终表现优异。研究表明，Scran 在批次校正（Buttner 等, 2019）和差异表达分析（预印本：Vieth 等, 2019）中表现优于其他归一化方法。

CPM、高计数过滤 CPM 和 scran 使用线性、全局缩放方法来归一化计数数据。还有非线性归一化方法可以考虑更复杂的非期望变异（Cole 等, 2019）。许多这种方法涉及对计数数据的参数建模。例如，Mayer 等人（2018）将计数数据拟合为负二项式模型，使用技术协变量（如读取深度和每个基因的计数）来拟合模型参数。模型拟合的残差作为基因表达的归一化量化。这样的方法可以将技术和生物数据校正（例如批次校正或细胞周期效应校正）与计数深度归一化结合。研究表明，非线性归一化方法在存在强批次效应的情况下优于全局缩放方法（Cole 等, 2019）。

我们不能期望一种归一化方法适用于所有类型的 scRNA-seq 数据。例如，Vieth 等人（2017）表明读取和计数数据最好使用不同的模型来拟合。实际上，Cole 等人（2019）发现不同的归一化方法在不同数据集上表现最佳，并建议使用他们的 scone 工具来选择特定数据集的适当归一化方法。此外，scRNA-seq 技术可以分为全长和 3' 富集方法（Svensson 等, 2017; Ziegenhain 等, 2017）。全长协议的数据可能受益于考虑基因长度的归一化方法（例如 Patel 等, 2014; Kowalczyk 等, 2015; Soneson & Robinson, 2018），而 3' 富集数据则不然。全长 scRNA-seq 数据的常用归一化方法是 TPM 归一化（Li 等, 2009），它源自整体 RNA-seq 分析。

正如可以对细胞计数数据进行归一化以使细胞间数据具有可比性一样，也可以对基因计数进行缩放以改进基因间的比较。基因归一化是将基因计数缩放为零均值和单位方差（z 分数）。这样可以在下游分析中使所有基因具有相等的权重。目前尚无共识是否应对基因进行归一化。Seurat 教程（Butler 等, 2018）通常应用基因缩放，而 Slingshot 方法的作者在其教程中则选择不对基因进行缩放（Street 等, 2018）。这两种选择的偏好围绕是否应在下游分析中对所有基因赋予相同权重，或者基因表达量的大小是否是基因重要性的有用代理。为了尽可能保留数据的生物学信息，本教程选择不对基因进行缩放。

在归一化之后，数据矩阵通常会进行 log(x+1) 转换。此转换有三个重要效果。首先，log 转换后的表达值之间的距离表示对数倍数变化，这是衡量表达变化的标准方法。其次，log 转换减轻了（但未完全消除）单细胞数据中的均值-方差关系（Brennecke 等, 2013）。最后，log 转换减少了数据的偏度，使其更接近下游分析工具通常假设的正态分布。尽管 scRNA-seq 数据实际上并非对数正态分布（Vieth 等, 2017），但这三个效果使 log 转换成为一种粗略但实用的工具。这种实用性在使用 log 转换进行的下游应用（例如差异表达测试（Finak 等, 2015; Ritchie 等, 2015）或批次校正（Johnson 等, 2006; Buttner 等, 2019））中得以体现。然而需要注意的是，归一化数据的 log 转换可能会在数据中引入伪差异表达效应（预印本：Lun, 2018）。当归一化大小因子分布在测试组间显著不同时，这种效应尤为明显。

陷阱与建议：

• 我们推荐使用 scran 对非全长数据集进行归一化。另一个选择是使用 scone 评估归一化方法，尤其适用于基于板的样本。对于全长 scRNA-seq 协议，可以使用整体方法校正基因长度。

• 对基因进行 0 均值和单位方差的缩放尚无共识。我们更倾向于不对基因表达进行缩放。

• 归一化数据应进行 log(x+1) 转换，以便与假设数据呈正态分布的下游分析方法一起使用。

数据校正与整合

上述归一化方法试图去除计数抽样的影响。然而，归一化数据可能仍然包含不期望的变异性。数据校正针对进一步的技术和生物协变量（如批次效应、掉失或细胞周期效应）进行调整。这些协变量并不总是被校正，而应根据下游分析的目标决定是否校正哪些协变量。我们建议将生物和技术协变量的校正分开考虑，因为它们用于不同的目的，且具有不同的挑战。

回归去除生物效应

校正技术协变量可能对揭示潜在生物信号至关重要，而生物协变量校正则有助于突出特定的感兴趣生物信号。最常见的生物数据校正是去除细胞周期对转录组的影响。可以通过对细胞周期得分进行简单线性回归来完成该数据校正，Scanpy 和 Seurat 平台中有实现此功能的工具（Butler 等, 2018; Wolf 等, 2018），或使用包含更复杂混合模型的专用包如 scLVM（Buettner 等, 2015）或 f-scLVM（Buettner 等, 2017）。细胞周期评分所需的标记基因列表可以从文献中获取（Macosko 等, 2015）。这些方法也可以用于回归去除其他已知生物效应，例如线粒体基因表达，后者被解释为细胞应激的指示。

在对数据进行生物效应校正之前，需要考虑几个方面。首先，校正生物协变量并不总是有助于解释 scRNA-seq 数据。去除细胞周期效应可以改善发育轨迹的推断（Buettner 等, 2015; Vento-Tormo 等, 2018），但细胞周期信号本身也具有生物学意义。例如，可以基于细胞周期评分识别增殖细胞群体（参见项目 GitHub 上的案例研究）。另外，生物信号必须在上下文中理解。鉴于生物过程在同一生物体中发生，这些过程之间存在依赖关系。因此，校正一个过程可能会无意中掩盖另一个过程的信号。最后，有人认为细胞大小的变异解释了通常归因于细胞周期的转录组效应（McDavid 等, 2016）。因此，通过归一化或专用工具如 cgCorrect（Blasi 等, 2017）校正细胞大小也部分校正了 scRNA-seq 数据中的细胞周期效应。

回归去除技术效应

用于回归去除生物协变量的回归模型变体也可以应用于技术协变量。单细胞数据中最显著的技术协变量是计数深度和批次。尽管归一化使基因计数在细胞间具有可比性，但计数深度效应通常仍存在于数据中。此计数深度效应可能既是生物的，也可能是技术伪影。例如，细胞大小可能不同，从而导致 mRNA 分子计数不同。然而，技术计数效应可能在归一化后仍然存在，因为没有任何缩放方法能够推测由于采样不足而未检测到的基因的表达值。回归去除计数深度效应可以改善依赖于发现细胞间过渡的轨迹推断算法的性能（参见项目 GitHub 上的案例研究）。在校正多个协变量（例如细胞周期和计数深度）时，应一次性对所有协变量进行回归，以考虑协变量之间的依赖关系。

另一种去除计数效应的替代策略是使用更严格的归一化程序，如降采样或非线性归一化方法（参见“归一化”部分）。这些方法在基于板的 scRNA-seq 数据集中可能尤为相关，因为细胞的计数深度变化较大，可能掩盖细胞间的异质性。

批次效应和数据整合

当细胞以不同组处理时，可能会出现批次效应。这些组可以由不同芯片上的细胞、不同测序通道中的细胞或在不同时间点采集的细胞组成。细胞经历的不同环境可能会影响转录组的测量，甚至转录组本身。批次效应的影响存在于多个层面：在实验中的细胞组之间、在同一实验室进行的不同实验之间，或在不同实验室的数据集之间。我们将前者与后两者区分开来。校正同一实验中样本或细胞之间的批次效应是批次校正的经典情境，源于整体 RNA-seq 分析。我们将其与多实验数据的整合（即数据整合）区分开来。批次效应通常使用线性方法校正，而数据整合则使用非线性方法。

最近对经典批次校正方法的比较显示，ComBat（Johnson 等, 2006）在低到中等复杂度的单细胞实验中表现良好（Buttner 等, 2019）。ComBat 基于基因表达的线性模型，将批次对数据均值和方差的影响均考虑在内（图3）。无论采用何种计算方法，最有效的批次校正方法是通过巧妙的实验设计预先消除批次效应（Hicks 等, 2017）。可以通过在实验条件和样本间混合细胞来避免批次效应。通过细胞标记（预印本：Gehring 等, 2018）或遗传变异（Kang 等, 2018）等策略，可以对实验中混合的细胞进行去混合。

图3. 批次校正前后 UMAP 可视化

细胞按来源样本着色。在批次校正前，批次间的分离非常明显，而在批次校正后则不明显。此处使用 ComBat 对来自 Haber 等（2017）的鼠肠上皮数据进行了批次校正。

相比批次校正，数据整合方法面临的额外挑战在于数据集之间的组成差异。在估计批次效应时，ComBat 使用批次内的所有细胞来拟合批次参数。这种方法会将批次效应与不同数据集中细胞类型或状态的生物学差异混淆。为克服这一问题，已开发出多种数据整合方法，例如典型相关分析（CCA；Butler 等, 2018）、互邻（MNN；Haghverdi 等, 2018）、Scanorama（预印本：Hie 等, 2018）、RISC（预印本：Liu 等, 2018）、scGen（预印本：Lotfollahi 等, 2018）、LIGER（预印本：Welch 等, 2018）、BBKNN（预印本：Park 等, 2018）和 Harmony（预印本：Korsunsky 等, 2018）。尽管数据整合方法也适用于简单的批次校正问题，但由于非线性数据整合方法自由度增加，需注意避免过度校正。例如，MNN 在更简单的批次校正情境中被证明不如 ComBat（Buttner 等, 2019）。需要进一步的比较研究来评估这些方法的普遍适用性。

表达恢复

另一类技术数据校正是表达恢复（也称去噪或插补）。单细胞转录组测量包含多种噪声源（Grün 等, 2014；Kharchenko 等, 2014；Hicks 等, 2017）。掉失是这种噪声中一个显著方面。推断掉失事件、用适当的表达值替换这些零值并减少数据集中的噪声是一些近期工具的目标（MAGIC: van Dijk 等, 2018；DCA: Eraslan 等, 2019；scVI: Lopez 等, 2018；SAVER: Huang 等, 2018；scImpute: Li & Li, 2018）。研究表明，进行表达恢复有助于提高基因-基因相关性的估计精度（van Dijk 等, 2018；Eraslan 等, 2019）。此外，此步骤可与归一化、批次校正和其他下游分析集成，如在 scVI 工具中实现的（Lopez 等, 2018）。虽然大多数数据校正方法以归一化数据为输入，但某些表达恢复方法基于负二项噪声分布，因此在原始计数数据上运行。在应用表达恢复时应注意，没有一种方法是完美的。因此，任何方法都可能在数据中过度或不足地校正噪声。事实上，表达恢复后已报道出现错误的相关信号（Andrews & Hemberg, 2018）。由于在实际应用中评估成功的表达恢复非常困难，这一情境对用户是否进行数据去噪提出了挑战。此外，目前可用的表达恢复方法在处理大数据集方面仍存在问题。目前尚无共识关于如何使用去噪数据（参见“处理数据阶段”部分）。一种谨慎的方法是仅将表达恢复用于数据的可视化展示，而非在探索性数据分析中生成假设。彻底的实验验证在此尤为重要。

陷阱与建议：

• 仅在进行轨迹推断时回归去除生物协变量，且需确保所去除的生物协变量不会掩盖其他感兴趣的生物过程。

• 共同回归去除技术和生物协变量，而非逐步进行。

• 基于板的数据集预处理可能需要回归去除计数、使用非线性归一化方法或降采样进行归一化。

• 当批次之间的细胞类型和状态组成一致时，建议使用 ComBat 进行批次校正。

• 数据整合和批次校正应采用不同方法。数据整合工具可能会对简单的批次效应过度校正。

• 用户应谨慎对待仅在表达恢复后出现的信号，探索性分析在没有此步骤的情况下可能表现更佳。

特征选择、降维和可视化

一个人类单细胞 RNA-seq 数据集可以包含多达 25,000 个基因的表达值。许多基因对于特定的 scRNA-seq 数据集没有信息性，且许多基因大多包含零计数。即使在 QC 步骤中过滤掉这些零计数基因后，单细胞数据集的特征空间仍可能超过 15,000 维。为了减轻下游分析工具的计算负担、减少数据噪声并可视化数据，可以使用多种方法来降低数据集的维度。

特征选择

减少 scRNA-seq 数据集维度的第一步通常是特征选择。在此步骤中，数据集被过滤，以保留那些对数据变异性“有信息性”的基因。因此，通常使用高变异基因（HVGs）（Brennecke 等, 2013）。根据任务和数据集的复杂性，通常选择 1,000 到 5,000 个 HVGs 进行下游分析（参见图 EV1 和数据集 EV1）。Klein 等人（2015）的一项初步研究结果表明，下游分析对 HVG 数量的具体选择具有鲁棒性。作者在将 HVG 数量从 200 调整到 2,400 之间时，在 PCA 空间中报告了类似的低维表示。基于此结果，我们倾向于选择更多的 HVGs 以确保信息量的丰富。

图 EV1. 用于不同规模数据集的高变异基因 (HVGs) 数量

数据通过对近期 scRNA-seq 分析论文的简要手动调查获得。绘制的数据，以及关于 scRNA-seq 技术、发表年份、参考文献和每个细胞的读取数的更多信息，可在数据集 EV1 中查阅。

一种简单但流行的 HVG 选择方法在 Scanpy 和 Seurat 中均有实现。该方法中，基因按其平均表达分箱，并在每个箱中选择方差与均值比最高的基因作为 HVGs。此算法的不同变体可处理计数数据（Seurat）或对数转换数据（Cell Ranger）。理想情况下，HVGs 应在技术数据校正后选择，以避免因批次效应等导致的高变异性基因被选中。其他 HVG 选择方法可参见 Yip 等（2018）的综述。

降维

在特征选择后，可以通过专用降维算法进一步减少单细胞表达矩阵的维度。这些算法将表达矩阵嵌入到低维空间中，旨在用尽可能少的维度捕捉数据的潜在结构。此方法适用于单细胞 RNA-seq 数据，因为这类数据本质上是低维的（Heimberg 等, 2016）。换句话说，细胞表达谱所在的生物学流形可以用远少于基因数的维度充分描述。降维的目标是找到这些维度。

降维方法的主要目标有两个：可视化和概括。可视化是尝试用二维或三维来最佳地描述数据集。这些降维结果用于散点图中的坐标，以获得数据的视觉表示。概括则不限定输出组件的数量，而是通过找到数据的内在维度来提取其主要成分，这对下游分析有帮助。虽然二维可视化输出不应用于总结数据集，但可使用概括方法通过主成分来可视化数据。然而，专用的可视化技术通常会提供更好的变异性表示。

降维是通过特征空间维度（基因表达向量）的线性或非线性组合生成的，尤其在非线性情况下，降维的解释性在此过程中被牺牲。一些常用降维方法的示例应用见图 4。随着可选方法的增加，本教程不详细审查这些方法，而是简要概述了可能帮助用户选择常见降维方法的实际考虑因素。关于单细胞分析降维的更详细综述可参见 Moon 等（2018）。

图4. scRNA-seq 数据的常用可视化方法

来自 Haber 等人（2017）的鼠肠上皮区域数据在以下首两个成分上可视化：(A) PCA，(B) t-SNE，(C) 扩散图，(D) UMAP 和 (E) 基于 ForceAtlas2 的力导向图布局。细胞按计数深度着色。(F) 前 31 个主成分（PC）的方差解释率图。该图的“拐点”位于 PC5 到 PC7 之间，用于选择分析数据集的相关 PC。

两种常用的降维技术（主要用于数据概括）是主成分分析（PCA；Pearson, 1901）和扩散图（Coifman 等, 2005），Haghverdi 等人（2015）使其在单细胞分析中流行。PCA 是一种线性方法，通过在每个维度中最大化残留方差来生成降维结果。尽管 PCA 在低维度上捕捉数据结构的能力不及非线性方法，但它是许多聚类或轨迹推断分析工具的基础。PCA 通常用作非线性降维方法的预处理步骤。一般情况下，PCA 通过其前 N 个主成分对数据集进行概括，N 的选择可以通过“拐点”启发法（见图 4F）或基于置换测试的 Jackstraw 方法来确定（Chung & Storey, 2015；Macosko 等, 2015）。PCA 的简单线性特性使得在降维空间中各区域的距离具有一致的解释性。因此，我们可以将感兴趣的量与主成分相关联以评估其重要性。例如，可以将主成分投影到技术协变量上，以检查 QC、数据校正和归一化步骤的效果（Buttner 等, 2019），或展示数据集中基因的重要性（Chung & Storey, 2015）。

扩散图是一种非线性数据概括技术。由于扩散成分强调数据中的过渡，因此在研究分化等连续过程时尤为适用。通常，每个扩散成分（即扩散图维度）突出不同细胞群体的异质性。

可视化

出于可视化目的，通常使用非线性降维方法（图4）。最常用的 scRNA-seq 可视化降维方法是 t 分布随机邻域嵌入（t-SNE；van der Maaten & Hinton, 2008）。t-SNE 主要关注捕捉局部相似性，代价是牺牲全局结构。因此，这些可视化方法可能会夸大细胞群体之间的差异，并忽略这些群体之间的潜在联系。此外，t-SNE 的困惑度参数选择也较为复杂，不同的取值会显著影响聚类数目（Wattenberg 等, 2016）。t-SNE 的常见替代方法包括一致性近似投影（UMAP；预印本：McInnes & Healy, 2018）或基于图的工具如 SPRING（Weinreb 等, 2018）。UMAP 和 SPRING 的力导向布局算法 ForceAtlas2 被认为能较好地呈现潜在拓扑结构（Wolf 等, 2019，补充注释4）。UMAP 在速度和扩展性方面表现出色，适合处理大量细胞数据（Becht 等, 2018）。因此，在缺乏特定生物学问题时，我们认为 UMAP 是探索性数据可视化的最佳选择。此外，UMAP 还可以在超过二维的情况下对数据进行概括。尽管我们尚未见到将 UMAP 用于数据概括的应用，但它可能成为 PCA 的替代方案。

一种替代经典细胞层级可视化的方法是基于分区的图形抽象（PAGA；Wolf 等, 2019）。此工具在使用聚类对可视化进行粗粒化的同时，能适当近似数据拓扑。结合上述任何一种可视化方法，PAGA 可生成粗粒化可视化，这在分析大量细胞的单细胞数据时特别有助于解释。

陷阱与建议：

• 根据数据集的复杂性选择 1,000 到 5,000 个高变异基因。

• 当基因表达值已归一化至零均值和单位方差，或使用模型拟合残差作为归一化表达值时，不适用基于基因表达均值和方差的特征选择方法。因此，必须在选择 HVG 之前确定需要执行的预处理。

• 应将降维方法分别考虑用于数据概括和可视化。

• 我们推荐使用 UMAP 进行探索性可视化；PCA 用于通用数据概括；扩散图作为轨迹推断概括的替代 PCA 方法。

• PAGA 与 UMAP 的组合适合可视化特别复杂的数据集。

预处理数据的阶段

尽管我们在上述内容中将 scRNA-seq 常用预处理步骤描述为一个连续流程，下游分析通常偏向使用不同级别的预处理数据，建议根据下游应用调整预处理。为帮助新用户理解此情况，我们将预处理划分为五个数据处理阶段：(i) 原始数据，(ii) 归一化数据，(iii) 校正数据，(iv) 特征选择后的数据，和 (v) 降维数据。这些数据处理阶段分为三层预处理层：测量数据、校正数据和降维数据。细胞和基因 QC 是必须执行的步骤，因此在此分类中被省略。尽管这些层的顺序表示了典型的 scRNA-seq 分析工作流程，也可以跳过某些层或对处理阶段顺序进行轻微调整。例如，对于单批次数据集，数据校正可能不必要。表 1 总结了每层预处理数据的适用下游应用。

Table 1. Stages of data processing and appropriate downstream applications

Pre‐processing layer	Stage of data processing	Appropriate applications
Measured	1) Raw	Statistical testing (Differential expression: marker genes, genes over condition, genes over time)
2) Normalized (+ log transformed)
Corrected	3.1) Corrected (technical correction)	Visual comparison of data (plotting)
3.2) Corrected (biological correction)	Pre‐processing for trajectory inference
Reduced	4) Feature selected	Visualization, trajectory inference
5) Dimensionality reduced (summarized)	Visualization, clustering, KNN graph inference, trajectory inference

表 1 中的预处理阶段分为三组：测量数据、校正数据和降维数据。我们将测量数据定义为原始数据及保留零值结构的处理数据。通过使用细胞特定的因子缩放计数数据，全局缩放归一化方法在进行 log(x+1) 转换后仍保留零表达值。相反，去除不期望的变异性会替换零表达值。校正数据层代表数据的“最干净”版本，是对潜在生物信号的最佳近似。我们称最终预处理层为降维数据。此数据层强调数据的主导特征，使用较少的特征集进行描述。

上述特征决定了预处理数据在特定下游应用中的适用性。作为最后的预处理阶段，降维数据是广泛适用的数据层的自然候选。然而，差异表达的测试在基因空间中才具有生物学解释性，而降维数据未（完全）呈现基因空间。降维数据的优势在于概括生物学信息并减少噪声，从而避免掩盖生物信号。因此，降维数据用于需要数据概括的探索性方法（如可视化、邻域图推断、聚类）以及计算复杂的下游分析工具（如轨迹推断）。实际上，许多轨迹推断方法在其工具中已经集成了降维。

个别基因的表达谱只能在基因空间中比较，这一空间存在于测量数据和校正数据中。表达谱的比较可以通过可视化和统计进行。我们主张将可视化和统计比较分配在不同的数据层上。对于基因表达的可视化检查，校正数据最为合适。如果原始数据用于可视化比较，用户需固有地理解数据中的偏差以正确解读结果。校正数据有助于这种解释。但在这里应分别考虑技术和生物协变量的校正。虽然校正生物协变量可能增强特定生物信号，但会导致对基础生物学的描述不准确，并掩盖其他可能相关的信号。因此，生物校正数据主要适用于专注于特定生物过程的分析工具，例如轨迹推断方法。

基因表达的统计比较最适合在测量数据层上进行。没有完美的数据校正方法可以实现去噪、批次校正或其他变异源的校正。因此，数据校正方法不可避免地会对数据进行过度或不足校正，从而以非预期方式改变至少一些基因表达谱的方差。基因表达的统计检验依赖于评估背景方差作为数据噪声的零模型。由于数据校正通常减少背景变异（见图 EV2），背景变异被过度校正的基因更可能被评估为显著差异表达。此外，某些数据校正方法（如 ComBat）将不符合实验设计的表达信号解释为噪声，并将其从数据中去除。除了低估噪声外，对实验设计信号的优化还可能导致效应量的高估。鉴于这些考虑，使用测量数据而非校正数据作为输入，构成了更为保守的差异测试方法。在使用测量数据时，应在差异测试模型中考虑技术协变量。

图 EV2. 基因表达数据在批次校正和去噪后的变异系数（CoV）变化

负值表示数据校正后 CoV 的降低。上排显示了 ComBat 批次校正后 (A) 鼠肠上皮 (mIE) 和 (B) 鼠胚胎干细胞 (mESC) 数据的 CoV 变化。下排显示了 DCA 去噪后 (C) mIE 和 (D) mESC 数据的 CoV 变化。mIE 数据来自 Haber 等（2017），mESC 数据来自 Klein 等（2015）。

上述观点得到一项近期关于 scRNA-seq 差异表达方法的比较支持，该研究仅使用原始数据和归一化数据作为输入（Soneson & Robinson, 2018）。该研究中使用的归一化数据仅涉及全局缩放方法。然而，目前许多非线性归一化方法模糊了归一化和数据校正的界限（参见“归一化”部分）。此类归一化数据可能不再适合用于差异测试输入。

陷阱与建议：

• 使用测量数据进行统计测试，使用校正数据进行数据的可视化比较，并使用降维数据进行基于基础生物数据流形发现的其他下游分析。

下游分析

在预处理之后，使用我们称为下游分析的方法来提取生物学见解并描述潜在的生物系统。这些描述是通过拟合可解释模型到数据中获得的。此类模型的示例包括：具有相似基因表达谱的细胞群代表细胞类型簇；相似细胞间基因表达的微小变化表示连续的（分化）轨迹；或具有相关表达谱的基因表示共调控。

下游分析可分为细胞层级和基因层级方法，如图 5 所示。细胞层级分析通常聚焦于两个结构的描述：簇和轨迹。这些结构进而可以在细胞层级和基因层级进行分析，形成簇分析和轨迹分析方法。总体而言，簇分析方法试图基于细胞分组来解释数据的异质性。相比之下，在轨迹分析中，数据被视为动态过程的快照，轨迹分析方法探讨这一潜在过程。

图 5. 下游分析方法概览

方法分为细胞层级和基因层级分析。细胞层级分析方法进一步细分为簇分析和轨迹分析分支，这些分支还包括基因层级分析方法。所有带蓝色背景的方法均为基因层级方法。

在此，我们描述细胞层级和基因层级的簇和轨迹分析工具，之后再详细介绍独立于这些细胞结构的基因层级分析。

簇分析

聚类

将细胞组织成簇通常是任何单细胞分析的第一个中间结果。簇允许我们推断成员细胞的身份。通过基于基因表达谱的相似性对细胞进行分组得到簇。表达谱相似性通过距离度量确定，通常以降维表示作为输入。例如，一种常见的相似性评分方法是在主成分分析（PCA）降维后的表达空间中计算欧几里得距离。生成细胞簇的方法有两种：聚类算法和社区检测方法。

聚类是一种经典的无监督机器学习问题，直接基于距离矩阵。通过最小化簇内距离或在降维表达空间中找到密集区域将细胞分配到簇。流行的 k-means 聚类算法通过确定簇质心并将细胞分配到最近的簇质心，将细胞划分为 k 个簇。质心位置经过迭代优化（MacQueen, 1967）。此方法需要输入预期簇数，通常未知且需启发式校准。k-means 应用于单细胞数据的距离度量方法多样，替代标准欧几里得距离的方法包括余弦相似性（Haghverdi 等, 2018）、基于相关性的距离度量（Kim 等, 2018）或 SIMLR 方法，该方法为每个数据集使用高斯核学习距离度量（Wang 等, 2017）。一项最新比较研究表明，基于相关性的距离在与 k-means 结合或作为高斯核的基础时可能优于其他距离度量（Kim 等, 2018）。

社区检测方法是图分割算法，依赖于单细胞数据的图表示。此图表示使用 K-最近邻（KNN）方法获得。细胞在图中表示为节点，每个细胞与 K 个最相似的细胞相连，通常在 PCA 降维表达空间中使用欧几里得距离获得 K 值。根据数据集的大小，K 通常设置在 5 到 100 个最近邻之间。结果图捕捉了表达数据的潜在拓扑结构（Wolf 等, 2019）。表达空间中密集采样的区域在图中表示为密集连接区域，这些密集区域通过社区检测方法检测。与聚类相比，社区检测更快，因为仅需考虑相邻细胞对属于同一簇。因此，此方法大大缩小了可能簇的搜索空间。

自 PhenoGraph 方法（Levine 等, 2015）以来，单细胞数据集的标准聚类方法已成为在单细胞 KNN 图上使用 Louvain 算法实现的多分辨率模块化优化（Newman & Girvan, 2004；Reichardt & Bornholdt, 2006；Blondel 等, 2008）。此方法是 Scanpy 和 Seurat 单细胞分析平台中实现的默认聚类方法，被证明在单细胞 RNA-seq 数据（Duò 等, 2018；Freytag 等, 2018）、流式和质谱细胞计量数据（Weber & Robinson, 2016）上表现优越。Louvain 算法概念上通过检测细胞之间的连接数量超过预期的细胞组来识别社区。优化的模块化函数包含一个分辨率参数，允许用户确定簇分区的规模。通过对 KNN 图进行子集化，还可以对特定簇进行二次聚类，从而允许用户识别细胞类型簇内的细胞状态（Wagner 等, 2016），但也可能导致仅由数据噪声产生的模式。

陷阱与建议：

• 我们推荐在单细胞 KNN 图上使用 Louvain 社区检测进行聚类。

• 聚类不必在单一分辨率上进行，特定细胞簇的二次聚类是关注数据集中更细致子结构的有效方法。

簇注释

在基因层级上，通过找到每个簇的基因标记对聚类数据进行分析。这些所谓的标记基因表征了簇，并用于为其注释一个有意义的生物学标签，代表簇内细胞的身份。由于任何聚类算法都会生成数据的分区，所识别簇的有效性只能通过成功注释所表示的生物学来确定。

虽然可能倾向于假设单细胞数据中的簇代表细胞类型，但决定细胞身份的变量有多个轴（Wagner 等, 2016；Clevers 等, 2017）。首先，什么构成细胞类型并不总是清晰。例如，“T 细胞”可能对某些人来说是一个满意的细胞类型标签，而另一些人可能会在数据集中寻找 T 细胞亚型并区分 CD4+ 和 CD8+ T 细胞（Wagner 等, 2016；Clevers 等, 2017）。此外，不同状态的相同细胞类型可能在不同簇中检测到。基于上述原因，最好使用“细胞身份”而不是“细胞类型”一词。在聚类和注释簇之前，用户必须决定感兴趣的注释细节级别，从而决定聚类分辨率。

识别和注释簇依赖于使用描述个体细胞身份预期表达谱的外部信息源。由于近期和持续的努力（例如小鼠大脑图谱（Zeisel 等, 2018）或人类细胞图谱（Regev 等, 2017）），参考数据库日益可用，大大简化了细胞身份的注释。在缺少相关参考数据库时，可以通过将数据派生的标记基因与文献标记基因比较来注释细胞身份（见项目 GitHub 上的案例研究），或直接可视化文献标记基因的表达值（图 6B）。需注意的是，后一种方法将用户限制在源自整体表达研究的传统细胞类型理解中，而非细胞身份。此外，有研究表明，常用的细胞表面标记在定义细胞身份方面存在局限性（Tabula Muris Consortium 等, 2018）。

图 6. 来自 Haber 等人（2017）的鼠肠上皮数据集的簇分析结果

(A) 使用 Louvain 聚类的细胞身份簇在 UMAP 表示中可视化。(B) 通过细胞身份标记表达识别干细胞 (Slc12a2)、肠细胞 (Arg2)、杯状细胞 (Tff3) 和潘氏细胞 (Defa24)。校正后的表达水平从低（灰色）到高（红色）显示。标记基因也可能在其他细胞身份群体中表达，例如杯状细胞和潘氏细胞。(C) 近端（上）和远端（下）肠上皮区域的细胞身份组成热图。深红色表示较高的相对细胞密度。

可以通过两种方式使用参考数据库信息来注释簇：使用数据派生的标记基因或完整的基因表达谱。可以通过在簇内的细胞与数据集中其他所有细胞之间进行差异表达（DE）测试来找到标记基因集合（参见“差异表达测试”）。通常关注在感兴趣簇中上调的基因。由于标记基因预计具有显著的差异表达效果，通常使用简单的统计测试（如 Wilcoxon 秩和检验或 t 检验）按表达差异对基因进行排序。各统计检验中排名靠前的基因被视为标记基因。簇可以通过将数据集中的标记基因与参考数据集中的标记基因进行富集测试、Jaccard 指数或其他重叠统计来注释。参考网站如 www.mousebrain.org（Zeisel 等, 2018）或 http://dropviz.org/（Saunders 等, 2018）允许用户在参考数据集中可视化数据集标记基因的表达，方便细胞身份注释。

在检测标记基因时需要注意两个方面。首先，标记基因的 P 值基于假设获得的细胞簇代表生物学真实情况。如果簇分配存在不确定性，则在统计测试中必须考虑簇分配与标记基因检测之间的关系。此关系的产生是因为簇和标记基因通常基于相同的基因表达数据确定。DE 测试中的原假设是假设基因在两组间的表达值分布相同。然而，由于在标记基因检测中两组是由聚类方法定义的，它们的基因表达谱在设计上存在差异。因此，即使在使用 splatter 生成的随机数据进行聚类时，我们也会找到显著的标记基因（Zappia 等, 2017）（参见附录补充文本 S3）。可以使用置换检验来考虑聚类步骤，从而获得适当的显著性衡量方法，相关内容在附录补充文本 S3 中详述。最近的差异表达工具也专门解决了这一问题（预印本：Zhang 等, 2018）。在当前设置中，P 值往往被放大，可能导致标记基因数量的高估。然而，基于 P 值的基因排序不受影响。假设聚类具有生物学意义，排名靠前的标记基因仍然是最佳的标记基因候选。在初始情况下，可以通过视觉检查来粗略验证标记基因。需要强调的是，P 值膨胀特定于通过无监督聚类方法定义的细胞身份簇。当通过单个基因的表达确定细胞身份簇时，P 值可以按预期解释用于其他所有基因。尽管这种单变量方法在特定情况下常见（例如胰岛素在 β 细胞中的表达或血红蛋白在红细胞中的表达），但通常不推荐使用。

第二，标记基因用于区分数据集中的簇，因此不仅依赖于细胞簇，也依赖于数据集的组成。如果数据集的组成未准确代表背景基因表达，则检测到的标记基因会偏向于缺失的成分。在低细胞多样性数据集中计算标记基因时，尤其需要考虑此方面。

最近，自动簇注释方法已出现。通过直接比较注释的参考簇的基因表达谱与单个细胞，工具如 scmap（Kiselev 等, 2018b）或 Garnett（预印本：Pliner 等, 2019）可以在参考数据集与待分析数据集之间转移注释。因此，这些方法可以同时执行注释和簇分配，而无需数据驱动的聚类步骤。由于实验条件下细胞类型和状态组成可能不同（Segerstolpe 等, 2016；Tanay & Regev, 2017），基于参考数据的聚类不应替代数据驱动方法。

聚类、簇注释、二次聚类和再注释的迭代过程可能耗时。自动簇注释方法极大地加快了这一过程。然而，自动和手动方法各有利弊，难以推荐其一优于另一种。速度提升通常伴随灵活性上的妥协。如上所述，参考图谱并不一定包含与待分析数据集完全相同的细胞身份。因此，不能放弃手动注释的标记基因计算。对于包含多个簇的大型数据集，当前的最佳实践是结合这两种方法。为了加快进程，可以使用自动细胞身份注释对细胞进行粗略标记，并确定需要二次聚类的区域。随后，计算数据集簇的标记基因并与参考数据集或文献中的已知标记基因集合进行比较。对于较小的数据集和缺乏参考图谱的数据集，手动注释足以完成任务。

陷阱与建议：

• 不要使用标记基因 P 值来验证细胞身份簇，尤其是当检测到的标记基因未能帮助注释簇时。P 值可能会被放大。

• 注意相同细胞身份簇的标记基因在不同数据集中可能有所不同，仅仅由于数据集细胞类型和状态组成不同。

• 如果存在相关的参考图谱，建议使用自动簇注释与基于数据派生标记基因的手动注释相结合来注释簇。

组成分析

在细胞层级上，可以从其组成结构的角度分析聚类数据。组成数据分析关注落入每个细胞身份簇的细胞比例。这些比例可能会因疾病而发生变化。例如，研究显示沙门氏菌感染会增加小鼠肠上皮中肠细胞的比例（Haber 等, 2017）。

研究单细胞数据中的组成变化需要足够的细胞数量以稳健地评估细胞身份簇的比例，以及足够的样本数量以评估细胞身份簇组成的预期背景变异性。由于合适的数据集最近才变得可用，尚未开发出专门的工具。在上述小鼠研究中，细胞身份计数使用泊松过程建模，将实验条件作为协变量，总检测到的细胞数作为偏移量。此时可以通过回归系数进行统计检验，以评估某种特定细胞身份的频率是否显著改变。然而，同一数据集中其他细胞身份的检验彼此并不独立。如果一个细胞身份簇的比例发生变化，则其他簇的比例也必须发生变化。因此，不能通过此模型评估整体组成是否显著改变。在缺少专门工具的情况下，组成数据的视觉比较可以为样本间组成变化提供信息（图 6C）。未来在该领域的发展可能会借鉴质谱细胞计量学（例如 Tibshirani 等, 2002；Arvaniti & Claassen, 2017；Lun 等, 2017；Weber 等, 2018）或微生物组文献（Gloor 等, 2017），其中组成数据分析已受到更多关注。

陷阱与建议：

• 请考虑样本间细胞身份簇比例变化的统计检验是彼此相关的。

轨迹分析

轨迹推断

细胞多样性不能仅通过聚类等离散分类系统充分描述。驱动观察到的异质性发展的生物学过程是连续的过程（Tanay & Regev, 2017）。因此，为了捕捉细胞身份之间的过渡、分化分支过程或生物功能的渐进、不同步变化，我们需要基因表达的动态模型。这类方法称为轨迹推断（TI）。

轨迹推断方法将单细胞数据视为连续过程的快照。此过程通过在细胞空间中寻找路径来重构，这些路径最小化相邻细胞之间的转录变化（图 7A 和 B）。沿着这些路径排列的细胞由伪时间变量描述。尽管该变量与从起始细胞的转录距离相关，但通常被解释为发育时间的近似（Moignard 等, 2015；Haghverdi 等, 2016；Fischer 等, 2018；Griffiths 等, 2018）。

图 7. Haber 等人（2017）鼠肠上皮数据的轨迹分析与图抽象

(A) 使用 Slingshot 推断的近端和远端肠细胞分化轨迹，远端谱系按伪时间从红色到蓝色显示。数据集中的其他细胞为灰色。(B) PCA 空间中的簇间 Slingshot 轨迹。簇的缩写如下：EP—肠细胞前体；Imm. Ent.—未成熟肠细胞；Mat. Ent.—成熟肠细胞；Prox.—近端；Dist.—远端。(C) 图 7A 中远端肠细胞轨迹的伪时间密度，不同颜色表示伪时间区间中的主要簇标签。(D) 数据集的抽象图表示投影到 UMAP 表示中。簇显示为带颜色的节点，出现在其他轨迹中的簇标注以便对比。“TA”表示过渡增殖细胞。(E) 使用“GAM” R 库展示在一般肠细胞轨迹中的基因表达动态。

自 Monocle（Trapnell 等, 2014）和 Wanderlust（Bendall 等, 2014）确立了轨迹推断（TI）领域以来，可用方法数量急剧增加。当前的 TI 方法在所建模路径的复杂性上各不相同，模型包括简单的线性或分叉轨迹，到复杂的图、树或多分叉轨迹。在 Saelens 等人（2018）的全面比较中，得出的结论是没有单一方法可以在所有类型的轨迹上表现最佳。相反，应根据预期轨迹的复杂性选择 TI 方法。比较表明，Slingshot（Street 等, 2018）在从线性到多分叉模型的简单轨迹上优于其他方法。如果预期更复杂的轨迹，作者推荐使用 PAGA（Wolf 等, 2019）。如果已知确切的轨迹模型，也可使用更为专门的方法以提高性能（Saelens 等, 2018）。通常，任何推断的轨迹都应使用替代方法加以确认，以避免方法偏差。

在典型工作流程中，当内置降维步骤时，TI 方法通常应用于降维数据或校正数据。由于细胞内通常同时发生多种生物学过程，因此可以回归去除其他过程的生物效应，以隔离预期轨迹。例如，T 细胞在成熟过程中可能经历细胞周期转换（Buettner 等, 2015）。此外，由于多种高性能 TI 方法依赖聚类数据，TI 通常在聚类之后进行。推断的轨迹中的簇可能代表稳定或亚稳定状态（见“亚稳定状态”部分；图 7B 和 C）。随后，可以将 RNA 速度叠加在轨迹上以增加方向性（La Manno 等, 2018）。

推断的轨迹不一定表示生物学过程。最初，这些轨迹仅表示转录相似性。很少有 TI 方法评估其模型中的不确定性（Griffiths 等, 2018）。因此，需要进一步的信息来验证是否确实捕捉到生物学过程，这些信息可以包括扰动实验、推断的调控基因动态及 RNA 速度的支持。

陷阱与建议：

• 我们推荐使用 Saelens 等人（2018）的综述作为指南。

• 推断的轨迹不一定代表生物学过程，建议收集更多证据以解释轨迹。

基因表达动态

验证推断的轨迹不是转录噪声的结果的一种方法是分析轨迹的基因层级。伪时间中平滑变化的基因表征了轨迹，可用于识别潜在的生物过程。此外，轨迹相关基因预计包含调控所建模过程的基因。调控基因帮助我们理解生物过程的触发机制，并代表潜在的药物靶点（Gashaw 等, 2012）。

早期寻找轨迹相关基因的方法涉及沿轨迹的细胞簇间差异表达（DE）测试（Haghverdi 等, 2016；Alpert 等, 2018），而现在则通过将基因表达与伪时间回归来检测沿轨迹变化的基因。为了强制表达沿此协变量平滑变化，伪时间通过拟合样条或额外的局部回归步骤（如 loess）进行平滑。回归框架在其噪声模型假设和用于描述表达随伪时间变化的函数类上有所不同。通过为基因对伪时间的依赖性进行模型选择可以得到潜在的调控基因。此伪时间上的 DE 检验与簇间 DE 测试同样受到轨迹推断方法的影响（见“簇注释”部分）。因此，在此设置中获得的 P 值不应视为显著性评价。

目前专门的基因时间动态工具较少。BEAM 是集成在 Monocle TI 流程中的工具（Qiu 等, 2017a），允许检测特定分支的基因动态。除此流程外，用户可以选择 LineagePulse（https://github.com/YosefLab/LineagePulse），该工具考虑掉失噪声，但仍在开发中，或使用 limma 包（Ritchie 等, 2015）或标准 R 库编写自己的测试框架。在线 Slingshot 教程（Street 等, 2018）和图 7E 中提供了此类示例。

鉴于目前可用的工具有限，尚无法确定研究基因时间动态的最佳实践。所有上述方法均可用于探索性地研究基因动态。未来，高斯过程可能为研究基因时间动态提供自然模型。此外，检测调控模块而非单个基因可能会提高信噪比并促进生物学解释。

亚稳定状态

细胞层级的轨迹分析通过伪时间研究细胞密度。假设细胞是无偏采样的，轨迹中的密集区域表明偏好转录状态。当将轨迹解释为时间过程时，这些密集区域可能代表亚稳定状态，例如在发育过程中（Haghverdi 等, 2016）。可以通过绘制伪时间坐标的直方图找到这些亚稳定状态（图 7C）。

细胞层级分析的统一

聚类和轨迹推断代表了单细胞数据的两种不同视角。这两种视角可以在粗粒化的图表示中统一。通过将单细胞簇表示为节点，并将簇之间的轨迹表示为边，可以同时表示数据的静态和动态特性。这种统一方法由基于分区的图抽象工具 PAGA 提出（图 7D；Wolf 等, 2019）。PAGA 使用一个统计模型来表示细胞簇之间的相互作用，并在细胞相似度高于预期的簇节点之间建立边。在最近的综述中，PAGA 被认为优于其他轨迹推断方法（Saelens 等, 2018）。这是唯一能够处理断开拓扑结构和包含循环的复杂图的方法，这使得 PAGA 成为可视化整个数据集拓扑的有效工具，尤其适合探索性分析。

基因层级分析

虽然我们迄今为止专注于描述细胞结构的基因层级分析方法，但单细胞数据的基因层级分析具有更广泛的应用范围。差异表达测试、基因集合分析和基因调控网络推断直接研究数据中的分子信号。与描述细胞异质性不同，这些方法将异质性作为理解基因表达的背景。

差异表达测试

在表达数据中常见的问题是确定是否有基因在两个实验条件之间表现出差异表达。差异表达（DE）测试是一个源于整体基因表达分析的成熟问题（Scholtens & von Heydebreck, 2005）。相较于整体差异测试的优势在于，我们可以在单细胞背景下通过在细胞身份簇内进行测试来考虑细胞异质性，从而了解特定细胞身份在特定实验条件下的转录反应（Kang 等, 2018）。

尽管两者旨在回答相同的问题，整体和单细胞 DE 工具在方法上存在差异。整体方法旨在从少量样本中准确估计基因方差，而单细胞数据不存在此问题。另一方面，单细胞数据包含独特的技术噪声伪影，如掉失和高细胞间变异性（Hicks 等, 2017；Vallejos 等, 2017）。这些伪影在专门为单细胞数据设计的方法中得以考虑（Kharchenko 等, 2014；Finak 等, 2015）。然而，近期一项大规模比较研究显示，整体 DE 测试包的表现与表现最佳的单细胞工具相当（Soneson & Robinson, 2018）。此外，当通过引入基因权重来将整体工具适应于单细胞数据时，这些工具据称优于其单细胞对应物（Van den Berge 等, 2018）。根据该比较，表现最佳的 DE 分析工具是结合 ZINB-wave 估计权重的 DESeq2（Love 等, 2014）和 EdgeR（Robinson 等, 2010）。需要独立的比较研究来确认包含加权整体 DE 测试方法的结果。

加权整体 DE 测试的性能提升以计算效率为代价。鉴于单细胞实验中细胞数量日益增加，算法运行时间成为选择方法时日益重要的考量。因此，单细胞工具 MAST（Finak 等, 2015）是加权整体 DE 工具的有效替代品。MAST 使用障碍模型来处理掉失，同时建模基因表达随条件和技术协变量的变化。在上述研究中，它是表现最佳的单细胞 DE 测试方法（Soneson & Robinson, 2018），并在单个数据集的小规模比较中优于整体和单细胞方法（Vieth 等, 2017）。MAST 的运行时间比加权整体方法快 10 到 100 倍（Van den Berge 等, 2018），而使用 limma–voom（Law 等, 2014）可进一步加速 10 倍。尽管 limma 是一种整体 DE 测试方法，但 limma–voom 在性能上与 MAST 相当。

由于 DE 测试应在未经校正的测量数据上进行，考虑混杂因素对稳健估计差异表达基因至关重要。尽管 DE 测试工具通常允许用户灵活地将混杂因素纳入模型，用户必须谨慎选择添加到模型中的变量。例如，在大多数单细胞实验设置中，样本和条件协变量相互混杂，因为很少可能在多个条件下获得单一样本。如果将样本和条件协变量都纳入模型，则无法明确分配与这些协变量相关的变异。因此，在条件检验时，不能以给定形式将样本协变量纳入模型。当校正多个分类批次协变量时，目视观察协变量的混杂群体变得愈发困难。在此情况下，测试模型设计矩阵是否为满秩很有帮助。即使设计矩阵不是满秩的，DE 测试工具通常会调整矩阵并在不发出警告的情况下运行，这将无法提供预期结果。

在这里描述的情景中，条件协变量由实验设置确定。因此，对该协变量（在同一簇内）的 DE 检验独立于聚类过程。此设置区分了条件间的 DE 测试和簇间的 DE 测试。条件间 DE 测试获得的 P 值表示预期的显著性度量，且需进行多重检验校正。为了减轻多重检验负担，可以从数据集中排除不感兴趣的转录物。尽管伪基因或非编码 RNA 可能具有信息性（An 等, 2017），但通常在分析中被忽略。

陷阱与建议：

• DE 测试不应在校正数据（去噪、批次校正等）上进行，而应在包含技术协变量的测量数据上进行。

• 用户不应依赖 DE 测试工具来校正带有混杂协变量的模型。模型指定应仔细进行，以确保设计矩阵为满秩。

• 我们推荐使用 MAST 或 limma 进行 DE 测试。

基因集合分析

基因层级分析方法经常产生大量难以解释的候选基因列表。例如，数千个基因可能在处理与对照细胞之间存在差异表达。通过根据共享特征将基因分组并测试这些特征是否在候选基因列表中过度表示，可以简化这些结果的解释。

基因集合信息可以在各种应用的整理标签数据库中找到。为解释 DE 结果，通常基于基因在共同生物过程中的参与将其分组。生物过程标签存储在数据库中，例如 MSigDB（Liberzon 等, 2011）、基因本体（Gene Ontology，Ashburner 等, 2000；The Gene Ontology Consortium, 2017）、以及通路数据库 KEGG（Kanehisa 等, 2017）和 Reactome（Fabregat 等, 2018）。可以使用大量工具测试基因列表的注释富集情况，这些工具在 Huang 等（2009）和 Tarca 等（2013）的综述中进行了比较。

单细胞分析领域的一个新发展是使用成对的基因标签进行配体–受体分析。这里，从受体和配体的表达推断出细胞簇间的相互作用。配体–受体对标签可从最近的 CellPhoneDB（Vento-Tormo 等, 2018）获得，并用于使用统计模型解释簇间高度表达的基因（Zepp 等, 2017；Zhou 等, 2017；Cohen 等, 2018；Vento-Tormo 等, 2018）。

基因调控网络

基因并非独立运作。相反，一个基因的表达水平由与其他基因和小分子的复杂调控相互作用决定。揭示这些调控相互作用是基因调控网络（GRN）推断方法的目标。

GRN 推断基于基因共表达的测量，如相关性、互信息，或通过回归模型（Chen & Mar, 2018）。如果两个基因在考虑所有其他基因作为潜在混杂因素的情况下仍表现出共表达信号，则称这些基因具有因果调控关系。推断基因调控关系与检测轨迹相关的调控基因有关。实际上，一些单细胞 GRN 推断方法使用机械微分方程模型来处理轨迹数据（Ocone 等, 2015；Matsumoto 等, 2017）。

尽管一些 GRN 推断方法专为 scRNA-seq 数据设计（SCONE：Matsumoto 等, 2017；PIDC：Chan 等, 2017；SCENIC：Aibar 等, 2017），但最近的比较研究表明，无论是整体还是单细胞方法在这些数据上表现不佳（Chen & Mar, 2018）。尽管如此，GRN 推断方法可能仍提供有价值的见解来识别生物过程的因果调控因子，但我们建议谨慎使用这些方法。

陷阱与建议：

• 用户应注意推断的调控关系中的不确定性。富含调控关系的基因模块比单独的边缘更可靠。

分析平台

单细胞分析工作流是独立开发工具的集合。为了便捷地在这些工具之间传输数据，围绕一致数据格式开发了单细胞平台。这些平台为构建分析管道提供了基础。目前的可用平台包括 R（McCarthy 等, 2017；Butler 等, 2018）和 Python（Wolf 等, 2018）中的命令行工具，以及具有图形用户界面（GUI）的本地应用程序（Patel, 2018；预印本：Scholz 等, 2018）或网络服务器（Gardeux 等, 2017；Zhu 等, 2017）。平台概览可见 Zhu 等（2017）和 Zappia 等（2018）。

在命令行平台中，Scater（McCarthy 等, 2017）和 Seurat（Butler 等, 2018）可轻松与 R Bioconductor 项目提供的多种分析工具接口。Scater 在 QC 和预处理中表现突出，而 Seurat 是当前最流行、最全面的平台，包含大量工具和教程。最近加入的成员是 scanpy（Wolf 等, 2018），一个成长中的 Python 平台，其在处理更大规模细胞数据方面表现更佳，利用了 Python 中日益增多的工具，这尤其适用于机器学习应用。

图形用户界面平台允许非专业用户构建单细胞分析工作流。用户通常通过预设的工作流引导分析，但也限制了用户的灵活性。这些平台特别适合探索性分析。Granatum（Zhu 等, 2017）和 ASAP（Gardeux 等, 2017）在集成工具方面有所不同，Granatum 包含更丰富的分析方法。作为网络服务器，这两个平台可随时使用，但计算基础设施会限制其处理大型数据集的能力。例如，ASAP 仅在 92 个细胞的数据集上进行了测试。基于网络的 GUI 平台的替代方案包括在本地服务器上运行的 FASTGenomics（预印本：Scholz 等, 2018）、iSEE（Rue-Albrecht 等, 2018）、IS-CellR（Patel, 2018）和 Granatum。这些平台和 GUI 封装器可随本地计算能力扩展。未来，人类细胞图谱门户的持续开发将带来更强大的数据可视化工具，能够处理大量细胞数据。

结论与展望

我们回顾了典型 scRNA-seq 分析工作流的步骤，并在案例教程中实现了这些步骤（https://www.github.com/theislab/single-cell-tutorial）。该教程旨在遵循当前的最佳实践，这些最佳实践基于已有的对比研究。虽然整合个别最佳实践工具并不保证得到最优管道，但我们希望该工作流能代表单细胞分析领域的最新现状，为新入门的研究者提供了合适的起点，也为人类细胞图谱在 scRNA-seq 分析中建立最佳实践的努力作出了贡献（预印本：Regev 等, 2018）。需要注意的是，现有的方法对比不可避免地滞后于最新方法的发展。因此，我们尽可能提及了尚未经过独立评估的新进展。随着未来新工具的开发和更多比较研究的进行，这里提出的具体工具推荐将需要更新，但数据处理阶段的总体考虑应保持不变。

其中两个特别值得关注的进展方向是深度学习工作流和单细胞组学整合，因为它们有可能对分析管道产生颠覆性影响。由于深度学习具有扩展到大数据的灵活性，已彻底改变了计算机视觉和自然语言处理等领域，且在基因组学中开始产生显著影响（Webb, 2018）。最早的 scRNA-seq 应用从降维到去噪逐渐出现（例如 scVis：Ding 等, 2018；scGen：预印本：Lotfollahi 等, 2018；DCA：Eraslan 等, 2019）。最近，深度学习被用于构建一个嵌入式工作流，能够拟合数据、去噪并在模型框架内执行下游分析，如聚类和差异表达（scVI：Lopez 等, 2018）。在此设置下，可以将噪声和批次效应估计包含在下游统计检验中，同时保留对数据中变异的准确估计。这类集成建模方法有望取代当前的分析管道，后者通常是单个工具的组合。

随着单细胞组学技术的进步，对集成组学分析管道的需求将不断增长（Tanay & Regev, 2017）。未来的单细胞平台将需能够处理不同的数据源，如 DNA 甲基化（Smallwood 等, 2014）、染色质可及性（Buenrostro 等, 2015）或蛋白质丰度（Stoeckius 等, 2017），并集成这些多种模态的工具。在此设置下，已无法仅使用单一的读数或计数矩阵作为教程的起点。然而，一些平台已开始适应多模态数据结构，以整合从未剪接和剪接读数数据计算的 RNA 速度（La Manno 等, 2018）。单细胞多组学整合可通过共识聚类方法、多组学因子分析（Argelaguet 等, 2018）或多组学基因调控网络推断（Colomé-Tatché & Theis, 2018）进行。具备这些功能的分析工作流将是下一阶段的发展。我们设想这样的多组学分析工作流将在我们为 scRNA-seq 奠定的基础上构建。