3.1 参考基因组及其注释

Genome Biology | 最细RNA-seq数据分析名词

大多数scRNA-seq实验是使用人或小鼠的组织、类器官或细胞培养物完成的。这些生物基因组的初稿是在大约20年前发表的，但组装和注释定期更新。有两种常用的组装文件来源：UCSC（命名规律为hg19、hg38、mm10等）和GRC（GRCh37、GRCh38、GRCm38）。UCSC 和 GRC 组装的主版本与主染色体匹配（例如，来自 hg38 的 chr1 = 来自 GRCh38 的 chr1），但在额外的重叠群contig和所谓的 ALT 位点上有所不同，它们在次要版本之间有变化（例如 GRCh38.p13）。更多信息可以在这里和这里获得。基因组组装通常以fasta文件的形式分发 - 一个简单的文本文件，包括序列名称和序列。
Donor供体和Accepter受体
基因组组装
基因组组装：基础知识与基本思路

基因组注释过程包括定义基因组（基因）的转录区域，以及注释具有外显子-内含子边界的精确转录本，并为新定义的特征分配类型 - 例如蛋白质编码、非编码等。下面的示例显示了一个由 5 个转录本组成的基因：3 个蛋白质编码（红色）和 2 个非编码（蓝色）。基因组注释通常以 GTF 或 GFF3 文件格式提供，它们按层次结构组织。每个基因都由一个唯一的基因ID定义;每个转录本都由唯一的转录本 ID 及其所属的基因定义。外显子、UTR 和编码序列依次分配给特定转录本。

参考基因组及基因注释
基因组序列注释
保姆级参考基因组及其注释下载方法（图文详解）

图 3.1：典型真核基因的转录本和内含子-外显子结构

人类和小鼠基因组注释的流行来源是 RefSeq、ENSEMBL 和 GENCODE。RefSeq是三者中最保守的，并且每个基因的注释转录本往往最少。RefSeq转录本ID以NM_或NR开头，例如NM_12345。ENSEMBL 和 GENCODE 彼此非常相似，可以互换使用。这些基因名称以 ENSG（人类）和 ENSMUSG（小鼠）开头;转录本分别以 ENST 和 ENSMUST 开头。

除了基因 ID 之外，大多数基因还分配了一个通用名称（“gene symbol”）;例如，人肌动蛋白 B 具有 ensembl 基因 ID ENSG00000075624 和symbol ACTB。人类基因名称由HGNC定期更新和定义。小鼠基因名称由类似的联盟MGI决定。

目前人类基因组的ENSEMBL/GENCODE注释包含大约60k个基因，其中20k个是蛋白质编码，有237k个转录本。大多数基因可以按类型粗略地分为蛋白质编码基因、长非编码RNA、短非编码RNA和假基因。在更高的分辨率下，定义了 40 多种生物型。基因生物型注释也经常在注释版本之间发生变化（见下文）。
Gene and transcript types (ensembl.org)

如何选择合适的转录本进行过表达研究

NCBI数据库一个基因有多个转录本时选哪个？

图 3.2：GENCODE 版本中基因类型变化的桑基图

3.2 bulkRNA-seq和全长scRNA-seq数据的处理

批量RNA-seq的原始reads处理通常分两步完成： read alignment和read counting。这两个步骤都包含重要的警告，这些警告可能会强烈影响单个基因的表达估计。可以针对参考基因组或转录组进行read alignment。由于动物基因组中的广泛剪接，必须使用剪接感知比对器对基因组进行read比对;两种最流行的现代工具是 STAR 和 hisat2。典型的read coverage显示在下面的 A 中;请注意，给定基因的 3' 和 5' 末端的read coverage相对均匀。有些reads与超过 1 个位置完美对齐;这些读取通常称为multimappers. 多映射比对。与转录组比对时，歧义要大得多，因为许多转录本彼此非常相似（例如，仅相差一个外显子）;然而即使在基因水平上，模糊性也是显而易见的（下图B）。

Quality control | Functional genomics II (ebi.ac.uk)

Read mapping or alignment | Functional genomics II (ebi.ac.uk)

转录组入门(6)： reads计数

