【综述】Nature Methods | 干货！一文读懂单细胞转录组分析的现状和问题！
Kharchenko, Peter V. The triumphs, and limitations of computational methods for scRNA-seq. Nature Methods (2021): 1-10. PMID: 34155396
导语
各位亲爱的小伙伴们，本期又是小师弟来和你们唠一唠生信的文章了，有没有谁想念我了呢
？哈哈，不知道大家在科研中是不是经常接触单细胞转录组数据呢？相信大伙儿都能找到各式各样的教程文章带你一步一步地做单细胞分析，然后跑出各式各样漂亮的图。但是你知道这背后每一步可能碰到的陷阱和障碍吗
？做科研要求我们不仅知其然还要知其所以然，知己知彼方能百战百胜。近期这篇发表于nature methods上的文章就以review的形式探讨了单细胞转录组数据分析中每一步所蕴含的假设以及可能会遇到的陷阱
。小师弟觉得十分有用，因此特地写了这篇推送文章来和大家一起解读一下这篇文章的要点（强烈建议小伙伴们访问原文地址进行详细品读）。为了方便大家快速了解本文内容以及选读感兴趣的部分，我在摘要之后以图文形式列出了本文的要点。
1
摘要
单细胞转录组测序技术在过去数十年间得到了飞速发展，与此相伴而生的是各类计算生物学方法的大量涌现。随着检测技术的精度和通量的不断提升，生物学各方面的复杂性在新算法的帮助下被逐步揭示，这些复杂性具体体现在细胞群体组成、基因调控以及细胞动态演变等各方面上。与此同时，计算生物学的快速发展使得人们需重新评估各种算法所依赖的统计模型，所要实现的实验目的，以及所要处理的巨大数据量。在此篇文章中，作者回顾了单细胞转录组数据分析的关键步骤，总结了不同方法所依赖的假设，并指出这些方法的成功之处和仍然存在的不足之处。由于单细胞转录组测序技术已成为生物学研究的主流技术之一，因此本文的讨论对于广大科研人员具有极大的参考意义。
2
要点总结
3
要点解读
1. 如何以统计的视角看待单细胞数据
单细胞转录组数据由于其技术的限制，其所存在的不确定性要大于常规RNA-seq数据，其中的原因主要包括两点：1）每个细胞的总RNA量仅有一小部分能被仪器所检测到；2）单细胞之间被捕获的RNA量存在巨大差异。因而测序技术通常通过提高单细胞的捕获数量来抵消随机性。
另一方面，研究者们也尝试通过概率模型来解释观测到的随机性。很多研究采用了负二项分布模型来对基因表达的分布进行建模。需要注意的是，许多研究认为随机性导致单细胞数据包含了过多的零值，因而这些研究的模型通常包含了一个零膨胀(zero-inflation)的部分来解释过多的零值。但一项既往研究的实验结果显示：使用UMI(uniquemolecular identifier)技术的单细胞数据已能够极大地缓解这个现象[1]，已没有必要使用包含零膨胀的模型。
2. 比较转录状态
在建立了模型之后，下一步就是利用模型对基因的转录状态进行比较，这就是常见的差异表达分析，其目的是比较一个基因在两个细胞群体间表达量的差异。通过统计模型可知基因在两个细胞群体中的分布形状，那么两个分布之间的重叠部分就是基因表达在两个群体中相同的概率（见下图），可理解为统计检验中的P值（红字为个人观点，有待商榷）。在一些研究中，人们使用参数模型（比如负二项分布，高斯分布等）来完成分布的估计和差异表达分析，然而当单细胞数量增多后，参数模型不仅计算缓慢并且效果不佳。因而在许多研究中，研究者倾向于采用计算速度较快的非参数分布来完成差异表达分析(例如 Wilcoxon秩和检验)。不过需要注意的是，如果细胞群体足够大，几乎所有的基因都会呈现统计学差异(见下图)，因此引入其他标准来辅助差异基因的筛选就变得十分重要(例如设定一个fold-change阈值，使用更小的细胞亚群进行差异表达分析等)。
除了进行差异表达分析外，比较转录状态还包括了计算细胞与细胞间的距离或相似性（此步骤是下游其他分析诸如细胞聚类，可视化以及轨迹分析的基础）。若将单细胞的每个基因视作一个维度，那么毫无疑问单细胞数据是一种高维数据，在高维空间中计算细胞间的距离将不可避免的面临“维数灾难”(curse of dimensionality)的问题。简而言之，在高维空间中，常用的距离度量方式(欧式距离)将无法区分距离近的细胞和距离远的细胞。所幸的是，可通过降维的方式来缓解这一问题，下一小节将会集中探讨单细胞数据降维的问题。
3. 高维度数据降维
单细胞数据的有效维度数量要比其原始维度低得多。何为有效维度？有效维度是表征单细胞数据变化程度所需的最低维度。从生物学意义上也可以很好地理解这个问题。基因组中的许多基因通常存在功能的高度相关性，因此原始的高维数据冗余程度极高，而降维则等同于寻找最简约的表征数据的方式。
第二大特点是单细胞数据的稀疏性。由于PCA的原理为求协方差矩阵的矩阵分解，因此只有当数据分布呈对称形状时PCA才能达到最佳效果。然而，如前所述，单细胞数据含有很多零值，导致数据分布形状向零值方向偏斜（见下图），整个分布呈极为不对称的形状。针对此，既往文献中采用了一些校正方式，可供参考[2-4]。
PCA之所以无法克服以上两大困难，原因在于其计算矩阵分解的过程中没有考虑到单细胞数据自身的统计分布。针对此，既往研究采用的思路为把数据分布建模和矩阵分解结合到一起，例如ZINB-WaVE就是其中的代表之一，详情可见[5]。
PCA归根结底仍只是一种高维到低维空间的线性投影(linear projection)。例如照片就是一种高维到低维的线性投影方式，把三维空间的物体投影到维度更低的二维空间（在此小师弟又要安利一下三体小说，关于维度的描述实在太棒
）。然而投影方式不一定要局限于线性投影，非线性投影很多时候能够更好地反应数据本身的结构，现今有关单细胞分析的文章中常用的非线性投影方式有UMAP、t-SNE以及autoender等。(小师弟个人观点，UMAP速度快，同时能够保留远距离细胞群体间的信息，对下游分析诸如轨迹分析和细胞聚类等很有帮助)。
4. 利用最近邻图构建流形
在先前的小节里已经提到了高维数据存在极大冗余性，因此需要采用降维的手段来分析数据。除此之外，我们还可以通过流形在低维空间中表征高维空间的数据。何为流形？流形简而言之就是一种低维空间的几何表面，其所代表的区域能够逼近单细胞数据的分布。现今学术界普遍采用的方式是用最近邻图(nearest-neighbor graph)来估计流形的形状。最近邻图的构建简单明了，只需通过一定手段计算细胞与细胞间的距离，再将距离作为边，细胞作为节点即可构建最近邻图（也称为最近邻“网络”，见下图）。最近邻图在单细胞数据分析中的应用十分广泛，前一小节提到的UMAP和t-SNE降维方法实际上就是最近邻图在二维空间上的可视化，与PCA相比，它们不受表达值量级和分布形状的影响，可以很好地表示高维空间中细胞与细胞间的局部距离和数据的总体结构。
5. 细胞聚类
细胞聚类也是单细胞数据分析中一个重要步骤。聚类的结果是将细胞分为了不同的群体，以方便下游的进一步分析，比如在群体之间进行差异表达分析。通常来说，其依据的基础是利用最近邻图来鉴定不同的细胞群体，使得群体中细胞与细胞的联系比群体间细胞与细胞的联系更为紧密。常用的聚类算法有Louvain聚类和Leiden聚类等。不过，这些算法所依据的只是细胞转录组数据的相似性，实则并没有考虑到真正的生物学意义。因此，对聚类结果的解释一定要结合实际生物学问题。这使得在实践中，对聚类的分辨率大小的选择变得十分重要。例如在有些研究中需要将T淋巴细胞视作一个大类群而另一些研究中则需要将T细胞细分为不同亚群，这时候就需要给予聚类算法不同的分辨率等级。一些层次聚类算法，例如walkTrap[6]就能很好捕捉到不同分辨率等级下的细胞类群，为分析提供了灵活性。
6. 细胞的动态演变过程
单细胞转录组数据是探究器官组织发育以及细胞类型分化演变的利器之一。现有的转录组测序技术可对细胞群体在某个时间点的状态拍一个“快照”。在对原始数据降维之后，在低维空间中可通过这个静态的“快照”推断细胞群体的演变轨迹，这就是常见的单细胞轨迹分析。Monocle是目前基于这一原理最常被使用的轨迹分析工具(最新版本为Monocle3)[7]。需要注意的是，并非所有的细胞都参与了演变或分化的过程，但是此类轨迹分析工具通常将所有细胞都作为输入用于构建轨迹。
此外，此类轨迹分析还隐含一大假设，即只要测序得到的细胞样本数量足够大，那么单一时间点所捕获的细胞“快照”能够涵盖细胞所有可能的状态。但是此假设对于某些生物学过程来说显然有误。例如，细胞群体对外部刺激的响应过程就无法使用单一时间点的“快照”衡量。虽然有一些轨迹分析工具可以结合多个时间点的单细胞数据[8,9]，但大部分主流工具仍然只能针对单一时间点进行轨迹分析。
最后，轨迹分析的一大缺陷还在于无法给出细胞群体演变的方向。于是本文作者在这个方向上做了努力，在一项已发表的工作中提出了RNA速度(RNA velocity)这一概念。简单来说，RNA速度衡量的是前体RNA和成熟mRNA之间的比例，如比例达到平衡态则基因表达在下一时刻可能不发生变化；如果前体RNA数量远大于成熟mRNA，则下一时刻基因有可能上调表达，反之则下调。利用这一特性可以预测该时间点每个细胞在下一时刻的状态，进而推断细胞群体演变的方向。（小师弟云：RNA速度也是一个很酷的概念，相应的分析工具叫做velocyto，大家可以在这里下载安装尝试一下：http://velocyto.org/。另外一个更新的基于RNA速度的工具叫scVelo，可从这里获得：https://scvelo.readthedocs.io）。
7. 单细胞转录组分析在未来将走向何方？
目前为止，单细胞转录组数据分析工具的主要努力方向在于解决数据的高噪声和高稀疏性。今后，随着测序技术的不断完善，研究的重心也会转移到对其他类型分子的测量和分析上：如DNA甲基化、染色质可及性以及蛋白质丰度等。此外，空间转录组也是未来的一大热门方向。如何整合这些不同的数据类型将会是计算生物学领域所要面对的挑战。
5
点评
本文是单细胞转录组分析领域的一篇极具指导意义的综述文章，基本涵盖了单细胞数据分析的所有重要步骤以及每一步蕴含的假设和可能踩到的“坑”，值得一读再读。同时，文中列出的许多分析工具对于广大科研人员的数据分析也是相当实用，大家可以参考我文末附上的参考文献。当然有时间的同学可以品读一下原文，由于推文篇幅所限，无法一一展现其中精华内容。
，我们下期再见啦
！
6
参考文献
1.Svensson,Valentine. Droplet scRNA-seq is not zero-inflated. NatureBiotechnology 38.2 (2020): 147-150.
2.Fan, Jean, et al. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nature methods 13.3 (2016): 241-244.
3.Eling, Nils, et al. Correcting the mean-variance dependency for differential variability testing using single-cell RNA sequencing data. Cell Systems 7.3 (2018): 284-294.
4.Hafemeister, Christoph, and Rahul Satija. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biology 20.1 (2019): 1-15.
5.Risso, Davide, et al. A general and flexible method for signal extraction from single-cell RNA-seq data. Nature communications9.1 (2018): 1-17.
6.Pons, Pascal, and Matthieu Latapy. Computing communities in large networks using random walks. International symposium on computer and information sciences. Springer, Berlin, Heidelberg, 2005.
7. Cao, Junyue, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566.7745 (2019):496-502.
8.Schiebinger, Geoffrey, et al. Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming. Cell 176.4 (2019): 928-943.
9.Tran, Thinh N., and Gary D. Bader. Tempora: cell trajectory inference using time-series single-cell RNA sequencing data. PLoS computational biology 16.9 (2020): e1008205.
注：详细信息请参见原文！
https://www.nature.com/articles/s41592-021-01171-x