单细胞转录组 —— 原始数据处理

单细胞转录组 —— 原始数据处理

前言

文章参考 Single-cell best practices 中原始数据处理部分。

我们这里所说的 sc/snRNA-seq 原始数据处理涉及多个步骤，从测序仪下机产生的 BCL 文件到原始的 FASTQ 文件，再到序列比对，最后以 reads 计数矩阵结束。

计数矩阵代表了每个细胞中每个基因所产生的不同分子数量的估计值，并可以根据每个分子推测的剪接状态对数据进行分层。

整体流程如下图

这个计数矩阵是所有下游分析的起点，包括数据归一化、整合和过滤方法，以及推断细胞类型、发育轨迹和表达动态的方法等。

因此对该计数矩阵进行稳健而准确的评估，是实现准确可靠的后续分析的基础和关键步骤。

在原始数据处理中的错误会导致后续分析中出现遗漏、削弱或扭曲原本存在的生物信息。

我们主要介绍原始数据处理的基本步骤，包括序列比对、细胞条形码（CB）识别和校正，以及通过解析唯一分子标识符（UMI）估算分子数。

注意

我们所说的原始数据一般指的是 FASTQ 文件，但是在 FASTQ 生成的上游也可能出现问题，如索引跳转，不过这些问题可以通过 in silico 方法和增强协议（如双索引）在很大程度上得到缓解。在本节中，我们将不考虑这些上游影响。

原始数据的质量控制

获得原始 FASTQ 文件后，可通过运行质量控制（QC）工具（如 FastQC）来评估读数质量、序列内容等，从而快速诊断 reads 的质量。

如果运行成功，会为每个输入文件生成一份 QC 报告。该报告总结了质量得分、碱基含量和某些相关统计数据，以发现源于文库制备或测序的潜在问题。

虽然现在的单细胞原始数据处理工具都会在内部评估一些对单细胞数据很重要的质量检查，如序列中 N 的含量和比对上序列的百分比等，但在处理数据前运行快速质量检查仍然是一个好习惯，因为它可以评估其他有用的质量控制指标。

序列比对

目前大部分序列比对算法都是基于 bulk 测序数据，但是单细胞测序的原始序列文件包含细胞的信息（如 CB、UMI）。

使用这些比对工具之前，需要对这些数据进行拆分提取，使得用于比对序列只包含原始的序列片段，同时要保证比对后还包含这些细胞和分子的信息。

这就需要使用单独的工具来执行比对前的处理，这种额外的开销可能会造成计算上的负担。此外，存储中间文件通常需要大量磁盘空间，而处理这些中间文件所需的额外输入和输出操作又进一步增加了运行时间开销。

为此，人们开发了几种专门用于比对 scRNA-seq 数据的工具，它们能自动和/或在内部处理这些额外的处理要求。

包括 Cell Ranger、zUMIs、alevin、RainDrop、kallisto|bustools、STARsolo 和 alevin-fry 等

提供了专门的处理方式用于比对 scRNA-seq 数据，可以解析 reads 中的 CBs 和 UMIs，同时也可以对 UMI 进行解析和去重。

这些工具在内部使用了各种不同的方法，有些工具生成传统的中间文件（如 BAM 文件）并随后对其进行处理，有些工具则完全在内存中工作，或使用简化的中间表示法以减轻输入/输出负担。

虽然不同比对方法的具体算法、数据结构以及时间和空间权衡可能会有很大不同，但我们可以对现有方法进行粗略分类：

• 比对的类型
• 参考序列的类型

不同比对类型

常用的比对算法主要分为三类：剪切比对、连续比对和轻量级映射。

比对方法使用各种启发式方法确定 reads 可能来源的潜在位置，然后在每个假定位置计算 reads 与参考序列之间的最佳核苷酸级别比对，通常使用动态规划算法。

动态规划算法有多种类型，适用于不同类型的比对，例如全局比对、局部比对和半全局比对。

全局比对适用于整个查询序列与参考序列对齐的情况，而局部比对则适用于查询序列的一个子序列与参考序列的一个子序列对齐的情况。半全局比对则是介于两者之间的情况，通常用于短读长比对。

此外，剪切比对（softclip）允许忽略或减轻在 reads 序列开头或结尾处出现的不匹配、插入或删除的惩罚。

