参考基因组及其注释
大多数scRNA-seq实验是使用人类或小鼠组织、器官或细胞培养物进行的。尽管这些基因组的初稿是在20年前发表的，但组装和注释的更新是相当定期的。有两个流行的组装文件来源： UCSC（他们的汇编被命名为hg19、hg38、mm10等），和GRC（GRCh37、GRCh38、GRCm38）。UCSC和GRC汇编的主要版本在主要染色体上是匹配的（如hg38的chr1=GRCh38的chr1），但在额外的contigs和所谓的ALT位点上是不同的，这些位点在次要版本之间是变化的（如GRCh38.p13）。更多的信息可以在这里和这里获得。基因组组装通常以fasta文件的形式发布--一个简单的文本文件，包括序列名称和序列。
基因组注释过程包括定义基因组的转录区域（基因），以及用外显子-内含子边界注释确切的转录物，并给新定义的特征分配一个类型--如蛋白质编码、非编码等。下面的例子显示了一个由5个转录本组成的基因： 3个蛋白质编码（红色）和两个非编码（蓝色）。基因组注释通常以GTF或GFF3文件格式提供，这些文件是按层次组织的。每个基因由一个独特的基因ID定义；每个转录本由一个独特的转录本ID和它所属的基因定义。外显子、UTR和编码序列又被分配给一个特定的转录本。
人类和小鼠基因组注释的流行来源是RefSeq, ENSEMBL和GENCODE。RefSeq是三者中最保守的，而且每个基因的注释转录物往往最少。RefSeq转录本的ID以NM_或NR开头，如NM_12345。ENSEMBL和GENCODE非常相似，可以为我们的目的互换使用。这些中的基因名称以ENSG（人类）和ENSMUSG（小鼠）开始；转录本分别以ENST和ENSMUST开始。
除了基因ID，大多数基因还有一个通用名称（"基因符号"）；例如，人类肌动蛋白B将有ensembl基因ID ENSG00000075624和符号ACTB。人类基因名称由HGNC定期更新和定义。小鼠基因名称由一个类似的联盟MGI决定。
目前人类基因组的ENSEMBL/GENCODE注释包含大约6万个基因，其中2万个是蛋白编码基因，还有23万个转录本。大多数基因可以按类型粗略地划分为蛋白质编码基因、长非编码RNA、短非编码RNA和假基因。在更高的分辨率下，定义了40多个生物类型。基因生物型的注释也经常在不同的注释版本之间变化（见下文）。
处理bulk RNA-seq和全长scRNA-seq数据
Bulk RNA-seq的原始读数处理通常分两步进行：读数对齐和读数计数。这两个步骤都包含重要的注意事项，会对单个基因的表达量估计产生很大影响。读数比对可以对照基因组或转录组参考文献进行。由于动物基因组中广泛存在的拼接现象，针对基因组的读数比对必须用拼接识别器来完成；两个最流行的现代工具是STAR和hisat2。典型的读数覆盖率显示在下面的面板A；注意读数覆盖率在给定基因的3'和5'端都是相对均匀的。有些读数与1个以上的位置完全一致；这些读数通常被称为多映射器。当对准转录组时，模糊性要大得多，因为许多转录本彼此非常相似（例如只相差一个外显子）；然而，即使在基因水平上，这种模糊性也是明显的（下面的面板B）。
在与基因组或转录组进行比对后，可以在基因或转录物水平上对读数进行总结。在基因组比对的情况下，最简单的策略是只计算映射到一个独特位置的读数（非多重映射者），并且只重叠一个基因。然而，这不可避免地会在基因表达估计中产生偏差（Pachter，2011）。更高级的策略包括将读数在它所对准的特征之间进行分割（例如，如果一个读数与3个旁系基因的对准程度相同，则每个旁系基因得到⅓的计数）。特定链的RNA-seq减少了在位于相反链上的重叠特征的情况下读数分配的模糊性。Subread软件包中的featureCounts就是一个有效实现上述所有计数方法的程序实例。
当使用转录组排列时，读数分配的模糊性对于简单的计数来说太大。因此，每个转录本和每个基因的丰度是用最大似然丰度估计来计算的，使用的是期望-最大化（EM）算法。这种方法使读数的不同部分被分配给它所映射的特征，并大大减少了与多映射器有关的偏差。被分配到转录本的读数（和读数部分）然后在基因水平上进行总结。实施这一策略的最广泛使用和支持的程序是RSEM。一般来说，这是批量RNA-seq定量的最准确的方法（Pachter，2011）。
上述传统方法（先对齐，再读数定量）的另一种方法是基于所谓的伪对齐方法。两个流行的工具，kallisto和salmon，使用非常相似的方法：
将参考转录组分成k-mers，并制作一个De Bruijn图；
将RNA-seq读数转换成k-mers；
使用k-mers将读数分配给一个转录本或几个转录本（"等价类"）；
在转录本或基因水平上总结所产生的计数。
期望最大化算法被用来寻找映射到多个转录本的读数的最佳分布。这两种工具都非常节省内存和CPU，而且相当准确，特别是对成对端或长单端读数。伪比对不产生比对的BAM文件，所以如果需要可视化（例如，当使用RNA-seq进行转录物注释时），也应单独进行比对。
关于批量RNA-seq的定量，有几件事需要注意。首先，一般认为测序的cDNA片段的数量与细胞中存在的RNA数量成正比。因此，当使用成对的读数时，每个读数对只被计算一次，因为它来自同一个cDNA片段。对于像人类和小鼠这样被很好地注释的基因组来说，使用单端读数进行RNA-seq是非常普遍的。其次，PCR重复片段通常在批量RNA-seq中被忽略，而且使用UMI也不会带来实质性的好处。一些独立的研究表明，去除重复序列或使用UMI并不能明显地提高批量RNA-seq的统计能力。
最后，虽然许多差异表达方法使用原始读数，但在进行聚类、PCA和其他类型的探索性分析时，使用样本内归一化是常规做法。这种归一化的最流行的方法是将原始读数转换为每百万转录本（TPM）。这种转换考虑了两个偏差： 1）不同样品的测序深度不同，与基因表达差异没有直接关系；2）长基因预计会比短基因产生更多的cDNA片段。因此，对于TPM的计算，原始读数首先除以有效转录本长度，即转录本长度-cDNA片段大小+1。之后，所得数字按线性比例加起来为一百万。因此，一个特定样本的所有TPM值的总和总是等于（大约）1,000,000。
基于液滴的scRNA-seq数据的读数比对和定量分析
一般考虑因素
单细胞RNA-seq数据与大量RNA seq有许多不同之处（见上文单细胞RNA-Seq章节介绍）。大多数现代scRNA-seq技术产生的读取序列包含三个关键信息：
确定 RNA 转录本的 cDNA 片段；
识别RNA表达的细胞的细胞条形码（CB）；
独特的分子标识符（UMI），允许折叠PCR重复的读数。
与批量RNA-seq相比，scRNA-seq处理的RNA量要小得多，而且要进行更多的PCR循环。因此，UMI条形码变得非常有用，现在在scRNAseq中已被广泛接受。文库测序通常使用成对的读数，其中一个读数包含CB + UMI（10x Chromium中的读数1），另一个包含实际的转录物序列（10x Chromium中的读数2）。
一个经典的scRNA-seq工作流程包含四个主要步骤：
将cDNA片段映射到一个参照物上；
将读数分配给基因；
将读数分配给细胞（细胞条码解复用）；
计算独特的RNA分子的数量（UMI重复计算）。
这个过程的结果是一个基因/细胞计数矩阵，它被用来估计每个细胞中每个基因的RNA分子的数量。
Read 使用Cell Ranger比对
Cell Ranger是处理10x Genomics Chromium scRNAseq数据的默认工具。它使用STAR比对，将 读数与基因组进行拼接感知的比对。之后，它使用转录本注释GTF将reads分成外音、内音和基因间，并通过reads是否与基因组对齐（有把握）来确定。如果一个reads至少有50%与外显子相交，那么它就是外显子；如果它是非外显子并与内含子相交，那么就是内含子；否则就是基因间（如下图）。在reads类型分配之后，对映射质量进行调整：对于与单个外显子基因座对齐但也与1个或多个非外显子基因座对齐的reads，外显子基因座被优先考虑，该reads被认为是有把握地映射到外显子基因座上，并给予最高的比对质量分数。
Cell Ranger通过检查它们与转录组的兼容性，进一步将外显子和内含子映射的reads与注释的转录本相匹配。读数的分类是基于它们是有义还是反义，基于它们是外音、内音还是它们的剪接模式与与该基因相关的转录本注释兼容。默认情况下，属于转录组的reads（上图中的蓝色）被转入UMI计数。
当检测的输入由细胞核组成时，很高比例的reads来自未拼接的转录本并与内含子对齐。为了计算这些内含子的reads，"cellranger count "和 "cellranger multi "管线可以在运行时加入include-introns选项。如果使用这个选项，任何以敏感性方向比对到单一基因的reads--包括上图中标有转录体（蓝色）、外显子（浅蓝色）和内含子（红色）的reads--都会被转入UMI计数。包括内含子的选项消除了对自定义 "pre-mRNA "参考的需要，该参考将整个基因体定义为一个外显子。重要的是，如果一个reads只与一个基因兼容，就被认为是唯一映射的。只有唯一映射的reads才会被转到UMI计数；多映射的reads会被Cell Ranger丢弃。在网络摘要的HTML输出中，继续进行UMI计数的reads集被称为 "自信地比对到转录组的reads"。
Cell Ranger的参考准备
在我们深入研究参考文献处理的细节之前，有必要注意一下默认的Cell Ranger人类和小鼠参考文献是如何准备的。初级基因组装配版本（即没有ALT位点）在所有的版本中都被用于对齐。注释的GTF文件被过滤，使用的脚本可以在这里找到。以下生物型被保留：蛋白质编码、长非编码RNA、反义以及属于BCR/TCR（即V/D/J）基因的所有生物型（注意，旧的Cell Ranger参考版本不包括后者）。所有假基因和小的非编码RNA都被删除。
有几个版本的Cell Ranger参考文献与软件预先打包在一起；2020-A是迄今为止最新的参考文献版本。Cell Ranger之前使用的所有单独的组装和注释组合都列在下面。使用每个参考文献生成的未经过滤的scRNAseq表达矩阵，预计包含的行数与 "过滤后的基因 "栏中的数值相等。此外，Cell Ranger还包含人类+小鼠的组合参考，这对同时涉及人类和小鼠细胞的实验很有用。
Chromium版本和单元格条码白名单
细胞条码序列是附着在珠子上的合成序列，可以识别单个细胞。独特的序列库被称为白名单，取决于Chromium库制备试剂盒的版本。白名单文件可以从Cell Ranger资源库中获得。Chromium使用的白名单有三个：737K-april-2014_rc.txt，737K-august-2016.txt，以及3M-february-2018.txt。第一个列表中的CB是14bp长，另外两个是16bp。下表提供了流行的10倍单细胞测序试剂盒的细胞条形码和UMI长度，以及适当的白名单文件：
| Chemistry | CB, bp | UMI, bp | Whitelist file |
| 10x Chromium Single Cell 3’ v1 | 14 | 10 | 737K-april-2014_rc.txt |
| 10x Chromium Single Cell 3’ v2 | 16 | 10 | 737K-august-2016.txt |
| 10x Chromium Single Cell 3’ v3 | 16 | 12 | 3M-february-2018.txt |
| 10x Chromium Single Cell 3’ v3.1 (Next GEM) | 16 | 12 | 3M-february-2018.txt |
| 10x Chromium Single Cell 5’ v1.1 | 16 | 10 | 737K-august-2016.txt |
| 10x Chromium Single Cell 5’ v2 (Next GEM) | 16 | 10 | 737K-august-2016.txt |
| 10x Chromium Single Cell Multiome | 16 | 12 | 737K-arc-v1.txt |
Cell Ranger使用以下算法来纠正白名单上的假定条码序列：
计算白名单中每个条形码在数据集中的观察频率；
对于数据集中不在白名单上的每一个观察到的条形码，找到与白名单相距1-Hamming距离的白名单序列。对于每个这样的序列：
计算观察到的条形码来源于白名单条形码的后验概率，在不同的碱基上有一个测序错误（基于碱基的质量得分）；
用具有超过0.975的最高后验概率的白名单条码替换观察到的条码。
更正后的条形码用于所有下游的分析和输出文件。在输出的BAM文件中，原始的未校正的条形码被编码在CR标签中，而校正后的条形码序列被编码在CB标签中。不能被分配到校正的条形码的读数将没有CB标签。如果你想知道为什么3M-february-2018.txt文件实际上包含600多万个独特的序列，可以在这里找到解释。
UMI计数（Counting）
通常所说的 "UMI计数 "包括读取计数，然后根据UMI序列进行PCR重复折叠。在计算UMI之前，Cell Ranger试图纠正UMI序列中的测序错误。有把握地映射到转录组的读数被放入共享相同的条形码、UMI和基因注释的组中。如果两组读数具有相同的条形码和基因，但它们的UMI相差一个碱基（即相差1个汉明距离），那么其中一个UMI可能是由测序中的替换错误引入的。在这种情况下，支持率较低的读组的UMI被纠正为支持率较高的UMI。
Cell Ranger再次按照条形码、UMI（可能经过校正）和基因注释对读数进行分组。如果两组或更多的读数有相同的条形码和UMI，但有不同的基因注释，那么有最多支持读数的基因注释将被保留用于UMI计数，而其他读数组则被丢弃。在最大支持读数并列的情况下，所有读数组都被丢弃，因为基因不能被有把握地分配。
在这两个过滤步骤之后，每个观察到的条形码、UMI、基因组合被记录为未过滤的特征-条形码（即基因-细胞）矩阵中的UMI计数。支持每个计数的UMI的读数也被记录在分子信息文件中。
细胞过滤
未经过滤的（"原始"）特征条码矩阵包含许多列，实际上是空的液滴。由于技术上的噪音，例如来自破碎细胞的环境RNA的存在，这些液滴中的基因表达计数不是零。然而，它们通常可以通过存在的RNA的数量与真正的细胞区分开来。在Cell Ranger中，有两种算法实现了这种细胞过滤，我们将其称为 "Cell Ranger 2.2 "和 "Cell Ranger 3.0 "过滤。
最初的算法（Cell Ranger 2.2）确定了 "条形码计数与每个条形码的UMI数 "图中的第一个 "膝点"。Cell Ranger 3.0引入了改进的细胞计数算法，能够更好地识别低RNA含量的细胞群，特别是当低RNA含量的细胞混入高RNA含量的细胞群时。例如，肿瘤样本通常包含大型肿瘤细胞与较小的肿瘤浸润淋巴细胞（TIL）混合，研究人员可能对TIL群体特别感兴趣。新算法是基于EmptyDrops方法（Lun等人，2018）。
基于伪对准的方法
伪比对（见上文）也可用于快速量化scRNA-seq数据集。目前，有两个软件套件实现了这种方法：kallisto/BUStools，以及Salmon/Alevin/Alevin-fry。为了保留模块化方法，这两个生态系统都引入了自己的格式，允许存储量化结果：kallisto/BUStools引入了BUS（Barcode, UMI, and Set）文件格式（Melsted et al, 2019），而Alevin/Alevin-fry正在使用RAD格式实现同样的目的（Srivastava et al, 2019）。
kallisto/BUStools和Alevin/Alevin-fry都实现了自己的细胞条形码和UMI纠错和解复用的算法--例如，Alevin不需要（但可以使用）细胞条形码白名单。然而，与基于比对的方法最大的区别在于伪比对的准确性较低，以及包含多映射读数。
kallisto/BUStools支持许多测序技术，包括CEL-seq、CEL-seq2和SMART-seq等低通量技术。支持的实验的完整列表可以用kb --list打印。Alevin目前只支持两个最流行的基于液滴的单细胞协议，即Drop-seq和10x Chromium。
总的来说，kallisto/BUStools和Alevin的效率非常高，可以用2-4Gb的内存处理人类或小鼠数据集，比Cell Ranger至少快一个数量级。这两种工具还能正确处理多映射读数，减少受影响基因的量化偏差。然而，一些出版物指出，基于伪比对的方法会将保留内含子的读数错误地映射到转录组中（Melsted等，2021；Srivastava等，2020）。众所周知，scRNA-seq实验，特别是单核RNA-seq，可以包含非常高比例的保留内含子的转录物。这种错误的分配使得数百个非表达基因看起来表达很弱，这可能会大大影响下游的分析，特别是标记物的选择（Kaminow等人，2021）。因此，目前仍在推动开发至少与Cell Ranger一样准确且速度更快的方法。
STARsolo和Alevin-full-decoy：高速度和高精确度
STARsolo是一个独立的管道，是本章之前提到的STAR RNA-seq对齐器的一部分。开发它的目的是为了产生与Cell Ranger非常相似的结果，同时保持计算效率。通常情况下，在相同的数据集上，STARsolo比Cell Ranger快几倍。STARsolo用于UMI折叠、细胞条形码解复用和细胞过滤的方法是有目的地重新实现了Cell Ranger使用的算法。从STAR 2.7.9a版本开始，STARsolo也能够正确量化多映射读数，使其成为快速准确的scRNA-seq处理的一个非常有吸引力的选择（Kaminow等人，2021）。STARsolo的额外好处是它灵活地实现了细胞条码和UMI搜索：知道读数内的相对位置，以及每个序列的长度，就可以处理大多数scRNA-seq方法产生的数据。
Alevin的开发者也承认上述内含子读数的问题，开发了一种特殊的解决方案，涉及到使用所谓的诱饵序列。在最近的STARsolo预印本中，带有全基因组诱饵的Alevin显示出与STARsolo或Cell Ranger非常相似的准确性（Kaminow等人，2021）。
关于非模型生物的说明
使用单细胞RNA-seq来描述不太为人所知的多细胞生物的特征正变得越来越流行，特别是作为关键物种的新基因组组装项目的一部分。这里有两件事需要注意。首先，一个正确组装和良好注释的线粒体序列是至关重要的，因为线粒体读数在每个scRNA-seq文库中占相当大的比例，并被广泛用于实验质量控制（见scRNA-seq简介）。最近的一项工作是将许多非模式脊椎动物的线粒体序列精心组装起来的汇编（Formenti等人，2021）。MITOS2是一个专门的服务器，可用于自动生成高质量的后生动物线粒体注释。
其次，重要的是要注意到，大多数从头测序的基因组的注释方法产生的基因模型不包括UTR序列。3'和5'scRNA-seq方法的读数分布都严重偏向基因的两端（如下图）。因此，使用没有UTR序列的基因注释将极大地扭曲量化和分析的结果。
Cell Ranger是10x Genomics提供的默认软件套件，它仍然是最广泛使用的读数对齐和量化工具。如果你缺乏生物信息学方面的经验，或者有很多其他用Cell Ranger处理的样品，请坚持使用它。我们鼓励你使用最新的Cell Ranger版本，以及它附带的最新注释文件（见上文3.3）。同时，STARsolo和Alevin-full_decoy提供了巨大的计算速度和正确的多映射处理，这减少了量化偏差，同时与Cell Ranger保持了非常高的兼容性。它们可能是对终端工具有经验的用户的最佳选择。最后，如果你正在处理一个注释不全的基因组，确保你的基因模型包括UTRs，而且你有一个组装好的、注释好的线粒体。
注：本文是翻译scRNA-seq官方文档，相关链接为：3 Processing Raw scRNA-Seq Sequencing Data: From Reads to a Count Matrix | Analysis of single cell RNA-seq data
2670
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
