单细胞RNA-seq预处理工具比较分析（bioRxiv）

前言

NGS系列文章包括NGS基础、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step) - Limma差异分析、火山图、功能富集）等内容。

2020年2月2日，苏黎世大学Mark Robinson团队在预印本bioRxiv上发表题为pipeComp, a general framework for the evaluation of computational pipelines, reveals performant single-cell RNA-seq preprocessing tools的研究内容，并建立了一个评估框架pipeComp，评估包括过滤、doublet检测、标准化、特征选择、降噪、降维和聚类等方面。

写在前面（结论）

由于该篇文章较长，所以先把结论放在这里，如果想去看下面作者的具体分析，可以继续往下看哦！作者在不同方面均提出较为具体的建议：

• （1）Filtering：建议使用Doublet检测和删除，并且可以使用诸如scDblFinder或scds之类的软件进行检测；基于细胞分布的过滤是无法捕获doublets的，应使用相对宽松的截止值（例如5个MAD或至少2个分布中的3个MAD）来排除质量较差的细胞；基于特征基因类型进行的过滤似乎没有益处。
• （2）Normalization and scaling：测试的大多数标准化方法均具有良好的性能，尤其是与归一化结合使用时，往往会对聚类产生积极影响；就亚群的可分离性而言，sctransform提供了最佳的整体性能，并且消除了文库大小和检测效率的影响；回归细胞协变量的常见做法，例如线粒体reads的检测率或比例几乎总是产生负面影响，从而导致与协变量的相关性增加并且聚类精度降低。因此，反不建议这么做。
• （3）Feature selection：Deviance提供了用于特征选择的最佳基因排名；包含的特征基因数量越多分类结果越好，我们的数据集中的4000个特征基因较为合理。
• （4）Denoising/imputation：去噪似乎有利于识别10x数据集中的亚群，但不适用于Smart-seq数据集；我们发现特别是ALRA（经过先前的标准化），DrImpute（经过预先处理）和DCA提供了最佳性能。
• （5）PCA：我们建议使用基于组分权重的irlba-based PCA ，如在Seurat中；建议不要使用常用的elbow方法（因为它过于保守）和jackstraw方法（性能低下，计算成本高）来确定要包括的组分，并建议使用基于改进的泊松混合模型的全局最大似然方法 。
• （6）Clustering：我们发现基于Scran的Walktrap聚类表现出良好的性能，但是计算量却增加了，这对于大型数据集可能是令人望而却步的；在先验知识可以指导分辨率选择的情况下，Seurat可以提供手动控制分辨率，而在缺少此类知识的情况下，基于Scran的walktrap聚类可以提供合理的估计。

可以看出这篇的结论与重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程 （原理、代码和评述）还是有不少不一样的地方。对照着阅读，在分析时根据自己的数据实际选择。没有标准的参数，只需要合理的解释。

背景

单细胞RNA测序（scRNAseq）及分析方法正在快速发展，并且已经出现了560多种软件工具，其中一半专门用于数据处理，例如聚类、排序、降维或标准化。随着新测序技术的发展，以及报道的细胞、基因和细胞种群数量的增加，可用工具的数量也随之增加。由于数据处理是任何scRNAseq分析的关键步骤，影响下游分析和解释，因此对可用工具的评估至关重要。因此，这种benchmark study可以帮助我们更好的进行工具选择。

作者收集了已知细胞组成的真实数据集（表1），并使用各种评估指标以多层次的方式调查围绕核心scRNAseq pipeline的各种参数和变化的影响。除了以前使用的具有真实细胞标签的benchmark数据集，还使用muscat模拟了两个数据集，这两个数据集具有基于真实的10x人和小鼠数据的分层次亚群结构（更多细节见表1和图1）。

图1:Overview of the benchmark datasets used.

表1

作者使用Seurat pipeline作为起点来进行以下方面的综合研究：

1）doublet识别；2）细胞过滤；3）标准化；4）特征选择；5）降维；6）聚类。

