单细胞RNA-seq数据分析最佳实践（上）

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

文章信息

Luecken MD, Theis FJ. Current best practices in single-cell RNA-seq analysis: a tutorial. Mol. Syst. Biol. 2019, 15: e8746.

摘要

single cell RNA-seq 提高了基因表达研究的分辨率，这项技术也带来越来越多的单细胞分析方法。这使得研究者难以驾驭这一多工具格局并从中搭建最新的工作流程来分析自己的数据。在这里，我们详细介绍了典型的单细胞 RNA-seq 数据分析步骤，包括预处理（质量控制、标准化、数据校正、特征选择和降维）以及细胞及基因水平的下游分析。我们根据独立比较研究为这些步骤制定了当前（2019年）最佳实践建议。我们已将这些最佳实践建议整合到工作流中，并将其应用于公共数据集，以进一步说明这些步骤在实践中如何工作。我们的案例研究可参见https://www.github.com/theislab/single-cell-tutorial。这篇综述将作为单细胞新手进入该领域的数据分析流程指南，并帮助现有的研究人员更新他们的分析流程。

关键词:分析流程开发；计算生物学；数据分析教程；单细胞 RNA-seq

概述

近年来，单细胞 RNA 测序 (scRNA-seq) 推进了我们对生物系统的认识。我们已经能够研究斑马鱼、青蛙和涡虫的细胞异质性 (Briggs et al,2018；Plass et al,2018；Wagner et al,2018)，并发现之前被掩盖的细胞群 (Montoro et al,2018；Plasschaert et al,2018)。该技术的巨大潜力促使计算生物学家开发一系列分析工具 (Rostom et al,2017)。尽管该领域正在努力确保单个工具的可用性，但单细胞数据分析中，新手的一个进入障碍（ a barrier of entry）是由于该领域相对不成熟而缺乏标准。在本文中，我们简述目前scRNA-seq 分析的最佳做法，为今后的分析标准化奠定基础。

标准化面临的挑战包括分析方法不断增加（截至 2019 年 3 月 7 日已达 385多种工具）和数据集规模爆炸性增长 (Angerer et al,2017；Zappia et al,2018)。我们正在不断寻找新的方法来使用我们所测得的数据。例如，最近的工具可预测分化中的细胞命运 (La Manno et al,2018)。分析工具的不断改进有利于产生新的科学洞察力，但这也使标准化更加复杂。

标准化的第二个挑战在于技术方面。scRNA-seq 数据的分析工具用各种编程语言，最突出的是 R 和 Python (Zappia et al,2018)。尽管跨环境的工具正在增长（预印：Scholz et al,2018），但编程语言的选择通常也是分析工具之间的一种选择。Seurat (Butler et al,2018)、Scater (McCarthy et al,2017) 或 Scanpy (Wolf et al,2018) 等热门平台提供了开发流程的集成环境，且包含大型分析工具。然而，这些平台仅限于使用各自编程语言开发的工具。通过扩展，语言限制也适用于目前可用的 scRNA-seq 分析教程，其中许多教程围绕上述平台（R 和 bioconductor 工具：https://github.com/drisso/bioc2016singlecell和https://hemberg-lab.github.io/scRNA.seq.Lun 等人，2016b；Seurat；scanpy.

考虑到上述挑战，我们并没有标准化分析流程，而是概述了当前的最佳实践和独立于编程语言的通用工具。我们指导读者完成 scRNA-seq 分析流程的各个步骤（图 1），介绍当前的最佳实践，并讨论分析陷阱提出开放性问题。由于工具的新颖性和缺乏比较，事实上无法确定最佳实践，因此我们列出了流行的可用工具。所概述的步骤从reads或计数矩阵开始，得出潜在分析终点，Lun et al (2016b) 涵盖了早期预处理步骤。整合现有最佳实践的详细案例研究可从我们的 github 获得，网址为：https://github.com/theislab/single-cell-tutorial/。在这里，我们在一个实际的示例工作流中应用了当前的最佳实践来分析公共数据集。分析工作流程用rpy2在 Jupyter-Ipython notebook中集成了 R 和 Python 工具。有了可用的文档，它很容易作为工作流模板进行二次修改。

图 1. 典型的单细胞 RNA-seq 分析工作流程示意图。原始测序数据经过处理和比对，得到计数矩阵，代表工作流程的开始。计数矩阵经过预处理和下游分析。使用 Haber et al (2017) 肠上皮细胞数据的最佳实践工作流程生成子图。

框1:实验性scRNA-seq工作流的关键元素

