【原】【数据分析】最好的单细胞RNA-seq分析实践教程（一）

GoDesign 2022-08-17 发布于北京

展开全文

单细胞RNA-seq分析使得基因表达分析达到了极高的分辨率。越来越多的人开始做单细胞分析，相关的分析工具也越来越多，然而在众多工具中找到合适的工具并建立适用于自己数据的数据分析方法也变得更加困难。2019年6月来自德国的学者Malte D. Leucken和Fabian J. Theis在Molecular Systems Biology上发表文章详述了单细胞RNA-seq分析的经典步骤，包括数据预处理（质量控制，归一化，数据校正，特征选择以及降维），细胞层面和基因层面的下游分析。本文通过独立比较不同算法给出最好的实践推荐。相关算法见：https://www.github.com/theislab/single-cell-tutorial。

——主要内容——

1 预处理与可视化

1.1 质量控制（本文内容）

1.2 归一化（本文内容）

1.3 数据校正与整合

1.4 特征选择，降维与可视化

1.5 预处理数据的阶段

1.6 聚类分析

2 分析平台介绍

据统计，截止到2019年7月，共有385个单细胞RNA-seq分析工具。工具持续不断的改善促进了新的科学发现，然而也使得标准化变得复杂。除了工具繁杂外，编程语言不统一也是制约标准化的一大因素，大部分工具是用R语言和Python写的。比较流行的Seurat, Scater和Scanpy平台提供了整合的分析环境和大量的分析工具。基于这些限制因素，本文主要给出目前最常用也是最好的分析工具，而不只是介绍标准化分析流程，同时不考虑编程语言。

图1. 典型的单细胞RNA-seq分析工作流程。

——预处理与可视化——

Raw data经过处理后得到read counts矩阵或counts矩阵，软件包括Cell Ranger, indrops, SEQC，以及read质控的zUMIs。以下所说的data均指count矩阵。

质量控制（Quality Control， QC）

在分析单细胞基因表达数据之前，我们需要保证所有的细胞标记数据均来自于活细胞。细胞QC主要基于三个协变量：每个barcode的counts数量（count深度），每个barcode的基因数量，每个barcode的线粒体基因数量。之后通过给定阈值过滤掉异常峰。这些异常barcodes可能是死细胞，即细胞膜膜破裂的细胞，也可能是成对的细胞。比如说，count深度很低的barcodes可被检测到的基因数量少，而线粒体基因数量却多，这是因为细胞质mRNA通过破裂的膜渗出了，只有线粒体内的mRNA还存在。相反的，细胞如果counts数量太高，那么有可能是两个细胞。因此，高count深度阈值主要用于过滤潜在的成对细胞。主流的成对细胞过滤方法包括DePasquale， Scrublet， Doublet Finder。

但是如果仅考虑这三个协变量，很有可能导致细胞信号的错误解释。比如，线粒体counts比例高的细胞可能在进行有氧呼吸过程。同样的，其他的QC过程也可能有其他的解释。低counts数的细胞也可能在静息状态，高counts数的细胞本身体积可能比较大。事实上，细胞间的molecular counts可能区别更明显。因此，过滤条件及阈值应综合考虑多个因素。未来，考虑多变量QC依赖性的过滤模型可以提供更灵敏的QC选择。

除了验证细胞的完整性，QC步骤也需要分析转录本。Raw count matrices包括超过20,000个基因，通过过滤在大部分细胞中不表达的基因，这一数量将大大降低。这一阈值通常用我们关心的细胞的最小聚类数量来定义，并留出一些dropout效应的余地。比如，过滤掉在不超过20个细胞中表达的基因可能会使得低于20个细胞的聚类类别无法被检测到。对于高dropout rates的数据集，这一阈值会使得检测大的聚类类别变得更复杂。这一阈值的选择需要考虑数据集中细胞的范围以及方便下游分析。

进一步的QC分析可以直接在count数据上进行。周围环境基因表达（Ambient gene expression）只非来源于barcode细胞，而是来源于其他裂解细胞的mRNA，也就是在构建细胞文库时由于细胞悬液被污染而来的mRNA。这些来自周围环境的counts对于后续下游分析有很大影响，包括标记基因识别或其他差异表达检测等。这一问题可以通过空载的单细胞液滴表达量数据来校正。SoupX工具可以用来进行这一步校正。在后续分析中直接忽略表达较强的环境基因也可以解决这一问题。

