scRNA-seq计算方法的优势和局限性
文章信息
文献标题：The triumphs and limitations of computational methods for scRNA-seq
发表时间：2020.06.21
发表杂志：Nature Methods（IF=30.822）
原文链接：https://doi.org/10.1038/s41592-021-01171-x
摘要
在过去的10年里，伴随着单细胞转录组测序技术的飞速发展，单细胞数据分析的计算方法也获得了相应的长足进步。随着实验技术的生产力和准确性的提升，新兴算法的开发也日益揭示了更复杂的生物学层面的信息，例如细胞类型的组成以及发育动态的基因调控等。同时，这种飞速发展也迫使我们不断地重新评估底层统计模型、实验目的以及数据处理的庞大体量。本文作者回顾了单细胞RNA测序（Single-cell RNA sequencing，scRNA-seq）分析的基本计算步骤，审视了不同方法背后的假设支撑，并强调了其中的成功典范、尚存的局限性和歧义之处。
引言
作者首先简要总结了scRNA-seq分析的基本步骤（图1和2）。scRNA-seq分析方法多种多样，但大多包含以下结构：
所测量信息的统计学模型； 数据在降维空间中的反映； 表达谱流形特征的近似反映（其中最简单且常见的近似方法是形成离散的转录亚群）。
细胞的统计学观点
与bulk RNA-seq相比，单细胞实验中每个细胞捕获的起始转录本太少，并且在不同的细胞之间，捕获的转录本量可相差不止一个数量级，因此，对细胞转录状态的测量存在更大的不确定性。整体上看，只要数据能够捕获足够的信息以识别相似的细胞，生信算法就能通过”借用“不同细胞的信息来推断这一群相似细胞整体的转录状态。换句话说，只要测序的细胞足够多，就能弥补单个细胞转录本捕获的不确定性（图3b）。当然，这同时需要付出更大的成本和劳动力。
基于单细胞数据的稀疏性和不确定性，从概率的观点探究基因表达状态可能更有现实意义，即在给定细胞中，考虑某一转录本有某一特定表达水平的似然性（图3c）。目前大多数scRNA-seq的统计模型是由负二项模型扩展而来，并包含更复杂的均值-方差关系，或是引入混合组分来处理scRNA-seq中大量出现的零值。其中，作者提到尽管scRNA-seq数据存在大量零值，但由于UMI的广泛应用，零值的实际比例并不高于mRNA抽样时期望的比例，零膨胀又一次被diss了……
现有的统计模型几乎都假设scRNA-seq是对单个细胞中mRNA分子的随机抽样，然而，尤其是当数据来自不同单细胞protocol时，真实的抽样率却呈现基因特异性偏倚。当涉及到跨样本差异表达分析时，如果样本来自同一建库测序方法，一般来说，这种基因特异性的抽样偏倚就影响不大；反之，跨平台样本的差异表达分析就会出现问题。
转录状态的比较
检验基因转录水平在任意两（群）细胞之间是否存在统计学差异的方法，是由scRNA-seq数据的统计模型决定的。例如，对给定的基因，如果能在每个细胞中估计其表达水平的后验概率分布，那么简单地计算其重叠部分（即表达水平在两细胞间相同这一事件的概率）即可完成差异表达检验（图3c）。目前针对scRNA-seq的差异表达检验方法包括参数或非参数模型。当比较的细胞群体较大（例如>100细胞）时，与参数模型相比，标准非参数检验（如Wilcoxon检验）在满足更少的前提假设下，其统计效能更好。相应地，统计效能的增加必然导致检验所得的差异表达基因数量增加（图3d）。这时，考量细胞亚群的定义以及差异表达的数量级（例如logFC）就显得更为重要了。
个人理解： 提高统计效能必然增加犯”一类错误“（假阳性）的风险，联系我们平时使用Seurat包的经验，FindMarkers()的默认检验方法就是Wilcoxon检验，得到差异分析结果后，往往还需要根据logFC和p_adj进一步筛选基因。
除了直接比较平均表达值水平，scRNA-seq还能比较诸如基因表达的可变性（Variability）或分布形状。作者提到，Martinez-Jimenez等人的研究指出在不同年龄CD4 T细胞中某些基因的平均表达水平相似，但其表达可变性却随着年龄的增长而增加。
最后，差异检验对批次效应（Batch effects）尤其敏感。大多数现有方法能控制简单的批次结构，但实验设计一旦涉及诸如多样本分类或匹配，就需要更精心设计的协变量控制。目前主流还是选择bulk RNA-seq的检验方法，而忽略细胞类型内部的细胞间差异。
除了差异表达分析，转录状态的比较还包括定量细胞在表达谱空间的相互距离。这一点贯穿在细胞分群、谱系推断或二维可视化等下游分析中。距离定量方法主要由两类：直接衡量基因转录差异的传统距离（如Euclidean、L1和Canberra），以及从统计偏差的角度衡量两细胞是否相等（如基于Poisson的相似性和Jensen–Shannon散度）。第一类方法可能无法区分不同维度上的变化，而第二种方法十分依赖测序深度和细胞覆盖度。此外，传统距离方法并不适应scRNA-seq固有的高维特点，在原始的基因表达空间中难以区分最近或最远点。幸好，通过在更低维空间中描述转录差异能有效应对”高维诅咒“，使细胞间距离变得可计算。
降维的需求
scRNA-seq数据的有效维度比所有可能的基因表达组合模式所允许的维度更低，部分原因在于scRNA-seq数据的视野有限：例如，对于描绘某一只有三个主要细胞类型的组织来说，大多数转录改变只需要用2个捕获细胞类型间差异的轴即可解释。根本上，更低的有效维度反映了细胞中基因调控逻辑所带来的转录协同性。例如，当10个基因受同一转录因子激活时，其表达特征用单一变量（而非10个变量）即可描绘。
scRNA-seq中最常用的降维方法是主成分分析（Principal component analysis，PCA），即寻找基因表达特征的有限个线性组合，使之尽可能多地解释数据中的方差。细胞间主要的转录方差则位于解释度最高的一部分PCs所描绘的低维超平面中。
scRNA-seq的count数据本身的特性使得PCA的应用更加复杂。首先，单个基因的方差很大程度上取决于其表达量（图3e），而top PCs会以牺牲整体转录模式为代价，重点关注高表达基因的细微波动。因此，我们需要比较观察到的转录本方差和根据表达量计算的方差期望值。最初，人们通过外源性引入spike-in对照来估计方差期望。由于会增加实验设计的复杂性，这种方法逐渐失宠了。作为外源性对照的替代品，大多数分析方法基于所有观察到的转录本来估计方差期望（图3g）。某给定基因的方差如果高于期望值，其表达模式就更可能区分主要细胞亚群（图3f-i），于是我们可以将PCA限定在这一小群高可变基因（Highly variable genes，HVGs） 中。这种统计模型也可以用于对每个基因的方差进行归一化，以计算超出期望的残差方差（SCTransform就是干这个的）。由于方差估计是基于所有被测细胞，高可变基因会更关注最剧烈的亚群差异（例如上皮细胞vs免疫细胞），而更微妙的差异（例如CD4和CD8 T细胞）可能需要对分群进行重分析。
影响PCA应用的另一点是scRNA-seq数据的稀疏性。PCA实质上是对协方差矩阵的优化分解，因此最适用于对称分布的数据。诸如log等方差稳定性变换（Variance-stabilizing transformation，VST） 常用于将count值的分布归一化为接近正态分布。然而，大量的零值使得各种变换结果都会出现一个偏态峰（图3j）。此时，top PCs会主要捕获细胞间的测序覆盖度（技术误差）。减轻这一误差的方法包括减少零值的权重或对每个细胞的总UMI count值进行回归（Regress out）。对后者来说，如果不同细胞类型本身的RNA量就有差异（生物学效应），直接对UMI count值回归反而引入新的混淆。此时，可以选择迭代分析的思路：先识别主要细胞类群，在亚群重分析时对测序深度回归。
PCA的局限性反映了一个基本事实：低维优化本身和scRNA-seq的底层统计学结构是脱节的。因此，诸如PCA等传统降维方法可能会由于捕获技术噪音、牺牲真实的生物学差异而遭到贬低。使用更一般的**因子分析（Factor analysis，FA）**方法能将因子分解与适当的统计模型相结合。例如零膨胀的因子分析（Zero-inflated factor analysis，ZIFA），将scRNA-seq数据理解为来自正态分布模型和零值插入模型的混合，以此接纳数据中大量的零值。除了正态模型，count模型（例如负二项分布）可能更准确地描述了scRNA-seq数据的统计学特征。尽管将表达矩阵分解为多个计数过程组分在统计学上仍有难度，一些方法通过相关的分层途径，利用单一计数过程对矩阵建模，计数过程的参数则由各个组分的线性组合驱动。例如ZINB-WaVE利用广义线性模型来驱动平均表达量和零膨胀负二项分布的drop-out参数。
转录本计数和降维空间的映射不一定局限于线性关系。类似地，描述scRNA-seq数据特性的统计参数（例如过离散量），也可能依赖于不同基因转录本的生物学特征（如GC含量、亚细胞定位）或细胞类型特征（如裂解易感性、酶成分）。然而，全面考虑这种拓展的潜在联系十分困难。为了映射更多复杂的非线性关系，一些科学家尝试使用自编码的神经网路（Autoencoder neural networks），以便学习复杂的非线性多维函数。例如，t-SNE本身是基于对高维空间中邻近细胞距离的经验优化，其坐标与基因表达状态之间缺少明确的解析函数（Analytical function），而神经网络编码的一系列变换则有足够的效能提供解析函数的某种近似（图4）。
简单来说，自编码器是一种神经网络，用于寻找某种函数，使之能以获得最优重建的形式，在原始数据和低维空间反映之间建立相互映射。与线性方法相比，采用非线性函数所得的降维空间能更有效地捕获细胞亚群的结构。然而，非线性映射对隐状态（Latent state）的解释可能有困难，使用非对称自编码器（Asymmetric autoencoder）将非线性编码和解释能力更高的线性解码相结合，是潜在的缓解方式之一。
scVI算法通过结合神经网络和概率模型，将表达谱检测的不确定性播散到隐空间（Latent space）。利用scRNA-seq的零膨胀负二项模型，scVI构建了一种变分自编码器（Variational autoencoder），在高维的转录空间和降维的隐空间之间拟合非线性映射。伴随网络（Companion networks）能追踪低维空间接收或发出映射的不确定性，并支持诸如聚类、轨迹推断和差异表达等下游分析，以便考虑scRNA-seq数据的不确定性。
利用邻居图（Neighbor graphs）实现细胞基因表达流形的近似
戏剧演员H.G. Wells在《隐形人》中抱怨，如果雪落在他身上，掉落的雪花会暴露他的形体。同样地，通过观察所测细胞的分布可以近似表示细胞表达特征的流形。常用方法之一是将邻近细胞在一定维度空间中相连，构成最近邻居图（Nearest-neighbor graph）。当细胞的抽样密度足够时，遍历此图即可跟踪流形的形状（图5a）。
邻居图由于其计算和解析优势，已经被用于大量scRNA-seq分析中。典型的图构建方法是通过高度优化的”近似最近邻“（Approximate nearest neighbor，ANN）搜索算法，或诸如共同邻居图（Shared-neighbor graph，SNN）等补充优化方法，将细胞和其k个最近邻（k nearest neighbors，kNN graph）相连。其他更精细的图构建方法则被开发用于整合多个数据集，以便克服多种批次效应，例如通过相互最近邻（Mutual nearest neighbors，MNN）算法对来自不同数据集的节点（细胞）建立额外的边（edges）进行匹配。解析邻居图结构的技术也已被用于scRNA-seq，例如通过对邻居图进行谱分析（Spectral analysis）来鉴定表达谱流形的主要特征。特别地，当图的Laplacian矩阵拥有最小特征值时，其特征向量与连通分量（Connected components）及图的主要结构轴有关。拉普拉斯特征映射可以有效地导出表达谱流形的精细低维反映，促进细胞异质性分析和可视化。
邻居图强化了包括亚群分析、动态轨迹推断或数据集映射在内的许多下游分析。当数据中含有大量亚群时，不同亚群的基因表达差异可能位于表达谱空间的不同平面，且具有剧烈变化的量级，此时如PCA等线性降维方法需要不止两个维度来反映亚群在数据空间的相对位置。相反，t-SNE和UMAP能在全局尺度上缩减距离的量级和方向，保留局部近邻关系。这些关系能在二维平面有效地捕捉细胞分群、连续轨迹或其他数据结构。
细胞聚类
在众多根据转录相似性聚类细胞的算法中，最常用的是通过邻居图鉴定倾向于相互连接的一个个细胞”社区“（community），如Louvain和Leiden算法。
细胞聚类使我们能有效地探索和解释数据，指导诸如差异表达等下游分析。然而，细胞聚类本身作为一种近似方式，注定无法携带除了转录相似性以外的生物学意义。细胞聚类可能区分不同细胞类型，捕捉同一细胞类型中微妙的转录改变，或是基于随机误差而过度划分某群完全一致的细胞。多数算法都允许用户经验性地调整分辨率（图5b,c）。事实上，聚类分辨率取决于对数据阐释的需要，例如某些案例中将T细胞统归入一个分群就足够了，而另一些案例可能需要将其进一步细分。包括Walktrap随机游走在内的层次聚类算法，能关联不同分辨率下的聚类结果，有助于描述大类群和细分亚群之间的关系。然而，当细胞周期或功能激活等来源的转录差异存在时，简单的分层可能就歇菜了（图5d-g）。
动态过程分析
在多数生物学背景如机体发育或刺激响应中，转录动态都是必然关注的问题。目前的scRNA-seq在设计之初是为了捕捉某一特定时间点的细胞快照，然而，底层的某些转录动态特征是可以推断的。譬如蚂蚁的照片，通过跟踪蚁群密度的细丝，可以描绘蚂蚁行进的路径（图6a）。最优路径的搜索通常有赖于假设特定拓扑结构，例如树形（主图问题，Principal graph problem，图5c,d）或曲线形（主曲线问题，Principal curve problem）。这种在低维空间追踪细胞密度的一般途径最初应用于monocle包，目前正被多种方法进一步细化，并逐步成为最常用的推断转录动态的途径，成功地用于捕捉细胞分化轨迹的分支或其他生物学背景下的动态过程。应用轨迹推断时应当注意，这些方法通常倾向于将所有输入的细胞用轨迹连接起来，而不考虑其是否真的参与该动态过程。因此，事先确定细胞亚群，严格限定用于轨迹推断的数据范围十分重要。同时，在解释结果时，时刻牢记轨迹推断的某些细节可能有不确定性。这种不确定性可能起源于表现某一特征的细胞亚群数量有限（例如某些分支只包含少数几个细胞）。大多数方法通过诸如PC轴、扩散成分轴甚至二维嵌入（UMAP等），在降维空间中拟合细胞的轨迹（monocle3等）；另一些方法通过连结细胞聚类而非单个细胞来简化这一问题（Slingshot、TSCAN或PAGA）。因此，应时刻牢记聚类或降维算法也可能对轨迹推断结果产生显著影响。
本质上讲，”蚂蚁踪迹“法依赖于观察在轨迹中每个点的细胞。换句话说，它假设潜在的动态过程是静止的、遍历的，因此样本量足够大的scRNA-seq数据能覆盖某群细胞所有可能的中间状态。然而，在伴随发育时间改变或外部刺激的情形下，单一数据集的细胞会呈现离散的聚类而非连续的轨迹。通过混合测序不同样本或多个时间点的细胞可能重新引入异步性，然而，除非有足够的特异性方法来设计多时间点实验，否则轨迹推断不能将取样时间考虑在内（啊……感觉某大佬不是很喜欢拿胚胎时间点拟合分化轨迹的吗？）。
基于密度的方法致力于对细胞经过的核心轨迹取近似，而对轨迹的细节留有大量不确定性，因为在相同密度模式下可能出现完全不同的情景（图6b）。最明显的歧义是细胞流的方向，多数分析方法依赖于外部知识来确定方向，例如使用鉴定过的marker基因确定轨迹的开头和终点。目前，Grun等人基于”干祖细胞比分化后的细胞倾向于表达更大范围的转录本“这一趋势，提出了一项细胞分化轨迹的系统评价标准。通过定量每个细胞基因转录的Shannon信息熵，该方法能自动分配轨迹方向使得细胞信息熵减少。这种熵启发式方法已用于多种研究发育的案例，但尚未推广到其他动态过程类型中。
除了细胞密度的”蚂蚁踪迹“，推断scRNA-seq快照中的mRNA生命周期动力学也能提供动态过程的额外信息。审查mRNA分子在其生命周期不同阶段（即新生转录、加工和降解）的平衡，可以推断出一个细胞在拍摄快照前后的动态。例如，如果一个细胞产生的新生mRNA量比补偿原有mRNA降解所需的量更多，我们可以推断，下一个时间点的转录本丰度将增加。常见的scRNA-seq方法能定量pre-mRNA（未剪接）和成熟mRNA（剪接过）之间的相对丰度，即RNA速率（RNA velocity），然后用于近似表示细胞在测序前后的一小段时间内可能的状态。RNA速率能揭示复杂的轨迹流形式，如分支、循环或逆流（图6c,f）。该途径不需要遍历性假设，然而，它仍有局限性。实际的动态过程必须与mRNA剪接动力学在时间尺度上具有可比性。目前的估计方法局限于具有足够长内含子区的转录本，因此可能引入明显的不确定性，并且需要在细胞的邻近关系中取平均来提高稳健性。不管怎样，RNA速率显示，通过对细胞的静态分子快照进行更详细的考量，可以推断出基因转录的时间动态。更细化的参数拟合技术能提高这种推断的准确性，同时，通过新的实验技术捕捉mRNA、表观遗传或蛋白质的生命周期，也可以提高RNA速率的推断能力。
展望未来：整合多个样本、分子模式和物理空间
今天，大多数在单细胞基因组学方面的分析致力于处理稀疏、嘈杂但富含信息的数据。随着技术的稳健性和可及性逐年增加，scRNA-seq已经成为生物学探索的重要工具。解决真实世界的问题需要更细化的实验设计。例如，疾病研究往往需要对不同分组的个体进行比较，而其他研究可能涉及纵向样本收集或多组织多样本的比较。scRNA-seq的进展催生了包括DNA甲基化、染色质可及性或蛋白质丰度等分子模态的单细胞组学。为了提供每个细胞在空间背景上的信息，空间转录组应运而生。
所有的进步都要求科学家开发能够考虑额外层面信息的分析工具，然而，从scRNA-seq数据中学到的经验教训在多组学层面仍然适用。例如，大量方法被用于“对齐”（Align）单细胞数据集，以便鉴定跨数据集相关的细胞类群。它们大多依赖于相同的归一化和降维步骤，并拓展了近邻关系图使之能映射来自不同数据的细胞。这些方法反映了早期批次校正问题的推广，并且可以用于整合不同模态的数据，如转录组和染色质可及性。尽管有新的计算方法开发出来追赶单细胞数据不断增加的体量和复杂度，但这些早期方法目前仍保持着相当的重要性。
写在后面
好家伙大佬的文献光翻译就用了三天……说实话，这篇综述的行文逻辑让人非常舒服，可以说一步一步地带我们系统回顾了scRNA-seq数据分析的基本思路。
作者略过了数据质控的问题，上来就直指scRNA-seq分析的根本——差异分析，即从数据中细胞（亚群）、样本等不同维度的差异； 说到差异分析，自然要提统计模型，于是作者总结了scRNA-seq数据与生俱来的统计学特征及目前的主流模型； 但是，数据的稀疏性和不确定性使得一般统计方法还远远不够，于是作者告诉我们，可以通过数据空间中衡量细胞间距离来评估相似性； 但是，传统计算距离的方法在高达几万个维度的基因表达空间中陷入”高维陷阱“，并且许多基因的表达规律是有相关性的，于是作者引出了降维的思路，回答”为什么要降维“以及”为什么可以降维“； 单细胞数据本身的稀疏性和噪音使得直接在count水平降维会引入新的混淆，于是作者说，要归一化，要特征选择； 仅用一般的线性降维可能无法完全将原数据特征反映在低维空间中，于是作者引出了目前的替代方案（因子分析、神经网络等）； 细胞们的基因表达特征可以视为一种流形，而降维后，通过评估细胞的远近能将这种流形近似表示出来（邻居图）； 回到最初的问题，数据分析的目的是寻找相似和差异，聚类分析可以将细胞在数据空间中的近邻（相似）关系直观地反映，并指导后面的差异表达分析和轨迹推断。
到今天，单细胞的数据分析大多已经流程化，然而，在一遍又一遍机械地重复那些既有流程时，我们作为生信工作者，应当时时重温数据分析中的底层逻辑，而不是最后沦为调包侠和调参侠。
生信人，不忘初心！