序列回帖与multi-mapped reads的处理

图 3.3：A：批量 RNA-seq 或 Smart-seq2 实验中的典型read coverage;B：将reads分配给特征时出现的歧义类型

在与基因组或转录组比对后，可以在基因或转录本水平上总结读取计数。在基因组比对的情况下，最简单的策略是只计算映射到唯一位置（非多映射）的reads，并且只重叠overlapping 一个基因。然而，这不可避免地会在基因表达估计中产生偏差（Pachter，2011）。更高级的策略包括在read比对到的特征之间拆分counts（例如，如果一个read与 3 个同源基因比对程度同样好，则每个获得 1/3 的计数）。链特异性RNA-seq可减少位于相反链上的特征重叠情况下reads分配的模糊性。一个有效实现上述所有计数方法的程序是 Subread 包中的 featureCounts。

链特异性的RNA-Seq

当使用转录组比对时，读取分配的歧义对于简单的计数来说太大。因此，通过采用期望最大化EM 算法的最大似然丰度估计来计算每个转录本和每个基因丰度。这种方法会将不同比例的reads分配给它映射到的feature，并大大减少与多映射reads相关的偏差。然后，在基因水平上总结分配给转录本的reads比例。最广泛使用和支持最充分的程序是 RSEM。一般来说，这是批量RNA-seq定量的最准确方法（Pachter，2011）。

上述传统方法（比对，然后read定量）的替代方法是基于所谓的 pseudo-alignment伪比对方法。两种流行的工具，kallisto 和 salmon，使用非常相似的方法：

• 将参考转录组拆分为 k-mers 并制作 De Bruijn 图;
• 将RNA-seq reads转换为k-mers;
• 使用 k-mers 将reads分配给一个或多个转录本（“equivalence class”）;
• 在转录本或基因水平上总结计数。

期望最大化算法用于查找映射到多个转录本的reads的最佳分布。这两种工具都非常节省内存和 CPU，并且非常准确，特别是对于双端或长单端reads。伪比对不会生成比对 BAM 文件，因此如果需要可视化（例如，当使用 RNA-seq 进行转录本注释时），也应单独进行比对。

图 3.4：转录组 De Bruijn 图的伪比对

关于批量RNA-seq定量，应注意以下几点。首先，通常假设测序的cDNA片段的数量与细胞中存在的RNA量成正比。因此，当使用双端reads时，每个read pair只计数一次，因为它源于同一个cDNA片段。对于人类和小鼠等注释良好的基因组，使用单端reads进行RNA-seq是很常见的。其次，PCR重复在批量RNA-seq中通常被忽略，UMI的使用也不会带来实质性的好处。几项独立研究表明，重复去除或使用UMI不会显著提高批量RNA-seq的统计功效。

最后，虽然许多差异表达方法使用原始reads计数，但在进行聚类、PCA 和其他类型的探索性分析时，通常使用样本内归一化。这种归一化最流行的方法是将原始reads转换为每百万转录本 （TPM）。这种转换解释了两个偏差：1）不同样本有不同的测序深度，与基因表达差异没有直接关系;2）长基因预计比短基因产生更多的cDNA片段。因此，对于 TPM 计算，首先将原始读取计数除以有效转录本长度，有效转录本长度定义为转录本长度 - cDNA 片段大小 + 1。在此之后，生成的数字将线性缩放，加起来为 100 万。因此，特定样本的所有 TPM 值的总和始终等于（大约）1,000,000。

3.3 基于液滴的scRNA-seq数据的比对和定量

3.3.1 一般注意事项

单细胞RNA-seq数据在许多方面与批量RNA测序不同（参见上文单细胞RNA-Seq简介一章）。大多数现代scRNA-seq技术生成的read序列包含三个关键信息：

• 识别RNA转录本的cDNA片段;
• 细胞条形码 （CB），用于识别表达 RNA 的细胞;
• 唯一分子标识符 （UMI），允许折叠属于 PCR 重复的read。