作者比较了不同方法并探究了较为合适的参数设置。

评估框架pipeComp及其应用原理

pipeComp包将pipeline定义为至少一个的、在前一个pipeline的输出上连续执行的函数列表(图2A)。此外，可以为每个步骤设置可选的基准函数，以提供标准化的、多层的评估指标。给定这样一个PipelineDefinition对象、一组可选参数(其中可能包括不同的子例程)和基准数据集，runPipeline函数将遍历所有参数组合，避免重复计算相同的步骤并动态编译计算(包括运行时)。可以使用pipeline中从该点向下的所有指标来评估给定参数的变化。这一点特别重要，因为end-point指标，例如用于聚类的adjusted Rand index(ARI)，并不完美。例如，虽然ARI评分的意义与真正的亚群的数量无关，但是到目前为止，所调用的聚类的数量是评分的最重要的决定因素：它离实际的亚群的数量越远，ARI越差。

图2: Overview of the pipeComp framework and its application to a scRNAseq clustering pipeline.

结果

（1）Doublet detection

Doublets被定义为在相同细胞barcode下测序的两个细胞（例如被捕获在同一液滴中），在scRNAseq数据集中相当常见，估计范围为1％至10％，具体取决于所使用的测序平台和细胞浓度（单细胞预测Doublets软件包汇总-过渡态细胞是真的吗？）。尽管相同细胞类型的doublets由于保留了基因之间的相对表达而在大多数下游分析中是无害的，但由不同细胞类型或状态形成的doublets可能会被错误分类，并可能扭曲下游分析。 在某些情况下，可以通过其异常高的reads和基因数来识别doublets，但这并非总是如此，于是出现了大量的用于识别doublets的工具。

作者首先使用10X数据集评估了这些方法检测doublets的能力，值得注意的是，基于SNP的分析通常将由相同细胞类型（但来自不同个体）的细胞创建的doublets描述为同型（与新型或异型doublets相反，即来自不同细胞类型的doublets）。这些同型双联体可能无法从单纯的基因计数中识别出来，通常被认为是无害的。尽管大多数方法可以准确地识别3个细胞系数据集（mixology10x3cl）中的doublets，但其他两个数据集却被证明更加困难（图3A）。scDblFinder是最快的方法（图3B），并且准确度也很高。在整个数据集中，被称为doublet的细胞更容易被分到其他的亚群中而不是单独聚为一类（图3C）。作者还发现，即使在设计数据不包含异型doublet的情况下，scDblFinder仍可以提高所有数据集聚类的准确性（图4）。

图3: Identification of doublet cells.

（2）Excluding more cells is not necessarily better

除了doublets以外，数据集可能还包括质量较低的细胞，消除这些细胞可以降低噪声。例如，对于含有高含量线粒体reads的液滴而言，这已被证明，通常是由于细胞降解和细胞质mRNA损失所致。通常的做法是基于某些属性，排除与大多数细胞不同的细胞。例如，可以通过scater的isOutlier函数执行，该函数针对给定的质控条件，测量每个细胞相对于所有细胞中位数的中值绝对偏差（MAD）。

然后作者研究了不同过滤条件的影响。根据细胞在某些特性（例如文库大小）上与整个细胞分布的距离来排除细胞，而这些特性倾向于在亚群中具有不同的分布。结果，整个分布中MAD数量的阈值可能导致对某些亚群的强烈偏差（图4A）。因此，作者研究了过滤准确性的提高与每个亚群所排除的最大细胞比例之间的权衡（图4B）。

