文献名称：教程：单细胞 RNA 测序数据的计算分析指南
文献期刊：nature protocols
发表时间：2021-1-9

摘要

单细胞 RNA 测序技术 (scRNA-seq) 是一种流行且功能强大的技术，可以分析大量单个细胞的整个转录表达谱。然而，这些针对这些单细胞实验产生的大量数据的分析，需要专门的统计和计算方法。本篇文章，我们概述涉及处理 scRNA-seq 数据的计算工作流程。我们讨论了一些最常见的分析任务和可用于解决核心生物学问题的工具。在本文和我们的配套网站 (https://scrnaseq course.cog.sanger.ac.uk/website/index.html) 中，我们提供了执行单细胞相关计算分析的最佳实践的指南。本教程针对对他们自身数据感兴的实验科研工作者以及寻求开发新计算方法的生物信息学家。

简介

scRNA-seq 已成为一种革命性的技术，用于表征复杂组织并回答bulk RNA-seq无法解决的问题。自 2009年发布第一个 scRNA-seq实验以来，已经发布了许多单细胞实验技术和商业平台。当前两大主要的 scRNA-seq技术平台是：
1.最常见的方法是使用微滴（microscopic droplets）或孔板（wells）来分离大量细胞，然后对文库进行测序深度相对较浅地测序。为了识别给定转录本来自哪个细胞，这些方法使用细胞条形码（cellular barcodes，连接到每个read的短核苷酸标签，对于液滴或孔板是唯一的）。当前主流的 10× Chromium 实验平台是典型的高通量、低测序深度的范例。该技术的一个重要优势是它支持唯一分子标识符 (UMI)。UMI 是在PCR扩增前附加到转录本上的短条形码，从而可以去除PCR聚合酶链反应重复拷贝，消除PCR扩增偏好性，获得更准确的基因表达水平估计。 该技术一个主要缺点是该平台仅允许对每个信使 RNA (mRNA) 的 5' 或 3' 端进行测序。
2.但是，许多研究采取了相反的策略，即分离相对较少的细胞，但对它们进行更深的测序。这些低通量、高测序深度的实验通常将细胞分离到单个孔中并应用 Smart-seq2技术。除了最近引入的 Smart-seq3 实验技术外，这些方法不支持 UMI，但它们通常表现出比基于液滴的技术更高的灵敏度，并且它们还允许对整个转录本进行分析。有关不同平台的详细概述，请参阅最新的综述和基准测试。

除了优化实验工作流程之外，最近的创新还大大减少了scRNA-seq中每个细胞的成本。因此，用于分析的细胞数量呈指数增长。鉴于生成的大量数据，单细胞数据分析需要有效的计算和统计方法。随着实验技术的迅速改进，用于处理数据的计算工作流程也得到了改进。本教程的目的是概述最常见的 scRNA-seq 数据分析类型。本文旨在作为我们为讲授 scRNA-seq 数据的计算分析而开发的课程材料 (https://scrnaseq-course.cog.sanger.ac.uk/website/index.html) 的配套材料。该网站于 2016 年首次推出，并不断更新以包含新方法并提供有关最佳实践的最新建议。

scRNA-seq分析的一个核心组成部分是基因表达矩阵，它表示每个基因种每个细胞观察到的转录本数量。工作流程可分为两个主要部分：1）生成基因表达矩阵；2）基因表达矩阵的分析（图1和表1）。虽然我们的在线教程包含这两个方面，但在这里，我们将重点介绍在获得基因表达矩阵后执行的分析类型。大多数基因仅在一部分细胞类型表达，但是，由于起始材料的含量较低，以及scRNA-seq实验中常用的测序深度较低，一些基因即使表达也无法检测到。导致基因表达矩阵中存在大量零值，这是有问题的，因为一些零值可能代表细胞中的实际低表达或零表达以及测量过程的变化。区分和适当建模这些观察到的零值来源的困难是计算分析的主要挑战之一。即使深度测序的数据集也可能有约50%的零值，而低测序深度的数据集可能有99%的零值。相比之下，在bulk RNA-seq数据集中，小于20%的数据条目为零值。

分析流程概述

数据质控

分析 scRNA-seq 的第一步是剔除不太可能代表完整单个细胞的细胞条形码。对于高通量方法，关键步骤是过滤掉不代表细胞的细胞条形码。最直接的方法是计算特定于数据集的细胞所需的最小 UMI 数量（阈值），高于阈值的将此类条形码视为细胞。最近开发的一些工具，例如 EmptyDrops，首先估计存在于空孔或液滴中RNA的背景水平，然后识别与背景显著偏离的细胞条形码，将此视为细胞。这种策略的优势在于，它能够检测相对于样本中其他细胞而言，RNA含量较低的细胞类型。

不幸的是，这些方法都不能将完整的活细胞与受损或垂死的细胞区分开来。必须进行第二轮数据质控，考虑检测到的基因数量、来自线粒体基因组的RNA比例以及每个细胞比对不上或多重比对read的比例。线粒体来源的基因比例高、检测到的基因少或未比对或多重比对read比例高的细胞通常会受损或死亡。具体阈值通常通过手动检查数据质控指标图来确定，因为最佳阈值取决于组织、细胞解离方案和其他技术因素。为关键指标定义异常值细胞（根据中值绝对偏差）能够直接构建数据集特定的阈值，但应谨慎应用，特别是对于包含高异质细胞类型的样本 。

除了一些表示背景噪声的细胞条形码外，细胞条形码也可能对应于多个细胞。通常，约5%的细胞条形码标记多个细胞，称为doublets。此外，最近的研究结果表明，多达20%的情况下，多个细胞条形码可能标记同一个单细胞，称为多重条形码（barcode multiplets）。scrublet和DoubletFinder 等工具从数据集本身模拟可能的双细胞，然后计算真实液滴条形码与模拟双细胞的相似性，并定义阈值以区分推断的双细胞和假设的单细胞。其他方法也可以成功应用（例如scds），但双细胞检测是一个复杂的问题，并且不能期望计算双细胞检测方法能够在所有实验设计中完美执行。

问题：多个cell barcode对应单个细胞；这种的处理方法是什么？

查看《Inference and effects of barcode multiplets in droplet-based single-cell assays》文献：
基于液滴的分流系统已成为单细胞基因组学研究必不可少的工具。与基于板的单细胞测定相比，基于液滴的方法，包括scRNA-seq和 scATAC-seq技术，可以在单个实验分析数千个细胞。检测细胞通量的显着增加是通过使用bead对细胞cDNA进行平行条形码来实现高度多样性的DNA条形码。至关重要的是，下游的计算分析假设一个条形码序列等同于一个细胞。
在这项工作中，我们提供了多条证据表明细胞通常通过（i）多个beads出现在同一个液滴内；(ii)单个bead内存在异质性的寡核苷酸序列（图1a）。这里，我们指的是在同一液滴中出现多个DNA条形码的情况，即“barcode multiplets”。我们发现多重条形码会显著影响单细胞分析和证明罕见的细胞事件（例如，细胞克隆型分析）。此外，我们提供了针对现有的单细胞数据集，特别是scATAC-seq平台，识别这些多重条形码的计算解决方案。最后，我们提供了缓解现有分析中这类多重barcode偏差的建议。

归一化

从测序实验中获得的有用read数将因细胞而异，必须纠正这种差异。对于scRNA-seq数据，这种影响是显著的，因为每个细胞的RNA的数量可能会因细胞周期阶段和其他生物因素而显著变化，即使在相同的细胞类型内也是如此。技术因素（例如，不同的液滴大小）可能会进一步增加测序深度的可变性。由测序深度不均一而产生的差异可以通过归一化来改善。

针对bulk RNA-seq数据，归一化相当于计算与样本测序深度相关的数量，通常称为“size factor”，并将所有基因的表达计数除以该值。原则上，类似的方法可以用于scRNA-seq，但大量的零意味着需要修改策略。scran包通过使用细胞池来估计size factor实现稳健的结果。或者，可以使用外部来源的RNA对照组或看家基因的RNA表达计数来估计size factor。

由于大量的零值，低表达基因的转录行为可能与高表达基因不同，以响应不同的测序深度。为了补偿这种行为，可以使用特定于每个基因表达水平的归一化策略。例如，SCnorm可用于低通量、高测序深度数据，sctransform可用于高通量、低测序深度的数据。2019年，开发了一种新的用于scRNA-seq基因表达计数缩放和推断的贝叶斯方法，称为bayNorm，其目的是在考虑mRNA捕获的影响后估计潜在的基因表达矩阵。

批次校正

与测序深度的差异类似，批次效应是必须考虑的技术混杂因素，这样才能出现真正的生物信号。批次效应是生物学中的一个常见问题，它们源于非生物学因素的差异，例如实验时间、进行实验的人或试剂的差异。如果没有适当考虑，批次效应可能会被误认为是真正的生物信号，但是，通过仔细的实验设计，它们可以完全避免。要将批次效应校正应用于数据集，不能混淆实验（即每个批次必须包含至少两个生物条件）。当在所有批次中处理所有生物条件时，批次效应校正最有效，称为“平衡设计”。不幸的是，当样本不能同时处理时（例如，如果细胞在收集后需要立即处理），通常不可能实现平衡设计。

传统的校正批次的方法，如ComBat，假设每个细胞的生物学状况是先验的，并利用此信息使用线性模型将生物效应与批量效应分离。然而，这种假设通常不适用于scRNA序列数据，因为单个细胞的细胞类型身份可能未知。为了应对这一挑战，mnnCorrect使用不同批次中细胞之间的相互最近邻来识别批次后的常见生物状况。这种相互最近邻方法也适用于为Seurat的典型相关分析（CCA）方法找到“锚”。这两种工具之间的主要区别是，mnnCorrect使用PCA从基因表达矩阵中删除批量效应，而CCA将细胞投影到一个共同的基因相关空间，并在该空间上执行校正。然而，即使是这些单细胞特定工具也假设跨批次共享生物条件，如果应用于一个有争议的实验，也会错误地去除真实的生物信号。

缺值填补和平滑

许多归一化策略不会更改零值，因此很容易假设它们代表缺失值，并从检测到的转录本通过数学推导获得的估计值进行填充。原则上，删除零值可以减少噪音，并使其更容易识别数据的潜在结构（例如，基因-基因相关性、细胞cluster、标记基因或发育轨迹）。

已经开发了几种工具来“填补”在scRNA-seq数据中发现的零值，包括scImpute, DrImpute和SAVER。DrImpute和scImpute的性能类似，而SAVER对数据的影响较小，产生的虚假信号更少。这些工具都依赖于在数据中找到可用于预测缺失值表达水平的数据结构。然而，这些方法假设数据集中的所有基因都由已识别的结构（已有的基因表达矩阵）确定，经常导致引入大量假阳性信号。

其他工具，如使用扩散模型的MAGIC和使用自动编码器的scVI，应用平滑算法来减少噪声。因此，这些方法采用数据驱动的方法，假设缺失值可以从具有类似基因表达谱的其他细胞中推断出来。与基于模型的零值填补方法类似，它们可以更容易地检测下游分析中的结构。这种算法的一个缺点是，基础模型可能会扭曲真实结构（例如，通过放大随机噪声），这可能会被误认为是生物模式。随着公开可用的单细胞图谱数量的增加，使用外部引用来插补缺失值变得可行。这种方法的示例是SAVER-X和netNMF-sc，它们可能没有相同的缺点，因为它们能够合并来自其他来源的相关信息。零值填补有助于改善scRNA-seq数据的可视化，但填补的数据确定的任何结构或模式（例如差异表达基因或轨迹）必须通过对预进行填补的数据应用适当的统计检测进行验证。

细胞周期分配

如果样本中含有活性细胞周期的细胞，这可能会导致生物学混杂因素，可能需要在下游分析中去除该类细胞。或者，细胞周期的阶段可能与正在研究的生物学问题有关。在任何一种情况下，都有必要将细胞分配到其适当的细胞周期阶段。有两种广泛使用的工具可用于识别细胞周期阶段：cyclone和Seurat。Cyclone分析在不同相对水平表达值的成对基因，以将细胞分配到G1、S或G2/M。无论是否归一化，cyclone都很精确，但它很难区分非细胞周期细胞。Seurat采用Tirosh等人提出的方法，根据已知标记基因G1/S和G2/M的平均归一化表达值对细胞进行评分。

一旦细胞被分配了一种细胞周期阶段，两种工具都使用一般线性模型来回归差异。此外，Seurat提供了一个选项，仅回归掉G1/S和G2/M细胞之间的差异，同时保留细胞周期的细胞和非细胞周期的细胞之间的差异。如果有人对细胞周期和非细胞周期的细胞亚群之间的差异感兴趣，后一种情况很重要。

高可变基因选择

在 scRNA-seq 实验中，每个基因代表一个维度，因此，对于小鼠或人类数据集，将有大约 20,000 个维度。然而，许多基因不会在给定的细胞或细胞类型中表达，并且在实验中检测到的数量取决于不同的实验技术。高通量、基于液滴的方法可以识别细胞中多达 5,000个基因，而更灵敏的方法可以检测较前面两倍的基因。然而，许多研究依赖于低深度测序，因此检测到的基因较少，有时每个细胞少于 1,000 个基因。基因数目过多会使分析变得困难，因为高维度的距离估计是不可靠的，即使在噪声水平较低的情况下也是如此。

高可变基因识别是相对于技术噪声具有很强的生物学信号。因此，在下游分析只限制在信息量最大的基因集上，可以减少维度的影响，减少噪音并简化分析。一些工具可能会识别固定数目的高可变基因，或尝试确定哪些特征基因包含大量的生物学信息。scRNA-seq 数据中的特征基因选择有两个复杂的因素：(i) 影响每个基因的技术噪声取决于该基因的平均表达；(ii) 小样本量难以估计方差。最广泛使用的特征基因选择策略是考虑高度可变的基因（即具有高于预期方差的基因）。对于使用 UMI 量化的数千个细胞的数据集，已表明噪声遵循负二项分布，可用于识别重要特征基因。 Seurat等工具使用非参数方法通过经验拟合方差和基因平均表达之间的关系来识别高可变的基因。相反，另一种特征选择策略是考虑具有高于预期数量的观察到的零值的基因。

许多特征选择方法的一个局限性是，它们考虑了整个数据集的整体可变性。因此，在罕见细胞类型中差异表达的基因可能无法检测到，因为这些细胞对总变异性的贡献很小。在这种情况下，其他指标，如量化转录本不均匀分布的基尼指数，可能更合适，如基尼聚类方法所示，该方法旨在识别小聚类。

降维和可视化

另一种降低基因表达矩阵高维负面效应的策略是对缩减的特征空间进行降维。有许多可用的方法 ，但最常用的策略涉及主成分分析 (PCA)，这是一种线性变换，可保留完整 PCA 空间中细胞之间的欧几里德距离，即使对于非常大的数据集也可以有效计算。为下游分析保留的组分数量将取决于数据集的复杂性，并且存在各种算法来识别适当的PC个数。由于这些通常在计算上运行起来很耗时，因此最常见的方法是绘制每个分量可解释的方差分数，然后直观地识别曲线急剧弯曲的点，通常称为“拐点”，并且只保留拐点以上的那些PC组分。

大多数 scRNA-seq 数据集都很复杂，它们的结构不能被两个或三个主成分捕获。因此，可视化算法用于创建一个二维图，汇总scRNA-seq数据集的大量的重要组成部分。当前的最佳实践方法是统一流形逼近与投影UMAP的降维方式。该算法使用细胞-细胞最近邻网络近似数据的拓扑结构，然后估计数据的低维嵌入，以最好地保留结构。UMAP 在很大程度上取代了 t 分布随机指标邻居嵌入降维方法（t-SNE）， 因为它能够更好地保存全局结构。最近的一项研究表明，这可以归因于UMAP实现中默认使用的初始化策略。然而，UMAP倾向于支持数据的完全性连接表示，而不是t-SNE支持的离散聚类。t-SNE和UMAP的一个缺点是它们都需要用户定义的超参数，结果可能对选择的参数数值敏感。此外，这些方法是随机的，提供良好的初始化可以显著改善这两种算法的结果。需要注意的是，无论是可视化算法保留了细胞之间的距离，因此，投影坐标的信息不能直接用于聚类或拟时序分析等下游分析。

无监督聚类

scRNA-seq 数据的无监督聚类是大多数分析的核心，因为它可以识别具有相似表达谱的细胞组。其中一些组可以代表不同的细胞类型，而其他组可以被认为是中间细胞状态（例如，细胞周期阶段），这取决于样本的生物学系统。早在 scRNA-seq 出现之前就已经开发了各种聚类方法，现有的工具都是经典方法的应用。一个例子是广泛使用的 k-means 算法 ，它构成了单细胞共识聚类 (SC3) 算法的基础。除了基本的 k-means 算法之外，SC3 还使用一种共识方法来平均多个聚类结果。另一个例子是用于网络聚类的 Louvain 算法 rithm，它成功地适应了Phenograph中的单细胞数据集，随后被 Seurat和 scanpy采用。基于Louvain算法的方法为细胞构建最近邻网络，然后识别网络中的不同社区。这些方法的优势在于它们的速度，即使对于非常大的数据集也是如此。独立的基准比较表明 SC3 和 Seurat 总体上表现相似，尽管对于单个数据集，其中一个或另一个可能表现更好，并且优于所有其他当前可用的方法。这些基准基于数据集，其中细胞类型身份可以通过转录组分析以外的其他方式（例如，已知表面标记的荧光激活细胞分选分析）建立。

每个聚类算法都有自己的一组参数，这些参数可以显著影响结果和生物学解释。例如，Louvain 算法有一个resolution参数会影响cluster的大小——较小的分辨率会生成较少数目的cluster。同样，对于k-means的方法，k 的值直接决定了聚类的数量。不幸的是，没有确定最佳参数的固定规则，用户通常必须根据手头的数据集做出明智的决定。重要的是要同时考虑数学和生物学方面，因为仅依靠一个标准的分析流程可能会导致不能提供最适合该数据集的结果，或者cluster的分群不合理。例如，可以计算cluster对输入参数的稳健性，并检查生成的cluster中已知存在于特定组织中的细胞类型。

拟时序分析

聚类分析将每个细胞分配给一个组（cluster），这在某些情况下是不合适的。例如，如果数据集代表一个发育过程或源自一个时间过程实验，那么将细胞视为从一个发育/时序连续体中绘制出来的更合适。这种可以表示空间位置、化学浓度或时间进程的连续轨迹通常被称为“伪时间”，每个细胞都可以被分配一个特定的位置。大多数工具无法确定细胞沿轨迹移动的方向或速度。相反，必须使用外部信息，例如时间过程实验的采样时间或发育轨迹的标记基因来推断这些数量。大量的伪时间推理方法已经发表，最近也有研究已经进行了基准测试。作者强调这些方法是互补的，它们为不同类型的数据选择哪种方法提供了指导。

大多数工具采用两种方法中的一种。第一种方法是基于流形方法的拟时序方法。使用降维技术来识别细胞所在的低维“流形”，并使用细胞-细胞图来描述流形的拓扑。使用这种策略的主流方法包括 Monocle和 DPT。第二种方法是基于cluster的聚类方法。在连接cluster并将单个细胞投影到cluster分支上之前，使用无监督聚类对细胞进行分组。此类方法的示例包括 TSCAN 和 Mpath。当进行轨迹分析的细胞密度不相等时，基于cluster的拟时序方法往往更准确——例如，来自一种状态的细胞可能比来自其他状态和大规模发育层次结构的细胞更频繁或更可靠地捕获。另一方面，当在过渡过程中对细胞进行均匀采样以及检查奇异过渡的细节时，流形方法表现最佳。

表示剪接和未剪接转录本的外显子和内含子read的相对丰度可用于推断 scRNA-seq 实验中的时间动态。 RNAvelocity 和 scVelo6等工具可以推断在对细胞进行采样时每个基因的表达是增加还是减少。尽管 RNAvelocity 使用简单的动态模型，但 scVelo 采用概率方法来解释单细胞数据中的不确定性。尽管这种方法受到测序深度和比对到内含子的读取数量的限制，但它可以推断每个细胞在表达空间中移动的方向以及对变化率的估计。结果可以在低维投影中可视化为指示每个细胞如何移动的箭头，类似于相平面。

差异基因表达

差异基因表达（DE）已成为bulk RNA-seq中最重要的应用之一，因为它提供了一系列在两个或多个生物条件之间受干扰的基因。scRNA-seq的DE分析更具挑战性，因为我们不仅仅是比较每个基因的单个值，而是需要比较基因表达水平的分布。单细胞数据特有的另一个挑战是，我们想要比较的细胞组不是先验定义的。相反，这些组通常是基于我们想要比较的基因表达水平来定义的，这违反了统计标准假设检验中的中心假设。事实上，已经证明，无监督聚类和差异表达分析可以导致人为的低P值。因为在定义细胞组别时使用了基因表达值，这可能会引入偏差，因为聚类和DE分析不再独立。

最近的一项比较得出结论，与目的构建方法相比，非参数Wilcoxon检验表现出色。作者还得出结论，为bulk RNA-seq开发的方法表现良好，尤其是当与基因表达矩阵的每个元素分配权重的策略相结合时。另一项基准研究得出了类似的结论，补充说基因表达量的标准化可以对结果产生重要影响。在专门为scRNA-seq量身定制的方法中，MAST，它使用高斯障碍模型将检测率差异和平均表达差异结合到一个测试中，已被报告具有最佳性能。

当单细胞实验包含多个生物学重复时，就会出现一个有趣的情况（例如，比较3个健康个体的细胞与3个患有糖尿病的个体的细胞）。当前的单细胞差异表达检验将每个单独的细胞视为生物学重复，不能解释共享的遗传背景或疾病状态。对于此类比较，当前的可能的操作是： (i) 计算每个细胞类型的每个个体的细胞间基因平均表达量，并将结果视为bulk RNA-seq样本；或 (ii) 执行所有个体 × 个体双重比较和过滤出单个个体所独有的差异基因结果。前一种方法类似于最近提出的 MetaCell 想法。 MetaCell 背后的想法是使用引导方法来识别原始数据集最稳定和可重现的特征。然而，随着 scRNA-seq 应用于更大的队列和比较研究，我们预计进一步的发展将导致更准确的统计模型用于更复杂的实验设计。

比较与整合数据集

随着 scRNA-seq 数据量的不断增长，重要的挑战是如何最好地整合单细胞数据集。批次效应是整合数据的主要挑战是数据来自不同实验室的实验，即使这些问题可以克服，整合这些数据集的重分析可能需要耗费大量的时间、精力和数据储存。相反，整合数据的一种替代策略是比较它们。当其中一个数据集非常大（例如细胞图谱）时，该策略特别有用。当给定一个或多个具有已知细胞类型的数据集时，scmap构建一个小索引。当给定一个新的查询数据集时，scmap根据转录谱快速确定识别心数据集的每个细胞与参考数据集的哪个细胞类型最接近。此外，scmap可以预测参考数据集中最紧邻细胞，这意味着当为细胞分配伪时间数值而不是离散cluster标签时，可以使用它。也就是说scmap可以做拟时序分析。最近有文献研究，对将细胞映射到参考数据集的不同方法进行了基准测试。另一种方法MetaNeighbor，旨在测试细胞类型在多个scRNA-seq数据集中是否一致。它通过计算跨数据集的细胞-细胞Spearman相关性来实现，允许MetaNeighbor验证细胞标签在多个实验数据中的重复性。

总结

scRNA-seq的计算分析是一个快速发展的领域。在很大程度上，这是由新平台和实验技术的开发推动的。然而，研究人员也提出了从数据中提取信息的新方法。未来几年可能会出现新的分析工具，进一步扩大scRNA-seq的使用。此外，我们还希望整合分析工作流的软件工具（如Seurat、scanpy和Bioconductor）会有所改进，使生物信息知识有限的用户更容易访问分析。

新型单细胞技术的快速发展，尤其是可以分析细胞的多组学方法和提供空间信息的方法，将需要新的计算方法来充分利用数据。此外，各种图谱项目产生的数据量不断增加，这将需要更适合于分析大细胞量的方法。