为了提高比对过程的效率和实用性，设计了多种复杂的修改和启发式方法。例如，Banded alignment 通过限制动态规划表格中的计算区域，将计算范围缩小到对角线附近的一个带状区域，从而降低了计算量和内存使用。

其他如 X-drop 和 Z-drop 启发式方法则在早期将不太可能成功的比对进行剪枝。

wavefront 比对等最新算法能够在较短时间和较小空间内确定最佳比对。

比对得分和实际比对的回溯通常以 CIGAR 字符串形式编码在 SAM 或 BAM 文件中。

CIGAR 字符串是一种压缩的字母数字表示方式，描述了比对的具体情况，例如 3M2D4M 表示三个匹配或不匹配，随后有两个碱基中的删失，再就是四个匹配或不匹配。

通过计算这些比对得分，能够过滤出与参考序列对齐良好的读取，避免低复杂度或假匹配。然而，这些方法计算量大，需要大量的计算资源和存储空间。

两类方法的对比如图所示

比对方法可进一步分为剪接比对和连续比对。剪接比对能够将 reads 序列比对到参考序列的多个不同片段上，允许在参考序列的相应部分之间存在较大间隔。通常应用于将 RNA-seq 数据比对到基因组。

连续比对则寻求将 reads 与参考序列的连续子串进行比对，具有较好的比对率，通常不允许存在大间隔。

轻量级映射方法如kallisto、RapMap 和 Alevin-fry，通过避免核苷酸级别的比对，提高了速度和映射吞吐量。

这些方法基于一组匹配的 k-mer 或其他类型的精确匹配以及它们在 reads 和参考序列中的相对位置和方向来确定比对位置。这大大提高了速度和映射效率，但也使得难以进行基于得分的评估，以判断读取与参考序列之间的匹配质量。

不同参考序列

处理比对算法的选择之外，也可以对参考序列做出不同的选择。比如

• 完整的参考基因组
• 注释好的转录组
• 扩展的转录组

但并不是所有的比对算法都支持不同参考序列的选择。例如，基于轻量级映射的算法目前还不能根据参考基因组对 reads 进行剪切比对。

比对到基因组

这种方法将目标生物体的整个基因组作为参考，通常考虑转录本的注释信息。如 zUMIs、Cell Ranger 和 STARsolo 采用这种方法。

由于很多 reads 都是来自剪接后的转录本，因此这种方法需要使用能够跨越一个或多个剪接点进行比对的算法。

主要优点包括：

• 能够处理来自基因组任何位置的序列，包括未注释的转录本。
• 构建全基因组索引后，比对内含子和已知剪接转录本的 reads 效率是一样的。
• 能够处理注释转录本、外显子和内含子之外的序列，提高基因检测和定量的灵敏度。

缺点是：

• 比对工具需要大量内存
• 剪接比对的复杂性高，计算成本高。
• 需要目标生物体有完好的基因组，这对非模式生物来说可能是个问题。

比对到转录组

为了减少剪接比对的计算开销，广泛采用的一种替代方法是仅使用注释好的转录本序列作为参考。

由于大多数单细胞实验来自有良好注释转录组的模式生物（如小鼠或人类），这种基于转录组的定量方法可能提供与基于基因组方法类似的覆盖度。

优点在于：

• 与基因组相比，转录组序列大小较小，减少了运行比对管道所需的计算资源。
• 无需解决复杂的剪接比对问题，而只需搜索连续比对

但也有一些局限性：

• 不能捕捉来自转录组外的序列，这在处理 snRNA-seq 数据时尤其不足。
• 即使在 scRNA-seq 实验中，来自转录组外的 reads 也可能占据很大比例，不应该丢弃。

比对到扩展的转录组

为了解决 reads 可能来自剪接转录本之外的问题，可以将剪接转录序列与其他预期产生 reads 的参考序列（如全长未剪接转录本或剪除的内含子序列）合并，对参考序列进行扩展。

扩展的转录组参考序列比全基因组所需的参考索引通常更小，同时保留了仅搜索连续读取比对的能力。

能够更快且占用更少内存，同时仍能处理来自剪接转录组之外的许多读取。

通过比对到扩展的参考序列，可以收集更多的位置，并在使用轻量级比对方法时，大大减少假比对。