与批量RNA-seq相比，scRNA-seq处理的RNA量要少得多，并且执行更多的PCR循环。因此，UMI 条形码变得非常有用，现在在 scRNAseq 中被广泛接受。文库测序通常采用双端reads，一个读段含有CB + UMI（在10x Chromium中read 1），另一个包含实际转录序列（在10x Chromium为read 2）。

经典的scRNA-seq工作流程包括四个主要步骤：

• 将cDNA片段映射到参考基因组;
• 将reads分配给基因;
• 为细胞分配reads（cell barcode demultiplexing）;
• 计算唯一 RNA 分子的数量（UMI deduplication）。

该过程的结果是基因/细胞计数矩阵，该矩阵用作每个细胞中每个基因的RNA 分子数量的估计。

3.3.2 Cell Ranger中的读取映射

Cell Ranger 是处理 10x Genomics Chromium scRNAseq 数据的默认工具。它使用 STAR 比对器，可对基因组的reads进行剪接感知比对。在此之后，它使用转录本注释 GTF 将reads分为外显子、内含子和基因间隔区，以及reads是否与基因组对齐。如果read至少有 50% 与外显子相交，则为外显子，如果为非外显子且与内含子相交，则为内含子，否则为基因间隔区（如下所示）。在read类型分配之后，进行map质量调整：对于与单个外显子位点对齐但也与 1 个或多个非外显子位点对齐的read，外显子位点被优先考虑，并且认为确实映射到外显子位点，给出最大的映射质量分数。

图 3.5.. Cell Ranger 中比对read的分类

Cell Ranger通过检查外显子和内含子与转录组的相容性，进一步将外显子和内含子可靠映射的读段映射到带注释的转录本。读段的分类基于它们是正义的还是反义的，以及它们是外显子的、内含子的，或者它们的剪接模式是否与与该基因相关的转录本注释兼容。默认情况下，转录组学的读段（上图中的蓝色）将转入 UMI 计数。

当检测的输入由细胞核组成时，很大一部分reads来自未剪接的转录本并与内含子对齐。为了计算这些内含子读数，“cellranger count”和“cellranger multi”流程可以使用选项 include-introns运行。如果使用此选项，则任何在正义方向上映射到单个基因的读段（包括上图中标记为转录组（蓝色）、外显子（浅蓝色）和内含子（红色）的read）都将被转移到 UMI 计数。include-introns 选项消除了对自定义“pre-mRNA”参考基因组的需要，该参考基因组将整个基因体定义为外显子。 重要的是，如果read仅与单个基因兼容，则认为该read具有唯一映射性。只有唯一的映射read才会被转移到 UMI 计数;多映射read被丢弃。在 HTML摘要 输出中，转移到 UMI 计数的reads集合称为“映射到转录组的可信reads”。

3.3.3 Cell Ranger 参考基因组准备

在深入研究引用处理的细节之前，请务必注意默认的人类和鼠标引用是如何准备的。原代基因组组装版本（即没有ALT基因座）用于所有版本的比对。注释 GTF 文件使用可在此处找到的脚本进行过滤。保留了以下生物型：蛋白质编码、长链非编码RNA、反义和属于BCR/TCR（即V/D/J）基因的所有生物型（请注意，较旧的参考版本不包括后者）。去除所有假基因和小的非编码RNA。Cell Ranger``Cell Ranger

Cell Ranger软件预打包了多个版本的参考基因组;2020-A是迄今为止的最新版本。下面列出了Cell Ranger以前使用的所有单独的组装和注释组合。使用每个参考基因组生成的未过滤的scRNAseq表达矩阵行数应与“过滤后基因”列中的值相等。此外，Cell Ranger还含有人+小鼠联合参考基因组，可用于涉及人类和小鼠细胞的实验。

3.3.4 Chromium 版本和CB白名单

细胞条形码序列是附着在识别单个细胞的珠子上的合成序列。唯一序列库称为白名单，取决于 Chromium 库制备试剂盒版本。白名单文件可从 Cell Ranger 存储库获得。Chromium 使用了三个白名单：737K-april-2014_rc.txt、37K-august-2016.txt、3M-february-2018.txt 。第一个 CB 长度为 14 个bp，另外两个长度为 16 个bp。下表提供了常用 10x 单细胞测序试剂盒的细胞条形码和 UMI 长度，以及相应的白名单文件：