从生物样本到可分析的单细胞数据需要经过多个步骤。典型的工作流程包括：单细胞解离、单细胞分离、文库构建和测序。对这些阶段的简要概述如下：单细胞实验的起始材料通常以生物组织样本的形式获得。

单细胞悬浮液的制备作为第一步，是在一个被称为单细胞解离的过程中产生的，其中组织被消化。为分析每个细胞中的 mRNA，必须分离细胞。单细胞分离根据实验方案的不同而不同。虽然基于平板的技术将细胞隔离到平板上的孔中，但基于液滴的方法依赖于在自己的微流体液滴中捕获每个细胞。在这两种情况下，都可能发生错误，导致多个细胞被捕获在一起（doublets or multiplets）、非活细胞被捕获或完全没有细胞被捕获（空液滴/孔）形成空滴的情况尤其常见，因为基于液滴的方法依靠低浓度的输入细胞流动来控制双联体率。每孔或液滴中都含有分解细胞膜和进行文库构建所必需的化学物质。胞内 mRNA 被捕获、反转录为 cDNA 分子并扩增的过程称为文库构建。当细胞隔离进行这一过程时，每个细胞的 mRNA 可以被一个孔或滴特定的细胞条形码标记。此外，许多实验方案也用唯一分子标识符 (UMI) 标记捕获的分子。测序前扩增细胞 cDNA，以增加其被测量的概率。UMIs 允许我们区分相同 mRNA 分子的扩增拷贝和从相同基因转录的不同 mRNA 分子的reads。

构建好文库后，使用细胞条形码进行标记，并根据协议进行UMIs标记。这些库汇集在一起(multiplexed)用于测序。序列产生reads数据，这些数据经过质量控制，再准备阶段根据指定的条形码(demultiplexing)和reads比对区分细胞。对于基于umi的协议，reads数据可以被进一步解复用以产生捕获的mRNA分子计数(count data)。

Pre-processing and visualization

对测序仪生成的原始数据进行处理，以获得分子计数（count 矩阵）或读数（reads矩阵）的矩阵，这取决于是否在单细胞文库构建方案中纳入了独特的分子标识符 ( unique molecular identifiers ，UMI)（有关分析前的实验步骤概述，请参见框 1）。Cell Ranger (Zheng et al,2017)、indrops (Klein et al,2015)、SEQC (Azizi et al,2018) 或 zUMIs (Parekh et al,2018) 等原始数据处理流程负责reads质量控制 (QC)，为其细胞barcode和 mRNA 来源分子（也称为解复用,demultiplexing）分配reads、基因组比对和定量。得到的reads或计数矩阵包含barcode x 转录本数量的高纬数据。此处使用术语barcode代替细胞，因为所有reads均为分配给相同的barcode可能与来自同一细胞的reads不一致。一个barcode可能错误地标记多个细胞（双联体）或可能不标记任何细胞（空滴/孔）。虽然reads和计数数据的测量噪声水平不同，但典型分析流程中的处理步骤相同。为了简单起见，我们将在本教程中将这些数据称为count矩阵。如果reads和count矩阵的结果不同，则专门指出reads矩阵。

Quality control

在分析单细胞基因表达数据之前，我们必须确保所有的细胞barcode数据都对应于活细胞。细胞 QC 通常基于三个 QC 变量进行：

• 每个barcode的计数数量（count depth ）
• 每个barcode的基因数量
• 每个barcode的线粒体基因计数分数 (Ilicic et al,2016；Griffiths et al,2018)

检查这些 QC 变量的分布，以确定是否存在通过阈值处理过滤掉的离群峰（图 2）。这些异常barcode可能对应于死细胞、膜破损的细胞或双联体。例如，低计数深度的barcode、很少检测到的基因以及线粒体计数的高分数都表明细胞的细胞质 mRNA 已经通过破损的膜漏出，只有位于线粒体中的 mRNA 仍然是保守的（图 2）。与之相反，非预期高计数和大量检测基因的细胞可能代表双联体。因此，高计数深度阈值常用于过滤掉潜在的双峰。最近的三种双联检测工具提供了更优雅和可能更好的解决方案 (DoubletDecon:preprint:DePasquale et al,2018；Scrublet:Wolock et al,2019；doublet Finder:McGinnis et al,2018)。