由于过滤会更改不同亚群的相对丰度，因此不适合使用ARI等全局聚类准确性度量。因此，作者使用Hungarian algorithm计算了每个亚群的精确度和召回率，并监控了他们的平均F1得分（F1得分：统计学中用来衡量二分类模型精确度的一种指标）。最初的观察结果是，尽管严格过滤往往会提高准确性，但它也可能变得有害，并且在没有非常严格的过滤并使亚群偏差最小化的情况下也可以实现大多数收益（图4B）。在单个的细胞聚类上应用相同的过滤标准（通过scran的quickCluster方法确定），几乎没有细胞被过滤掉。表明对总细胞群的过滤倾向于丢弃具有更极端特性（例如高文库大小）的亚群细胞，而不是低质量的细胞。最后，通过改变分布，doublet去除步骤结合过滤有时会导致被排除细胞的比例净减少，同时保持或提高准确性。作者建议使用doublet去除，然后再进行相对温和的过滤，例如在“default”条件下进行过滤。

图4: Effect of filtering on cell subpopulation structure and clustering.

（3）Filtering features by type

线粒体reads与细胞降解有关，并且有证据表明核糖体基因可以影响聚类，在分析中隐藏其他生物学结构。作者发现核糖体基因的去除强烈地降低了聚类的质量，表明它们代表了亚群之间的真正生物学差异。删除线粒体基因并且只用蛋白质编码基因对聚类的影响非常小。

（4）Normalization and scaling

作者调查不同标准化策略的影响。除了Seurat中包含的标准对数标准化以外，还测试了scran基于pooling的标准化、sctransform's variance-stabilizing transformation、基于稳定基因的标准化以及SCnorm。

作者评估了第一主成分在何种程度上仍与文库大小和检测到的基因数量保持关联，从而消除了与细胞亚群的混杂协变量（图5A）。归一化趋向于消除与这些特征基因的大部分相关性，并且大多数方法都能够消除大部分影响，使用sctransform时相关性最低。然而让人意想不到的是回归协变量往往会增加与其的关联。

作者进一步评估了标准方法对亚群可分离性的影响（图5B-C）。由于诸如ARI之类的聚类准确性指标受聚类数量的强烈影响，因此作者用轮廓宽度(silhouette width)和互信息（mutual information，MI）对其进行了补充。作者发现大多数方法（包括完全不进行标准化处理）在大多数亚群中表现良好。归一化往往会减少某些亚群的平均轮廓宽度，并增加一些难以区分的亚群的轮廓宽度，通常（但并非总是）有利于最终聚类的准确性。回归协变量在所有指标上的表现都系统地较差。sctransform系统地胜过其他方法，即使将其开发为可用于具有唯一分子标识符（UMI）的数据，它在Smart-seq protocol（Koh和Kumar数据集）中也表现良好。

图5: Evaluation of normalization procedures.

（5）Feature selection and Dimensionality reduction

一个标准的聚类pipeline通常涉及到一个高变化基因的选择步骤，其复杂性在于（sc）RNAseq的数字性质和均值-方差关系。

Seurat的早期方法涉及到对均值表达水平进行标准化的离散估计，而更近的版本(≥3.0)则依赖于另一种标准化的方差度量。虽然均值-方差关系的调整消除了对高表达基因的许多偏差，但这种均值-方差关系实际上有时可能反映生物学相关性，并有可能有助于对细胞类型进行分类。另一种常见的特征选择方法是使用均值表达最高的特征。最近有人建议使用偏差，而sctransform则基于变换方差提供自己的基因排序。

假设选择方法应该理想地选择亚群之间的变化性比内部高的基因，作者首先评估每种方法选择的（真实）亚群解释具有高比例方差或差异的基因的程度。在选择由亚群解释的具有高比例差异的基因时，标准化的变异性度量在系统上比非标准化的同类差。然而，在解释偏差的百分比方面，标准化的方法往往更优(图6A)。Deviance被证明是亚群解释的高变化基因的优先选择方法，但在选择具有高deviance的基因方面表现不佳。

接下来，作者评估了不同特征基因选择方法如何影响聚类准确度（图6B）。总的来说，作者发现偏差和方差的非标准化估计可以在数据集和标准化方法之间提供最佳结果。增加选择的特征数量可以系统地提高聚类准确性，通常在4000个特征基因后达到平稳。