Cell Ranger 使用以下算法根据白名单更正假定的条形码序列：

• 统计白名单上每个条码在数据集中观察到的频率;
• 对于在数据集中但不在白名单上的每个条形码，查找距离为 1 汉明距离的白名单序列。对于每个这样的序列：
  • 计算观察到的条形码源自白名单条码的后验概率，在不同的碱基处存在测序错误（基于碱基 Q 分数）;
  • 将观察到的条码替换为后验概率超过 0.975中 最高的白名单条码。

校正后的条形码用于所有下游分析和输出文件。在输出 BAM 文件中，原始未校正的条码用 CR 标签编码，更正后的条码序列用 CB 标签编码。无法分配更正条形码的读数将没有 CB 标签。 如果您想知道为什么 3M-february-2018.txt 文件实际上包含超过 600 万个独特的序列，可以在此处找到解释。

3.3.5 UMI计数

通常所说的“UMI计数”包括read counting，然后基于UMI序列的PCR重复折叠。在对 UMI 进行计数之前，Cell Ranger 会尝试纠正 UMI 序列中的测序错误。被可靠地映射到转录组的reads被放入共享相同条形码、UMI 和基因注释的组中。如果两组reads具有相同的条形码和基因，但它们的 UMI 相差一个碱基（即相距 1 的汉明距离），则其中一个 UMI 可能是由测序中的替换错误引入的。在这种情况下，受支持较少的读取组的 UMI 将更正为支持较强的 UMI。

Cell Ranger 再次按条形码、UMI（可能已校正）和基因注释对reads进行分组。如果两组或多组reads具有相同的条形码和 UMI，但基因注释不同，则保留支持reads最多的基因注释用于 UMI 计数，而丢弃其他reads。如果最大读取支持率相同，则所有读取组都将被丢弃，因为无法可靠地分配基因。

在这两个过滤步骤之后，每个观察到的条形码、UMI、基因组合被记录为未过滤的特征-条形码（即基因-细胞）矩阵中的UMI计数。支持每个UMI count的reads数也记录在分子信息文件中。

3.3.6 细胞过滤

未过滤（“原始”）特征条形码矩阵包含许多列，这些列实际上是空液滴。由于技术噪音例如来自破碎细胞的环境RNA的存在，这些液滴中的基因表达计数不为零。然而，它们通常可以通过存在的 RNA 量与真正的细胞区分开来。Cell Ranger为这种细胞过滤实现了两种算法，我们将称之为“Cell Ranger 2.2”和“Cell Ranger 3.0”过滤。

原始算法 （Cell Ranger 2.2） 在“barcode count vs UMIs per barcode”图中识别了第一个“拐点”。Cell Ranger 3.0 引入了一种改进的细胞检出算法，该算法能够更好地识别低 RNA 含量细胞群，尤其是当低 RNA 含量细胞混入高 RNA 含量细胞群时。例如，肿瘤样本通常含有大肿瘤细胞与较小的肿瘤浸润淋巴细胞 （TIL） 混合，研究人员可能对 TIL 群体特别感兴趣。新算法基于EmptyDrops方法（Lun等人，2018）。

![[Pasted image 20231115093505.png]]
图 3.6 Cell Ranger 2.2 和 3.0 滤波算法识别的拐点图和空液滴截止值

3.3.7 基于伪对齐的方法

伪比对（见上文）也可用于快速定量scRNA-seq数据集。目前，有两个软件套件实现了这种方法：kallisto/BUStools 和 Salmon/Alevin/Alevin-fry。为了保留模块化方法，两个生态系统都引入了自己的格式，允许存储定量结果：kallisto/UStools引入了BUS（CB，UMI和Set）文件格式（Melsted等人，2019），而Alevin/Alevin-fry使用RAD格式（Srivastava等人，2019）。

kallisto/BUStools和Alevin/Alevin-fry都实现了对CB和UMI的error correction和demultiplexing 算法 - 例如，Alevin不需要（但可以使用）CB白名单。然而，与基于比对的方法的最大区别在于伪比对的精度较低，并且包含multimapping reads。