图 2. Haber et al (2017) 的小鼠肠上皮数据集过滤决策的质量控制指标图。(A) 每个cell的计数深度直方图。较小的直方图在计数深度低于 4,000 时放大。根据在约 1,200 个计数处检测到的峰值，此处应用的阈值为 1,500。(B) 每个细胞检测到的基因数的直方图。在大约 400 个基因处可见一个小的噪声峰。这些细胞使用描述的阈值（红线）700 个基因过滤掉。计数深度分布从高到低计数深度。该可视化与 Cell Ranger 输出中显示的 logClog 图相关，该输出用于过滤空液滴。它显示了一个肘部的计数深度开始迅速减少约 1500 计数。(D) 通过线粒体读数部分染色的基因数量与计数深度的关系。线粒体读取片段仅在检测基因很少的特别低计数细胞中高。这些细胞被我们的计数和基因数阈值过滤掉。联合可视化计数和基因阈值显示联合过滤效果，表明较低的基因阈值可能已经足够

单独考虑这三个细胞 QC 变量中的任何一个都可能导致对细胞信号的误解。例如，具有较高线粒体计数的细胞可能参与呼吸过程。同样，其他 QC 变量也有生物学解释。低count和（或）基因的细胞可对应静止细胞群，高count的细胞体积可能更大。事实上，细胞之间的分子计数可能存在强烈差异（参见项目 github 的案例研究）。因此，当单变量阈值决策时，应联合考虑细胞 QC 变量（图 2D），这些阈值应尽可能设置为允许的，以避免无意中过滤掉活细胞群。考虑到多变量细胞 QC 的依赖性，筛选模型可能提供更敏感的 QC 选项。

含有异质混合细胞类型的数据集可能显示多个 细胞QC 变量峰值。例如，图 2D 显示了具有不同 QC 分布的两个细胞群。如果之前没有进行过滤步骤（注意 Cell Ranger 也进行细胞 QC），那么只有每个barcode峰的最低计数深度和基因应该被认为是非活细胞。进一步的阈值指导原则是使用所选阈值过滤掉的细胞比例。对于高计数过滤，该比例不应超过预期的双联率。

除了检查细胞的完整性，细胞 QC 步骤也必须在转录本水平上进行。原始计数基质通常超过 20,000 个基因。通过过滤掉在少数细胞中不表达的基因，可以大幅减少这一数量。设置此阈值的一个准则是使用最小cell群，并留下一些dropout 效应（dropout effects. ）的余地。例如，筛选出少于 20 个细胞中表达的基因可能会使检测少于 20 个细胞的细胞团变得困难。对于高脱落（dropout ）率的数据集，这个阈值也可能使较大簇的检测复杂化。阈值的选择应根据数据集中的细胞数量和预期的下游分析进行调整。

可直接对计数数据进行进一步 QC。Ambient gene expression（环境基因表达）指不是来自barcode细胞，而是来自其他溶解细胞的count，这些细胞的 mRNA 在文库构建之前污染了细胞悬液。这些增加的环境计数会扭曲下游分析，如标记基因鉴定或其他差异表达检测，尤其是当样本之间的水平变化时。在基于液滴的 scRNA-seq 数据集中校正这些影响是可能的，由于大量的空液滴，可用于模拟环境RNA表达谱。最近开发的SoupX(预印本:Young & 使用这种方法直接纠正计数数据。在下游分析中忽视强环境基因的实用方法也被用来解决这个问题(Ange- lidis et al, 2019)。

进行质量控制以确保数据质量足以用于下游分析。由于无法先验确定足够的数据质量，因此根据下游分析性能（例如，聚类注释）进行判断。在分析数据时，可能需要多次重新审查质量控制参数。通常，从允许的质控阈值开始，在执行更严格的质控之前研究这些阈值的影响是有益的。这种方法对于包含异质性细胞群的数据集特别重要，其中细胞类型或状态可能被错误解释为低质量离群细胞。在低质量数据集中，严格的 QC 阈值可能是必要的。可通过试验 QC 指标确定数据集的质量（见附录补充文本 S2，卑微小王手头并没有补充文档，从略）。在这种迭代 QC 优化中，应该注意数据窥视（data peeking.）。不应调整 QC 阈值以改善统计检验的结果。相反，可根据数据集可视化和聚类中的 QC 变量分布来评价 QC效用。

问题和建议: •通过基因数量、计数深度和线粒体reads分数的异常峰来执行细胞QC。考虑这些共同的影响而不是单独的考虑它们。 •尽可能地容忍QC阈值化，如果下游聚类无法解释，则重新QC。 •如果QC变量在样品之间的分布不同，则应针对每个样品分别QC，以解释样品质量差异，如Plasschaert等(2018)。