Alevin-fry 算法文章中提出，这种方法即使在处理标准 scRNA-seq 数据时也非常有价值，并建议构建增强的 splici（转录本+内含子）参考序列进行比对和定量。

splici 参考序列使用剪接转录组序列以及对应于注释基因内含子合并区间的序列构建。每个参考序列都标注其剪接状态，并在比对过程中使用考虑剪接状态的方法进行 UMI 解析。

这种方法的主要优点包括：

• 允许使用轻量级比对方法，同时大大减少假比对。
• 能检测到来自剪接和未剪接来源的 reads，提高后续分析的灵敏度。
• 统一了单细胞、单核和 RNA 速率的预处理流程。

细胞条形码矫正

基于液滴的单细胞分离系统，如 10x Genomics 系统，已成为研究细胞异质性原因及后果的重要工具。

在这种分离系统中，每个捕获的细胞的 RNA 物质在一个水基液滴内提取，并与带条形码的磁珠一起封装。

这些磁珠通过独特的寡核苷酸（称为 cell barcode，CB）标记每个细胞的 RNA，这些条形码随后与从 RNA 逆转录的 cDNA 片段一起测序。

磁珠包含高多样性的 DNA 条形码，使得细胞的分子内容可以进行平行条形码标记，并在后续数据处理中将测序 reads 分配到各个细胞中。

并不是所有与参考序列比对的测序片段都会被用于定量和条形码校正，一个常用的过滤标准是比对方向。

具体来说，某些化学反应策略规定，比对的读取应该仅来自于特定方向的转录本。例如，在 10x Genomics 3’ Chromium 中，期望的是 reads 与转录本的正义链比对，但其实也存在反义链上的 reads。

因此，如果测序方法遵循这样的链特异性协议，用户可以选择忽略或过滤反向互补方向比对的 reads。

条形码的错误类型

单细胞分析的标签、测序和去重方法通常表现良好。但是在基于液滴的文库中观察到的细胞条形码数量可能显著不同于最初封装的细胞数量。

主要有以下几个错误来源：

• 双重/多重（Doublet/Multiplet）：一个条形码可能与两个或更多细胞相关联，导致细胞计数不足。
• 空液滴：一个液滴中可能没有封装细胞，但环境中的 RNA 被条形码标记并测序，导致细胞计数过多。
• 测序错误：在 PCR 扩增或测序过程中可能出现碱基错误，导致细胞计数的过多或不足。

许多计算工具使用多种诊断指标来过滤非真实或低质量的数据。例如，许多方法用于去除环境 RNA 污染、检测 Doublet 和基于核苷酸序列相似性校正条形码错误等。

细胞条形码的校正策略

1. 基于已知条形码列表的校正

某些化学反应（如 10x Chromium）从已知的潜在条形码序列库中提取 CBs。因此，任何样品中观察到的条形码集应是该已知列表的子集。

常用策略是假定与已知列表中的某个序列完全匹配的每个条形码都是正确的，并校正所有其他条形码以匹配已知列表中的条形码。

但也存在一个问题，当一个错误条形码在已知列表中有多个可能的校正结果时，会导致校正结果不准。

2. 基于拐点的方法：

如果潜在条形码集未知，或希望直接从观测的数据中进行校正，可以采用这种方法。这种方法基于观察到样品中真实或高质量条形码通常与最多的 reads 相关联。

通过构建条形码的累积频率图（通常会出现一个拐点），可以识别频繁出现的条形码并将其视为有效条形码列表，其余条形码根据这个列表进行校正。

3. 基于预期细胞计数的过滤和校正：

这种方法适用于条形码频率分布没有明显拐点或存在技术伪影的情况。

用户提供预期的细胞数量估计，然后根据频率排序条形码，获取接近预期细胞数量的频率，并将其附近的条形码视为有效条形码。

4. 强制细胞数量的过滤：

用户直接提供排序频率图中条形码被视为有效的索引。频率大于或等于所选索引处频率的所有条形码被视为有效，并根据这个列表进行校正。

UMI 解析

在细胞条形码校正之后，reads 要么被丢弃，要么被分配到校正后的 CB。接下来，我们希望在每个校正后的 CB 中量化每个基因的丰度。

