单细胞RNA-seq数据分析最佳实践（中）
分享是一种态度
作者 | 周运来
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树，
一枚大型测序工厂的螺丝钉，
一个随机森林中提灯觅食的津门旅客。
文章信息
前文回顾
Pre-processing and visualization
Normalization
计数矩阵中的每个计数代表细胞 mRNA 分子的成功捕获、逆转录和测序（框 1）。由于每个步骤固有的变异性，相同细胞的计数深度结果却可能不同。因此，当基于计数数据比较细胞间的基因表达时，任何差异可能仅由采样效应（ sampling effects.）引起。通过例如缩放（ sampling effects）计数数据以获得正确的细胞间相对基因表达丰度来解决这一问题。
bulk RNA数据已有许多标准化方法 (preprint:Pachter,2011;Dillies et al,2013)。虽然其中一些方法已应用于 scRNA-seq 分析，但单细胞数据特有的变异来源如技术脱落（technical dropouts ）（取样导致的零计数，双零问题）促使开发出了针对 scRNA-seq 的标准化方法 (Lun et al,2016a；Vallejos et al,2017)。
最常用的规范化协议是 count depth scaling，也称为每百万计数或 CPM 规范化。该方案来自bulk 表达分析，并使用与每个细胞计数深度成比例的所谓大小因子对计数数据进行标准化。该方法的变体使用不同的因子或数据集中每个细胞的中位计数深度缩放。CPM 标准化假设数据集中的所有细胞最初包含相同数量的 mRNA 分子，计数深度差异仅由于取样产生。该假设与下采样（downsampling）方案相同，下采样方案是从数据中随机取样读取或计数，使所有细胞的计数预先规定的数量或更少。在下采样丢掉数据的同时，也增加了技术脱落率，而 CPM 和其他全局缩放标准化方法则没有。因此，下采样可以提供类似计数深度下细胞表达谱的更真实表示。
由于单细胞数据集通常由具有不同大小和分子计数的异质细胞群组成，因此更复杂的标准化方法通常是合适的。例如，Weinreb et al (2018) 使用了 CPM 的简单延伸，计算它们的大小因子时，排除了任何细胞中占总计数至少 5% 的基因。这种方法考虑到了少数高表达基因的分子计数变异性。基于 Scran 合并的尺寸因素估计（pooling-based size factor estimation）方法允许更多的细胞异质性 (Lun et al,2016a)。细胞合并后，根据基因的线性回归估算大小因子，以避免技术脱落效应。该方法将变异性限制在细胞间差异表达基因的 50% 以下，并且在独立比较中始终是性能最佳的标准化方法。经证明，Scran 在批次校正（ batch correction） (Buttner et al,2019) 和差异表达分析 (preprint：vith et al,2019) 方面的性能优于其他检测的归一化方法。在与原作者的小规模比较中，该方法也显示出稳健的尺寸因子估计值 (Vallejos et al,2017)。
CPM、高计数过滤 CPM 和 scran 使用线性、全局缩放标准化计数数据。还存在非线性归一化方法，可解释更复杂的异质性(Cole et al,2019)。许多方法涉及到计数资料的参数化建模。例如，Mayer et al (2018) 使用技术变量（如测序深度和每个基因的计数数量）拟合负二项模型，以拟合模型参数。模型拟合的残差作为基因表达的标准化定量。这种方法可以将技术和生物数据校正（例如批次校正或细胞周期效应校正）与计数深度归一化相结合。已证明非线性归一化方法优于全局缩放方法，尤其是在具有强批次效应的情况下 (Cole et al,2019)。因此，非线性归一化方法对于基于平板的 scRNA-seq 数据尤其相关，这些数据往往在平板之间存在批次效应。此外，与基于液滴的数据相比，基于平板的数据可显示每个细胞计数深度的较大变化 (Svensson et al,2017)。虽然非线性归一化方法或替代方法（例如下采样）似乎更适合这些条件，但需要进行比较研究来确认该假设。在本教程中，我们倾向于将标准化和数据校正（批次校正、噪声校正等）步骤分开，以强调数据的不同处理阶段（参见预处理数据部分的阶段）。因此，我们着重研究全局尺度归一化方法（global scaling normalization）。
我们不能期望单一的标准化方法适合所有类型的 scRNA-seq 数据。例如，vith et al (2017) 表明reads和计数数据可通过不同模型进行最佳拟合。事实上，Cole et al (2019) 发现没有一种归一化方法对不同的数据集表现都是最佳的，并认为应使用其 scone 工具为特定数据集选择适当的归一化方法。此外，scRNA-seq 技术可分为全长和 30 种富集方法 (Svensson et al,2017；Ziegenhain et al,2017)。来自全长方案的数据可能受益于考虑到基因长度的标准化方法（例如 Patel et al,2014；Kowalczyk et al,2015；Soneson）
细胞计数数据可以归一化，使细胞间具有可比性，同样，基因计数也可以按比例调整，以改善基因间的比较。基因归一化构成基因计数的标度，使其均值和单位变异(z值)为零。这种比例的变化影响了所有的基因下游分析的权重。是否对基因进行归一化目前尚无共识。虽然流行的Seurat教程(Butler et al, 2018)通常应用基因缩放（scaling,），但Slingshot方法的作者在他们的教程中选择不缩放基因(Street et al, 2018)。这两种选择之间的偏好围绕着是否所有的基因都应该在下游分析中得到同等的权重，或者一个基因的表达量是否代表了该基因的重要性。为了尽可能多地保留数据中的生物信息，我们选择在本教程中避免对基因进行筛选。
归一化后，数据矩阵通常是对数 (+ 1) 转换的。这种转变有三个重要作用。
首先，对数转换的表达值之间的距离代表对数倍数变化，这是衡量表达变化的经典方式。
其次，对数转换减轻（但不消除）单细胞数据中的均值方差关系 (Brennecke et al,2013)。
最后，对数转换降低了数据的偏斜度，以适用于假设数据呈正态分布的下游分析工具。
虽然 scRNA-seq 数据实际上不是对数正态分布 (Vieth et al,2017)，但这三种效应使对数转换成为一种粗糙但有用的工具。差异表达检测 (Finak et al,2015；Ritchie et al,2015) 或批次校正 (Johnson et al,2006；Buttner et al,2019) 的下游应用强调了这种有用性，这些应用将对数转换用于这些目的。但是，应该注意的是，归一化数据的对数转换可在数据中引入虚假差异表达效应（预印：Lun，2018）。当归一化大小因子分布在试验组之间存在强烈差异时，该效应尤其明显。
问题和建议:
•我们建议使用scran对非全长数据集进行标准化。
另一种方法是通过scone评估基于平台的数据集的标准化方法。全长scRNA-seq协议可以使用bulk 方法修正基因长度。
•对于将基因的均值和单位方差缩放到0没有共识。
我们宁愿不缩放基因表达（We prefer not to scale gene expression.）。
•规范化的数据应该是log(x+1)-转换后用于假设数据是正态分布的下游分析方法。
Data correction and integration
如上所述的标准化试图消除计数采样的影响。但是，归一化数据仍可能包含不希望的变异性。数据进一步的校正针对技术和生物学变量，如批次、脱落或细胞周期效应。这些变量并不总是进行校正。相反，决定考虑哪些变量将取决于预期的下游分析。我们建议分别考虑生物和技术变量的校正，因为这些变量用于不同的目的，并且存在独特的挑战。
Regressing out biological effects
校正技术变量（covariates）对于揭示潜在生物信号至关重要，校正生物变量对于挑选出关注的特定生物信号更加重要。最常见的生物数据校正是去除细胞周期对转录组的影响。该数据校正可通过 Scanpy 和 Seurat 平台 (Butler et al,2018；Wolf et al,2018) 或具有更复杂混合模型（如 scLVM (Buettner et al,2015) 或 fscLVM (Buettner et al,2017)）的专门包装中实施的细胞周期评分的简单线性回归进行。用于计算细胞周期评分的标记基因列表来自文献 (Macosko et al,2015)。这些方法也可用于回归其他已知的生物学效应如线粒体基因表达，其被解释为细胞应激的指征。
在校正生物学效应数据之前，应考虑几个方面。
首先，校正生物学变量并不总是有助于解读 scRNA-seq 数据。虽然去除细胞周期效应可改善发育轨迹的推断 (Buettner et al,2015；Vento-Tormo et al,2018)，但细胞周期信号也可提供生物学信息。例如，可根据细胞周期评分确定增殖细胞群（参见 github 项目的个案研究）。
生物信号必须在语境中理解。鉴于生物过程发生在同一生物体内，这些过程之间存在依赖性。因此，纠正一个过程可能无意中掩盖另一个过程的信号。
最后，有人认为，细胞大小的变化解释了通常归因于细胞周期的转录组效应 (McDavid et al,2016)。因此，通过标准化校正细胞大小，或专用工具如 cgCorrect (Blasi et al,2017)，也部分校正了 scRNA-seq 数据中的细胞周期影响。
Regressing out technical effects
用于回归生物学变量的回归模型变量也可应用于技术变量。单细胞数据中最显著的技术变量是计数深度和批次。尽管标准化比例计数数据使细胞之间的基因计数相当，但计数深度效应通常保留在数据中。这种计数深度效应既可以是生物的，也可以是技术的。例如，细胞可能大小不同，因此 mRNA 分子计数也不同。然而，归一化后的技术计数效应可能仍然存在，因为没有缩放方法可以推断由于采样不佳而未检测到的基因的表达值。回归出计数深度效应可以提高轨迹推理算法的性能，它依赖于找到cell之间的转换（参见 project github 的案例研究）。当校正多个变量（例如，细胞周期和计数深度）时，应在一个步骤中对所有变量进行回归，以考虑变量之间的依赖性。
另一种基于回归的消除计数影响的策略是使用更严格的标准化过程，如下采样或非线性标准化方法（参见标准化部分）。这些方法可能特别适用于基于平板（plate-based ）的 scRNA-seq 数据集，其中每个细胞计数深度的较大变化可以掩盖细胞之间的异质性。
Batch effects and data integration
当细胞以不同的分组处理时，可能发生批次效应。批次效应可以由不同芯片上的细胞、不同测序泳道中的细胞或不同时间收获的细胞组成。细胞经历的不同环境会对转录组的测量或转录组本身产生影响。产生的影响存在于多个层面：实验中的细胞组之间、在同一实验室进行的实验之间或来自不同实验室的数据集之间。
在这里，我们区分前两个和后两个场景。
在相同的实验中校正样品或细胞之间的批次效应是经典的来自 bulk RNA-seq 的批次校正（Batch effects）。
我们将这与多次实验的数据整合区分开来，我们称之为数据整合（data integration）。虽然批效应通常使用线性方法校正，但一般使用非线性方法进行数据整合。
最近对经典批次校正方法的比较显示，ComBat (Johnson et al,2006) 在低至中等复杂度的单细胞实验中也表现良好 (Buttner et al,2019)。ComBat 由基因表达的线性模型组成，其中在数据的平均值和方差中均考虑了批次贡献（图 3）。不考虑计算方法，**批量校正的最佳方法是通过巧妙的实验设计预先消除影响并完全避免影响 **(Hicks et al,2017)。通过合并实验条件和样品中的细胞，可避免批次效应。使用诸如细胞标记 (preprint:Gehring et al,2018) 或通过遗传变异 (Kang et al,2018) 等策略，可能分离实验中合并的细胞。
与批次校正相比，面临的另一个挑战是整合不同的数据集。估计批效应时，ComBat 使用一批中的所有细胞来拟合批次参数。这种方法将混淆批处理效应与细胞类型或数据集之间不相同状态之间的生物学差异。为克服该问题，开发了典型相关分析 (CCA；Butler et al,2018)、相互最近邻 (MNN；Haghverdi et al,2018)、Scanorama（预印：Hie et al,2018）、RISC（预印：Liu et al,2018）、scGen（预印：Lotfollahi et al,2018）、LIGER（预印：Welch et al,2018）、BBKNN（预打印：Park et al,2018）和 Harmony（预打印：Korsunsky et al,2018）等数据整合方法。数据整合方法虽然也可以应用于简单的批次校正问题，但考虑到非线性数据集成方法的自由度增加，我们建议警惕过度修正。例如，在更简单的批次校正设置中，MNN 的表现优于 ComBat (Buttner et al,2019)。需要对数据整合和批次校正方法进行进一步比较研究，以评估这些方法的应用范围。
Expression recovery
另一种类型的技术数据校正是表达恢复（expression recovery）（也就是去噪或插补）。单细胞转录组的测量包含各种噪声 (Gru net al,2014；Kharchenko et al,2014；Hicks et al,2017)。这种噪音的一个特别突出的方面是dropout。推断dropout事件，用合适的表达值替换这些零，减少数据集中的噪声，一直是几个最新工具（MAGIC：van Dijk et al,2018；DCA：Eraslan et al,2019；scVI：Lopez et al,2018；SAVER：Huang et al,2018；scImpute：Li& Li, 2018).已证明执行表达恢复可改善基因相关性估计 (van Dijk et al,2018；Eraslan et al,2019)。此外，该步骤可与归一化、批次校正和 scVI 工具中实施的其他下游分析整合 (Lopez et al,2018)。虽然大多数数据校正方法以归一化数据作为输入，但一些表达式恢复方法是基于预期负二项分布噪声，因此在原始计数数据上运行。应用表达恢复时，应考虑到没有一种方法是完美的。因此，任何方法均可能导致数据中的噪声过度校正或校正不足。事实上，表达恢复的结果报告了假相关信号（Andrews & Hemberg, 2018).考虑到在实际应用中评估成功的表达恢复的难度，这个场景对考虑是否消噪的用户来说是一个挑战。此外，对于当前可用的表达式恢复方法，大型数据集的可伸缩性仍然是一个问题。鉴于这些考虑，目前对于如何使用消噪数据尚未达成共识（见处理数据章节的阶段）。谨慎的方法是仅将表达恢复用于数据的直观显示，而不是在探索性数据分析过程中应用之。这里彻底的实验验证尤为重要。
问题和建议:
•回归出(Regress out )生物变量只是为了轨迹推断，其他生物过程被没有回归出的生物协变量所掩盖。
•同时考虑回归技术和生物变量，而不是针对某项回归。
•基于平台的数据预处理可能需要回归计数，通过非线性归一化方法进行归一化或向下采样。
•当cell类型和批次之间的状态成分一致时，我们建议通过ComBat进行批次校正
•数据整合和批次校正应通过不同方法进行。数据集成工具可能过度纠正简单的批处理效应。
•用户应谨慎对待仅在表达恢复后发现的信号。最好不使用该步骤进行探索性分析。
Feature selection, dimensionality reduction and visualization
一个人类的单细胞 RNA-seq 数据集可以包含多达 25,000 个基因的表达值。这些基因中的许多基因对于给定的 scRNA-seq 数据集不会提供有价值的信息，许多基因将大部分包含零计数。即使在 QC 步骤中过滤掉这些零计数基因后，单细胞数据集的特征空间也可以有超过 15000 个维度。为了减轻下游分析工具的计算负担，减少数据中的噪声，并使数据可视化，可以使用几种方法降低数据集的维数。
Feature selection
scRNA-seq 数据集降维的第一步通常是特征选择。在此步骤中，将筛选数据集仅保留可提供数据变异性信息的基因。因此，通常使用高度可变基因 (HVG) (Brennecke et al,2013)。根据任务和数据集的复杂性，通常选择 1,000 至 5,000 个 HVG 用于下游分析（见图 EV1 和数据集 EV1）。Klein et al (2015) 的初步结果表明，下游分析对 HVG 数量的准确选择具有稳健性。在 200 到 2400 之间变化 HVGs 数量时，作者报告了 PCA 空间中类似的低维表示。基于这一结果，我们宁愿选择更高数量的HVGs。
Dimensionality reduction
特征选择后，单细胞表达矩阵的维数可以通过专门的降维算法进一步降低。这些算法将表达式矩阵嵌入到一个低维空间中，目的是在尽可能少的维度中捕获数据中的底层结构。这种方法的工作原理是单细胞 RNA-seq 数据本身就是低维的 (Heimberg et al,2016)。换句话说，细胞表达谱所在的生物流形可以用比基因数量少得多的维度来充分描述。降维旨在找出这些维度。
降维方法主要有两个目标：可视化和提取主要变化。可视化是尝试以二维或三维的方式对数据集进行最佳描述。这些缩小的尺寸用作散点图上的坐标，以获得数据的直观表示。相反，对于描述数据中存在的变异性，较高的组分变得不太重要。摘要技术可以通过找到数据的固有维数，将数据简化为其基本组成部分，因此有助于下游分析。虽然 2 维可视化输出不应用于汇总数据集，但可以使用汇总方法，使用领先的缩减组件对数据进行可视化，专门的可视化技术通常可以更好地表示变异性。
通过特征空间维度（基因表达载体）的线性或非线性组合减少维度。特别是在非线性的情况下，在这个过程中牺牲了降维的可解释性。一些常用的降维方法的应用示例如图 4 所示。随着越来越多的方法可供选择，详细回顾这些方法超出了本教程的范围。我们简要概述了可能帮助用户在常用降维方法之间进行选择的实际考虑。Moon et al (2018) 提供了单细胞分析降维的更详细综述。
两种流行的降维技术（主要是总结方法）是主成分分析 (PCA；Pearson，1901) 和扩散图（diffusion maps ）（Coifman 等，2005），Haghverdi 等 (2015) 推广用于单细胞分析。主成分分析是一种线性方法，通过最大化每个进一步维度中捕获的残差来生成降维。尽管 PCA 并不像非线性方法那样能够捕获很少维度的数据结构，但它是目前许多可用的聚类或轨迹推断分析工具的基础。PCA 是一种常用的非线性降维预处理方法。通常，PCA 通过其前 N 个主成分汇总数据集，其中 N 可以通过肘部启发式（见图 4F）或基于置换测试的 jackstraw 方法（Chung Storey, 2015; Macosko et al, 2015).PCA 的简单线性的优点是在减少的维空间距离在这个空间的所有区域有一致的解释。因此，我们可以将感兴趣的数量与主成分相关联来评估它们的重要性。例如，主成分可以投影到技术干扰协变量上，以研究 QC 的性能、数据校正和标准化步骤 (Buttner et al,2019)，或显示基因在数据集中的重要性 (Chung et al,2019)。由于扩散成分强调的是数据中的转换，它们主要用于连续过程（如差异）感兴趣的情况。通常，每个扩散组分（即扩散图维度）突出显示不同细胞群的异质性。
Visualization
出于可视化目的，使用非线性降维方法是标准实践（图 4）。scRNA-seq 可视化最常用的降维方法是 t 分布随机邻域嵌入（t-SNE；van derMaaten & Hinton, 2008)。t-SNE 降维以全局结构为代价来获取局部相似性。因此，这些可视化可能夸大细胞群体之间的差异，并忽略这些群体之间的潜在联系。另一个困难是选择其复杂度参数，因为 t-SNE 图可能显示不同数值的簇 (Wattenberg et al,2016)。t-SNE 常用的替代方法是UMAP（预本：McInnes & Healy, 2018) 或者基于图的工具，如 SPRING (Weinreb et al,2018)。UMAP 和弹簧力导向布局算法 ForceAtlas2 可以说是底层拓扑的最佳近似值（Wolf et al,2019,Supplemental Note 4)。该比较中 UMAP 的不同之处在于其扩展至大量细胞的速度和能力 (Becht et al,2018)。因此，在没有特殊生物学问题的情况下，我们将 UMAP 视为探索性数据可视化的最佳实践。而且，UMAP 还可以在两个以上维度汇总数据。虽然我们不知道 UMAP 在数据汇总中的任何应用，但它可能证明是 PCA 的合适替代方法。
细胞水平上经典可视化的替代方法是基于分区的图形抽象 (PAGA；Wolf et al，2019)。该工具已被证明可以充分近似数据的拓扑结构，同时使用集群粗粒化可视化。结合上述任何一种可视化方法，PAGA 都会产生粗粒度的可视化，这可以简化单细胞数据的解释，尤其是用大量细胞的时候。
问题和建议:
•我们建议根据数据集的复杂性选择1000到5000个高度可变的基因。
•当基因表达值被归一化为零均值和单位方差时，或当模型拟合的残差被归一化表达值时，不能使用使用基因表达均值和方差的特征选择方法。因此，在选择HVGs之前，必须考虑要进行什么预处理。
•应分别考虑降维方法进行总结和可视化。
•我们推荐使用 UMAP 进行探索性可视化；使用 PCA 进行一般性总结；使用扩散图替代 PCA 进行轨迹推断总结。
•具有 UMAP 的 PAGA 是可视化特别复杂数据集的合适替代方案。
Stages of pre-processed data
虽然我们已经将 scRNAseq 中常见的预处理步骤概述为上述工作流程，但下游分析通常倾向于采用不同水平的预处理数据，建议根据下游应用调整预处理。为了向新用户阐明这种情况，我们将预处理划分为 5 个数据处理阶段：
原始数据，
标准化数据，
校正数据，
特征选择数据，
降维数据。
这些数据处理阶段分为三个预处理层：
测量数据、
校正数据和
缩减数据（降维）。
应始终进行细胞和基因 QC，因此从此处省略。而这些处理层的顺序代表了 scRNA-seq 分析中的典型工作流程，也可以跳过某个处理层或者在处理阶段的顺序上有轻微的改变。例如，对于单批数据集，可能不需要进行数据校正。在表 1 中，我们总结了每一层预处理数据的适当下游应用。
表 1 预处理阶段分为实测数据、校正数据和缩减数据 3 组。我们将测量数据定义为原始数据和保留零结构的处理数据。通过使用cell特定因子缩放计数数据，全局缩放规范化方法即使在 log (+ 1)转换之后也保留 0 表达值。相反，纠正不需要的变异性数据替代零表达值。校正后的数据层代表数据的最干净版本，是基础生物信号的最接近近似值。我们称最后的预处理层为缩减数据。该数据层强调数据的主要方面，可以使用简化的功能集进行描述。
上述特征决定了预处理数据对于特定下游应用的适用性。作为最后的预处理阶段，缩减数据将是广泛适用的数据层的候选。然而，差异表达检测仅在基因空间中进行生物学解释，未（完全）用简化数据表示。约简数据的作用在于生物学的总结和噪声的减少，可能掩盖生物信号。因此，缩减数据用于需要数据总结（可视化、邻域图推理、聚类）的探索性方法和计算复杂的下游分析工具（轨迹推理）。的确，许多轨迹推理方法在工具本身中加入了降维过程。
单个基因的表达谱只能在基因空间中进行比较，在测量和校正数据中捕获。表达谱的比较可以通过可视化和统计学进行。我们认为应该对不同的数据层进行可视化和统计比较。基因表达的目测检查，校正数据最为合适。如果提供原始数据进行可视化比较，则要求用户固有地理解数据中的偏倚，以解释结果，校正数据有助于这种解释。然而，此处应单独考虑技术和生物学变量的校正数据。对生物变量的校正可能增加特定生物信号的强度，也将掩盖可能相关的其他信号。因此，生物校正数据主要适用于关注特定生物过程（如轨迹推理方法）的分析工具。
基因表达的统计学比较在测量数据层上最合适。没有完美的数据校正方法可用于消噪、批次校正或其他变异来源的校正。因此，数据校正方法不可避免地对数据进行了过高或过低的校正，因此以一种非预期的方式改变了至少一些基因表达谱的方差。基因表达的统计检验依赖于评估背景方差，作为数据中噪声的无效模型。由于数据校正倾向于减少背景变异（图 EV2），背景变异被数据校正方法过度校正的基因将更可能被评估为显著差异表达。此外，某些数据校正方法（例如 ComBat）将不符合实验设计的表达信号解释为噪声，随后从数据中删除。除了低估噪声外，这种实验设计信号的优化会导致高估效应量。鉴于这些考虑，使用测量数据作为输入，而不是使用校正数据，构成了对差异试验更保守的方法。使用测量数据，在差异检验模型中可以并且应该考虑技术变量。
上述观点得到了最近一次scRNA-seq差异分析方法比较的支持，该方法仅使用原始数据和规范化数据作为输入(Soneson & Robinson, 2018)。本研究使用的归一化数据仅围绕全局标度方法。然而，目前许多可用的非线性归一化方法模糊了归一化和数据校正之间的界限(参见“归一化”一节)。这种标准化的数据可能不再适合作为差异分析的输入。
问题和建议:
•使用测量数据（measured data）进行统计检验，根据发现的生物数据流形，对数据进行可视化比较时使用修正数据（corrected data ），对其他下游分析时使用简化数据（ reduced data ）。
往期精彩
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
单细胞天地欢迎你