QC主要确保数据的质量足够用于后续分析。“足够的数据质量”无法先验确定，只能够通过后续分析来判断（比如cluster注释）。因此，在分析数据时，有必要重复修正QC的决策。一般在最开始分析时，先采用较为宽松的QC阈值，观察对于后续结果的影响，再进行更为严格的QC。尤其是有多个不同细胞群落时更需要注意，因为有些细胞类型有可能会被误认为低表达量的异常细胞。值得注意的是，QC的阈值不应该以提高统计分析的结果而优化，而是应该通过数据集中QC协变量的可视化和聚类来评估。

图2. 质量控制指标图示，以对于小鼠肠上皮数据集的预处理过滤决策为例。

归一化

在count矩阵中的每一个count都表示一个有效的转录本。Count深度对于同样的细胞来说可能会因为实验过程中其内部变化而有所不同。因此，当我们对比细胞间的count数据时，其差异可能仅仅由于采样效应。归一化则可以解决这一问题，也就是通过归一化数据得到细胞间的矫正后的基因表达数据。

Bulk基因表达有很多归一化方法。虽然其中的一些方法可以被用于scRNA-seq分析，但是单细胞数据的特异性差异来源，比如dropout（由于采样不均导致的counts为0），则促进了单细胞数据归一化的发展。

最常用的归一化方法是count深度归一化，也就是”counts permillion”或者CPM归一化。这个方法假设数据集中的每个细胞在最开始时有相同数量的mRNA分子，count深度差异仅由于采样导致。这一假设与降采样方法类似，也就是随机对reads或counts进行采样从而达到之前设定的采样数量，甚至少于这一数量。虽然降采样扔掉了一些数据，但是它提高了dropout rates，但是CPM及其他归一化方法则没有提高dropout rates。因此，降采样方法可以提供一种更加真实地细胞表达谱的表示。

由于单细胞数据集中经常包括不同的细胞群体，其细胞数量与分子counts均不同，一些更为复杂的归一化方法更适合于处理这一问题。包括high-count filtering CPM和Scran方法。这些方法都是线性的，全局的归一化方法。

非线性方法可以考虑到一些更为复杂的多于变化。大部分这些方法涉及对count数据的参数化建模。比如说，有研究利用负二项分布对count数据进行建模，利用包括read深度和每个基因的counts数量这样的协变量来拟合参数。模型拟合的残差表示了基因表达量的归一化。这样的方法可以结合技术上与生物上的数据校正。非线性归一化比全局归一化方法的效果要好，尤其是在batch effects较强的时候。因此，非线性归一化主要应用于基于平板的scRNA-seq数据，平板间会存在batch effects。另外，平板数据在count深度上的变化会比液滴数据要大。

我们不能指望有适用于所有scRNA-seq数据的归一化方法。比如全长scRNA-seq数据更适用于考虑基因长度的归一化方法，而3’ enrichment数据则不需要。全长数据中比较常用的方法是TPM归一化方法，这种方法起源于bulk数据分析。

就像细胞的count数据通过归一化使他们在细胞间可比一样，基因的count数据也可以通过归一化使它们在不同基因间对比。基因归一化包括将基因counts缩放为均值为0方差为1（z scores）。这一变化保证了所有基因在下游分析中权值相同。目前对于是否对基因进行归一化还没有定论。但是最流行的Seurat教程通常应用基因变换，Slingshout方法的作者则不进行这种变换。是否进行变换的关键点在于后续分析时，基因是否有必要变得权值相同，又或者这个基因表达的量级是否是一个重要的信息。为了获得更多的信息，本教程倾向于不对基因做变换。

经过归一化后，数据矩阵通常变为log(x+1)。这个变换有三个重要的影响。第一，log变换的表达值的距离通常表示log fold change；第二，log变换改变了单细胞数据的均值-方差关系；第三，log变换降低了数据偏差，从而接近数据正态分布的许多下游分析工具的假设。虽然scRNA-seq数据实际上并不是对数正态分布，但这三种影响使得log变换成为粗略的却有用的工具。这一用处主要体现在后续的差异性表达分析或是batch校正。不过也需要注意，归一化数据的log变换可能会向数据引入伪差异性表达效应。归一化大小因子在不同测试集中分布差异较大时，这一效应尤其明显。