关于降维，作者将重点放在尚未进行比较的广泛使用的方法上：Seurat的PCA、scran的denoisePCA和GLM-PCA。鉴于Seurat的默认PCA通过每个组分的方差对细胞嵌入（cell embeddings）进行了加权，因此作者还评估了每种方法对该加权的影响。

选择降维方法的影响远大于归一化或特征选择的影响（图7）。GLM-PCA倾向于增加已经定义好的亚群的平均轮廓宽度，但是事实证明，Seurat的PCA程序在所有指标上都比较出色。与GLM-PCA一样，scVI的线性解码器（LD）也不明确地依赖归一化计数。对于具有更多亚群的数据集，在所有度量上，其性能均低于使用非加权Seurat PCA的pipeline。总体而言，按方差对主要成分进行加权（如Seurat）的方式对轮廓宽度和ARI得分具有积极影响。

图6: Evaluation of feature selection methods.

图7: Evaluation of common dimensionality reduction methods.

（6）Estimating the number of dimensions

降维的一个常见步骤是选择适当维度以用于下游分析。由于欧几里得距离随着非区分维数的增加而减小，因此通常需要选择足够的维数以保留大多数信息，排除可能的技术噪声。总体而言，维度数量的增加导致亚群数量的减少。目前已经提出了不同的方法来选择合适的维数，包括从视觉上解释方差的“Elbow”（拐点）到更复杂的算法。

作者评估了intrinsicDimension软件包中实现的维数估计的性能，以及诸如“elbow”法之类的通用过程、某些特定于scRNAseq的方法（如JackStraw过程）或scran的denoisePCA、和Fisher分离度分析的最新应用。根据Seurat的加权PCA，研究团队比较了各种维数估计在数据集中的能力。作为真实维数的第一近似值，作者计算了每个主成分的方差，该方差在大多数数据集中的前几个成分之后急剧减少（图8A）。

图8B显示了上述方法的维度估计与基于亚群的维度估计之间的差异（图8A）。值得注意的是，这些方法在计算时间方面存在很大差异（图8B），并且我们看不到估计的准确性和方法的复杂性之间的关系。大多数方法在各种聚类测度上均表现良好，但基于转换后的泊松混合模型（maxLikGlobal，使用20个nearest neighbors）的全局最大似然提供了维数估计，能够最好地分离亚群（图8C）并得到最佳的聚类准确性（图8D）。

图8: Estimating dimensionality

（7）Clustering

作者对不同的聚类方式进行比较。总体而言，与Seurat聚类相比，walktrap方法似乎优于fast greedy算法，并且通常不易出现错误分类，尽管后者提供了对分辨率的更多控制。基于Walktrap的聚类与Seurat之间的主要区别在于计算时间（图9B）。作者发现使用Annoy近似到nearest neighbors在某种程度上减少了计算时间，并且没有明显的准确度的损失。但是，它仍然比Seurat慢得多。对于所有方法，Zhengmix8eq和simMix1数据集的一些难以区分的亚群在所有指标方面的分类仍然非常不准确。

图9: Evaluation of clustering methods.

（8）Further extensions to the pipeline: imputation/ denoising

作者根据对分类的影响评估了各种imputation或降噪技术。由于初步分析表明，所有方法对标准化数据的性能均相同或更好，因此研究团队在将数据过滤和标准化之后、归一化和降维之前应用了这些方法。尽管某些方法（例如DRImpute process and alra norm）确实改善了一些难以捉摸的亚群的可分离性，但是没有一种方法对所有亚群的平均轮廓宽度有系统的积极影响（图10A）。

当限制自己的聚类分析产生”适当”的聚类数目时，所有经过测试的方法都比没有imputation步骤的情况改善了分类（‘none’ ，图10B）。某些方法（例如，ENHANCE）始终表现不佳。10X数据集的典型特征是较低的细胞覆盖率和特征基因检测率，DrImpute和DCA表现出最佳性能。相反，基于Smart-seq技术的标准化counts反而不利于聚类精度。

图10: Evaluation of imputation/denoising methods

作者：Tiger 编辑：生信宝典