kallisto/BUStools支持多种测序技术，包括 CEL-seq、CEL-seq2 和 SMART-seq 等低通量测序技术。支持实验的完整列表可以用 kb --list打印。Alevin 目前仅支持两种最流行的基于液滴的单细胞方案，即 Drop-seq 和 10x Chromium。

一般来说，kallisto/BUStools和Alevin非常高效，允许处理具有 2-4 Gb RAM 的人类或小鼠数据集，并且至少比 Cell Ranger快一个数量级。这两种工具还可以正确处理multimapping reads，从而减少受影响基因的定量偏差。然而，在一些出版物中已经注意到，基于伪比对的方法错误地将保留的内含子的reads映射到转录组（Melsted et al, 2021; Srivastava et al, 2020）。众所周知，scRNA-seq实验，特别是单核RNA-seq可以包含非常高比例的具有保留内含子的转录本。这种错误的分配使数百个未表达的基因看起来表达较弱，这可能会对下游分析产生重大影响，尤其是marker选择（Kaminow 等人，2021 年）。因此，仍然需要开发至少与 Cell Ranger一样准确且快得多的方法。

3.4 STARsolo 和 Alevin-full-decoy：高速和高精度

STARsolo是一个独立的管道，是本章前面提到的STARRNA-seq比对器的一部分。它的开发目标是生成与 Cell Ranger非常相似的结果，同时保持计算效率。通常，STARsolo比Cell Ranger在同一数据集上快几倍。 STARsolo使用的UMI 折叠、cell barcode demultiplexing和细胞过滤的方法有目的地重新实现Cell Ranger使用的算法。从版本 2.7.9a 开始，还能够正确定量多映射read，使其成为快速准确的 scRNA-seq 处理的选择（Kaminow 等人，2021 年）。另一个好处是它灵活地实现了cellular barcode 和 UMI search：知道read中的相对位置以及每个序列的长度，就可以处理大多数scRNA-seq方法生成的数据。`

Alevin开发人员还承认了上述 intronic reads 问题，开发了一种特殊的解决方案，涉及使用所谓的 decoy sequences。在最近的STARsolo预印本中，Alevin全基因组诱饵显示出与STARsolo 、Cell Ranger非常相似的准确性（Kaminow 等人，2021 年）。

3.5 关于非模式生物的注意事项

使用单细胞RNA-seq来表征鲜为人知的多细胞生物正变得越来越流行，特别是作为关键物种从头基因组组装项目的一部分。这里有两点需要注意。首先，正确组装和良好注释的线粒体序列至关重要，因为线粒体reads构成了每个scRNA-seq文库的很大一部分，并广泛用于实验质量控制（参见scRNA-seq简介）。最近的一项努力为许多非模式脊椎动物精心组装的线粒体序列汇编（Formenti 等人，2021 年）。MITOS2 是一种专门的服务器，可用于自动生成后生动物的高质量线粒体注释。

其次，需要注意的是，大多数从头测序基因组的注释方法生成的基因模型不包括UTR序列。3'和5' scRNA-seq方法在基因两端的读长分布都存在严重偏差（下图）。因此，使用没有UTR序列的基因注释将极大地扭曲定量和分析的结果。

纯二代测序从头组装基因组

图 3.7：3' 和 5' 10x scRNA-seq 实验中的典型read coverage

3.6 简要总结和处理建议

Cell Ranger是 10x Genomics 提供的默认软件套件，它仍然是使用最广泛的read比对和定量工具。如果缺乏生物信息学方面的经验，或者有许多样品被Cell Ranger处理，请坚持使用。我们鼓励您使用最新版本及其附带的最新注释文件（参见上面的 3.3）。同时，STARsolo和Alevin-full_decoy提供出色的计算速度和正确的多映射处理，从而减少定量偏差，同时保持与Cell Ranger的高兼容性 .对于有终端工具经验的用户来说，它们可能是最佳选择。最后，如果你正在处理一个注释不良的基因组，请确保你的基因模型包括UTR，并且你有一个组装良好和注释清楚的线粒体。