由于扩增偏差，必须根据 UMI 进行去重，以评估样本真实的分子数量。去重的目的是识别在 PCR 扩增前每个细胞中每个原始分子的序列片段数。

这个过程的结果是为每个细胞中的每个基因分配一个分子计数，随后在下游分析中作为该基因的原始表达估算值。

UMI 解析的必要性

在理想情况下，正确的（未改变的）UMI 标记的 reads，可以唯一比对到一个共同的参考基因，并且UMI 与原始分子之间存在双射关系。

UMI 去重过程就是：相同 UMI 的 reads 是来自同一个分子的 PCR 重复。每个基因捕获和测序的分子数量就是比对到该基因上的不同 UMI 数量。

但实际遇到的问题使得上述简单规则不足以识别 UMI 的基因来源，需要开发更复杂的模型。

主要有两个挑战：

• UMI 错误：在 PCR 或测序过程中引入的 UMI 序列出现错误（如碱基错配等），可能会使估计的分子数量增多
• 多重比对：当一个 reads 或 UMI 比对到参考序列的多个位置或多个基因时，导致基因的分子估算不准

基于图的 UMI 解析

目前开发了许多方法来解决 UMI 解析的问题，我们将重点介绍一种基于 UMI 图的框架。

每个连通分量代表一个子问题，其中某些 UMI 的子集被折叠（即，可以认为是相同来源的）。

许多流行的 UMI 解析方法可以通过简单地修改图的精炼方式和在图上执行的折叠或解析过程来解释。

在单细胞数据的背景下，UMI 图是一个有向图，其节点集和边集定义了 UMI 之间的关系。reads 的等价关系基于它们的 UMI 和比对信息。

我们所说的 reads 等价是，当且仅当它们具有相同的 UMI 标签并比对到相同的参考位置。

定义 UMI 图后，可以应用许多不同的解析方法。一种解析方法可能是简单地找到连通分量、聚类图、贪婪算法等。

UMI 解析的最后一步是使用解析后的 UMI 图量化每个基因的丰度。

对于丢弃多基因等价类的方法，生成当前处理细胞中基因的分子计数向量，方法是计算每个基因标记的等价类数量。

在这种模型中，折叠的等价类到基因的分配是潜变量，去重的基因分子计数是主要参数。通常，将基因唯一等价类的证据用于概率性地分配多基因等价类。

上述 UMI 解析和定量过程通常对每个细胞（校正后的 CB）单独执行，以创建所有细胞中所有基因的完整计数矩阵。

然而，高通量单细胞样本中每个细胞信息相对稀少，限制了 UMI 解析，反过来又限制了基于模型的解决方案的有效性。

计数矩阵质量控制

在生成计数矩阵后，进行质量控制（QC）评估是非常重要的。质量控制通常包括多个不同的评估，以帮助评估测序数据的整体质量。

这些评估包括基本的全局指标，如整体比对率、每个细胞观察到的不同 UMI 的分布、UMI 去重率的分布、每个细胞检测到的基因数量的分布等。

此外，还存在一些工具，如 Loupe browser, alevinQC 或 kb_python report、可以帮助组织和可视化这些基本指标。

空胞检测

质量控制的第一个步骤是确定哪些细胞条形码对应于高可信度测序的细胞。

在基于液滴方法中，某些条形码与环境 RNA 而非捕获细胞的 RNA 相关联，这是因为液滴未能捕获到细胞。

这些空液滴仍然会产生测序的 reads，但这些 reads 的特征与正确捕获的细胞条形码相关的读取有明显不同。

有许多方法可以评估条形码是否可能对应于空液滴。一种简单的方法是检查条形码的累积频率图，其中条形码按与之关联的不同 UMI 的数量降序排序。该图通常包含一个拐点，可以作为区分正确捕获细胞和空液滴的可能点。

虽然拐点的方法直观且通常能估计出合理的阈值，但也有一些缺点。例如，并非所有累积直方图都显示出明显的拐点，并且设计能够稳健地自动检测这些拐点的算法是非常困难的。

最后，仅凭与细胞条形码关联的总 UMI 计数可能并不是判断条形码是否关联空或受损细胞的最佳信号。

因此，开发了许多专门用于检测空或受损液滴或一般认为低质量细胞的工具，如 DIEM、dropkick 和 miQC 等。

这些工具结合了多种细胞质量的度量，包括不同 UMI 的频率、检测到的基因数量和线粒体 RNA 的比例，通常通过应用统计模型将高质量细胞与假定的空液滴或受损细胞分类。

这意味着通常可以对细胞进行评分，并根据估计的后验概率选择最终的过滤结果，确定细胞不是空的或受损的。

虽然这些模型通常适用于 scRNA-seq 数据，但在 snRNA-seq 数据中可能需要应用多个额外的过滤或启发式方法，以获得稳健的过滤结果。

双胞检测

除了确定哪些细胞条形码对应于空液滴或受损细胞，还需要识别那些对应于双重体或多重体的细胞条形码。

当一个液滴捕获两个（双胞）或多个（多胞）细胞时，这些细胞条形码在表示的 reads 数量和 UMI 数量以及基因表达谱方面会产生偏斜的分布。

许多工具也已被开发出来，用于预测细胞条形码的双胞状态。如 Scds、Solo 和 DoubletDecon

一旦检测到双胞和多胞，可以在后续分析中移除或进行相应调整。

计数矩阵

当完成初步的原始数据处理和质量控制后，继续进行后续分析时要注意，细胞与基因的计数矩阵只能近似反映原始样本中的测序分子。

在原始数据处理流程的多个阶段，采用了一些启发式方法和简化措施来生成这个计数矩阵。

例如，数列比对并不完美，细胞条形码的校正也存在不足。特别是 UMI 的准确解析非常具有挑战性，与多重比对相关的 UMI 问题经常被忽略，同时多个引物位点，特别是在未剪接分子中，可能违反常假设的一分子对应一 UMI 的关系。

鉴于这些挑战和为应对这些挑战而采用的简化措施，如何最好地表示预处理数据以供后续工具使用仍然是一个活跃的研究领域。

例如，一些量化工具实现了可选的 EM 算法，该算法以概率方式分配与多个基因映射的读取相关联的 UMI。然而，这可能导致非整数计数矩阵，这可能不太符合某些下游工具的预期。

另一种选择是将 UMI 解析到基因组而不是单个基因，保留预处理输出中的多重比对信息

值得注意的是，类似的表示方式在转录本水平上已经被使用了十多年，用于 bulk RNA-seq 转录本定量工具的内部表示，并且这种转录本水平的表示方式也被提议用于全长 scRNA-seq 数据的聚类和降维分析以及 scRNA-seq 数据的差异表达分析。

在这种情况下，生成的计数矩阵的维度将不再是 $C\times G$（其中 𝐺 是量化注释中的基因数量），而是 $C\times E$（其中 $E$ 是给定样本中所有细胞的不同基因分组的数量，通常称为等价类标签）。

通过将这些信息传播到输出计数矩阵，可以避免应用启发式或概率 UMI 解析方法，从而避免在下游分析中使用的计数数据中的数据丢失或偏差。

当然，为了使用这些信息，下游分析方法必须期望以这种格式的信息。此外，这些下游方法通常必须有一种方法最终将这些计数解析到基因水平，因为缺乏对基因分组丰度的直观生物学可解释性。

尽管如此，这种表示方式在向下游分析工具传递更完整和准确的信息方面提供的好处是显著的，并且正在开发利用这种表示方式的工具（如 DifferentialRegulation）。

讨论

具体工具的选择在很大程度上取决于手头的任务和可用计算资源的限制。

如果要进行标准的单细胞分析，基于轻量级映射的方法是一个不错的选择，因为它们比现有的基于比对的工具更快（通常快得多）、更节省内存。

如果要进行 snRNA-seq 分析，levin-fry 尤其是一个有吸引力的选择，因为它可以节省内存，而且即使转录组参考序列扩大到包括未剪接的参考序列，它的索引仍然相对较小。

另一方面，如果需要比对到转录组之外的序列，或者 reads 的基因组比对位置对于相关下游工具或分析（例如，使用 sierra 等工具进行差异转录本分析）是必要的，则推荐使用基于比对的方法。

根据 Brüning 等人的研究，在基于比对的流程中，STARsolo 比 Cell Ranger 更受青睐，因为前者比后者快得多，所需的内存也更少，同时还能产生几乎相同的结果。