来源 | Bioart
撰文 | 伊凯
责编 | 兮
修改 | 生信宝典
自北京大学汤富酬教授（当时为英国剑桥大学格登研究所(Gurdon Institute) Azim Surani实验室博士后）等人于2009年在Nature Methods上发表首个单细胞测序（single cell sequencing）方案以来【1】，这项革命性技术已历经十年的飞速发展；分子生物学、微流控（microfluidics）技术和纳米技术等关联技术的长足进步催生了数十种全新的单细胞测序方案，使测序细胞数目呈现指数级增长 （生信宝典注：指数级增长的转折点是郭国骥老师的工作）（下图）【2】。同时，通过谷歌搜索趋势分析可以发现，对单细胞测序这一词条的相对搜索频率在全球范围内一直呈稳定上升趋势，甚至在2018年超过了同样仅有十余年应用史的重要分子生物学测序方法——染色质免疫共沉淀测序（ChIP-seq）（下图）。
更重要的是，已在常规混合组织测序（bulk sequencing）分析中得到普遍应用的多组学测序技术和单细胞测序技术的成功结合使得在单细胞层面获取DNA甲基化、组蛋白修饰、染色质可接近性等转录组之外的信息从而精确标定细胞状态与命运成为可能（下图）【3】。不过，尽管近年来出现了几项同时测量单细胞多模态信息的技术，但在绝大多数情形下，一次单细胞测序实验往往仍只能覆盖一类单模态信息 （模态应该指：转录组、基因组、修饰组吧）。
因此，在单细胞测序数据量井喷和多模态信息独立并发产生的背景下，数据整合就成为了一个亟待解决的核心挑战。具体而言，对单细胞测序数据的整合可分为两个层次：1）如何将同一组细胞类型基于不同技术平台或由不同批次实验产生的测序数据进行修正与合并？；2）如何将同一组细胞类型的多组学特征信息整合从而完整描绘细胞状态？【4】。显然，这两项挑战在常规混合测序分析中往往也存在，例如RNA测序中的所谓“批次效应”（batch effect）；不过，其解决方案却无法直接迁移至单细胞测序分析中，原因主要在于组织混合测序中对属于同一生理/病理状态的样本的细胞类型组成不变假设在单细胞测序中并不成立。
作为对上述挑战的重要回应之一，来自纽约基因组中心（New York Genome Center）的Rahul Satija课题组于去年4月在Nature Biotechnology上发表了题为Integrating single-cell transcriptomic data across different conditions, technologies, and species的研究论文，提出了以典型关联分析（canonical correlation analysis, CCA）方法为核心的单细胞测序数据整合方案【5】。典型关联分析是一类被广泛运用至图像分析、信号处理和基因组学等领域的经典统计学技术，其基于凸优化和特征值分解方法搜索高维度特征的线性组合以使两组或多组数据间具有最小的整体相关性从而实现数据降维和数据内涵结构的准确捕捉【6】。在该案例中，研究人员假设具有同样细胞类型的两个独立的单细胞RNA测序数据集在本质上可视作对同一（批）转录组的两次抽样，因而必然具有较高的数据整体相关性。因此，对其进行CCA分析不仅使两组数据矩阵从对应于全基因组数千个基因的高维度数据被压缩至对应于前二十项典型关联组分的低维度；更重要的是，这一线性映射还使得去除了批次效应之后的对应于细胞状态的生物学特征信息得以特异性地保留。如下图所示，基于处理前的原始单细胞表达谱所进行的t-SNE降维分析表明单细胞的聚类特征主要由批次效应决定，而经过处理后的不同批次单细胞则很好地交融在了一起，其聚类特征仅由真实的生物学状态所决定。
值得注意的是，在同期的Nature Biotechnology上，来自英国剑桥大学的John Marioni课题组提出了另一项基于截然不同的统计学方法的单细胞测序数据整合方案【7】。不同于以数据整体关联最大化为核心思想的CCA方法，该课题组在近邻（nearest neighbor, NN）算法的基础上提出了互近邻（mutual nearest neighbor, MNN）搜索方法，通过度量两组独立单细胞测序数据集中的数据点间的表达谱余弦距离（cosine distance）找出两个数据集中所有存在的唯一互相最接近的数据点配对；进而，作者假设这些具有互相近邻对应关系的单细胞之间的表达谱差异仅由批次效应造成，故可利用其平均差异计算出批次效应量，在其中一个数据集中减去这一值后即得到了修正后的可直接整合和对比的两组数据集（下图）。
北大汤富酬教授和文路研究员在同期的评论文章中对上述两项工作给予了高度评价，强调了对大规模单细胞测序数据进行有效整合的重要性【8】。
2019年6月6日，Satija课题组在Cell上发表了题为Comprehensive Integration of Single-Cell Data的长文，报道了基于全面更新的统计模型的单细胞测序数据整合方案，并利用这一技术成功实现了单细胞转录组、蛋白组、表观组和空间信息数据集之间的信息迁移。这一最新方法已被整合至基于R语言的单细胞测序分析库Seurat的最新发行版中。
在该研究中，作者提出利用两个或多个独立单细胞测序数据集中存在的少部分属于同一细胞类型的数据点对整个数据集进行“锚定”进而整合。具体而言，该新算法综合了前述两类分别基于CCA降维和MNN搜索互近邻的方案，具有更高稳健性和整合准确率：首先，该算法以在低维空间中匹配成功的互近邻单细胞对为锚点（下图A、B、C），通过考察每个锚点配对在参考数据库和被测数据集中的近邻的交叠率对相应的锚点的可靠性进行评价，形成一个锚点稳健性分数向量（下图D）；其次，该算法测量被测数据库中的非锚点单细胞与锚点单细胞的距离，在乘以锚点稳健性分数后得到一个锚点权重矩阵，即被测数据库中的每一个非锚点单细胞受到每一个锚点单细胞在数据修正上的影响程度；最后，将被测数据库中的每个数据点的表达谱减去该数据点对应的锚点权重值与锚点差异值的乘积，即得到细胞类型特异的批次效应修正后的可直接与参考数据集进行整合的表达谱（下图E）。
为了验证上述数据修正与整合算法相较于现有方法（尤其是CCA和MNN）的优势，研究人员分析了这些方法对一个包含8个独立亚数据集的胰岛单细胞RNA测序数据集的整合效果。经由UMAP降维和可视化，研究人员发现新算法（Seurat v3）成功将来源于不同数据集的单细胞聚类至同一聚落中，而之前的CCA算法（Seurat v2）和MNN算法，亦或是另一算法scanorama都存在相当程度的整合后批次效应（下图）。
为了定量分析上述由聚类模式展示的不同算法在数据修正和整合上的表现差异，研究者试图从两个角度对其进行考察：一是整合后来源于不同批次的数据点的混杂程度，二是整合后同一批次中的数据点的原有聚类模式的保留。这两个参数分别表征了数据整合对批次效应的去除程度和对真实生物学效应的保存程度。如下图所示，在两项不同的数据整合任务中，新算法都具有最优的综合表现。
在已确证新整合算法的显著优越性之后，研究人员通过三个具体案例分别展示了其在单细胞转录组、表观组、蛋白组和空间位置信息数据集之间进行信息迁移的强大本领。在一个例子中，研究人员希望通过锚定一对具有匹配细胞类型的单细胞RNA测序数据集和单细胞染色质可及性测序（scATAC-seq）数据集从而将由转录组特征决定的细胞类型聚类结果转移至ATAC测序数据中，这一操作背后的逻辑在于：尽管ATAC测序可以提供基因表达上游的关于基因调控元件开放程度的全新信息，然而其目前还尚不能单独达到足以媲美基于单细胞RNA测序对细胞进行非监督聚类的效果。因此，研究人员的计划是首先将已经由基因表达谱识别出的细胞类型标签直接迁移至不存在这一信息的ATAC测序数据点上，再进一步利用ATAC测序结果本身对细胞聚类结果进行微调，或寻找新的无法由转录组信息单独识别的独立细胞群。由于ATAC测序不是对单个基因的表达水平进行定量，因此为使得两个数据集之间可成功整合，研究人员首先利用ATAC测序得到的关于基因启动子和增强子等调控元件的可及性水平间接推算出相应基因的“激活程度”，然后便可如前所述对两个数据集进行准确锚定，随后再直接利用代表被测数据集中非锚点单细胞受锚点单细胞影响程度的锚点权重矩阵与参考数据集中的细胞聚类标签矩阵相乘，实现聚类标签的转移，即得到被测数据集（ATAC测序数据集）中单细胞的基于基因表达谱的细胞类型。依照这一方法，研究者发现新算法不仅取得了显著高于现有算法的表现（下图C），更重要的是通过其实现的整合操作导致了多个新的细胞亚群的发现（下图E）。这一结果表明，通过在测量不同细胞模态信息的单细胞测序数据集之间共享细胞类型标签，便能够充分利用多组学信息在对细胞状态进行精确刻画上的优势。
如果以离散变量形式存在的细胞聚类标签可以经由这一全新的数据集锚定整合算法在不同数据集间进行迁移，那么其它以连续变量形式存在的数据，例如单细胞膜蛋白表达量或基因表达水平本身，是否也可以实现这一操作呢？对这个问题的答案将显著影响现有的单细胞测序数据的分析模式。一个例子是现存最大的人骨髓单细胞测序数据集——由人类细胞图谱计划（Human Cell Atlas, HCA）提供的近30万个健康人群骨髓单细胞基因表达谱——缺失了膜蛋白表达水平的相关信息，这显然阻碍了对这一以免疫细胞群为核心的大规模数据集的深入挖掘分析。为填补这一数据空白，研究人员对3万多个人类骨髓细胞进行了单细胞全转录组RNA测序和25个膜蛋白表达水平的定量分析，构建了一个具有较好代表性的人骨髓细胞参考数据集。接着，研究者重复了前述数据集锚定、整合及信息迁移步骤（在数学上，即以膜蛋白表达水平的25维连续变量矩阵代替了细胞类型标签的二元矩阵），成功为HCA中的数十万个本不具有膜蛋白表达水平信息的单细胞填补了这一空白。对这一数据迁移结果的可靠性的信心来源于基于原本3万多个膜蛋白表达水平信息齐备的单细胞的测试结果，其中绝大部分蛋白的预测表达量均与真实值呈显著强相关（下图A）。另外，研究者还发现在数据锚定和整合过程中，并不需要利用全转录组表达谱，实际上在仅利用200个基因时就已经具有相当好的蛋白表达预测结果（下图B）。这些结果表明新算法不仅能实现独立单细胞测序数据集间的离散型数据迁移，还能实现连续型数据迁移，这无疑极大地扩展了数据整合的应用场景。
最后，研究人员还探究了利用数据整合算法将单细胞RNA测序的全转录组基因表达水平信息迁移至具有单细胞空间信息但仅有1千个左右基因表达值的RNA荧光原位杂交（RNA FISH）数据集中，从而实现全空间尺度的全转录组定量分析的可能性（下图A）。由之前的结果可知，两个独立单细胞数据集之间的锚定往往仅需要极小部分的基因参与，因此FISH数据集中基因表达维度较低的问题并不会显著影响数据整合结果。事实上，经由前述相同的分析步骤，研究人员成功地将全转录组3万多个转录本的表达水平迁移至了FISH数据集中的单细胞，从而在全空间尺度上恢复了整个基因表达谱的完整信息。这一结果的可靠性可由同时存在于两个数据集中的基因在FISH中的真实分布和预测分布的显著相似性得出（下图B）。
总之，Satija课题组开发的这一普遍适用于各类单细胞测序数据集的锚定整合方案，通过精妙的算法设计和严谨的独立测试，在消除批次效应、保留真实生物学效应、离散及连续型多模态信息完整迁移等方面均取得了优异表现。在当今单细胞多组学测序数据量爆发式增长但对于单个实验室而言测序成本仍然较高、测序流程仍然复杂的背景下，以Seurat为代表的数据整合方法无疑为充分挖掘及联合分析已有数据和新生数据中的生物学靶标提供了强大助力。
原文链接：
https://doi.org/10.1016/j.cell.2019.05.031
参考文献
1.Tang, F. et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat. Methods (2009). doi:10.1038/nmeth.1315
2.Wen, L. & Tang, F. Single-cell sequencing in stem cell biology. Genome Biology (2016). doi:10.1186/s13059-016-0941-0
3.Shema, E., Bernstein, B. E. & Buenrostro, J. D. Single-cell and single-molecule epigenomics to uncover genome regulation at unprecedented resolution. Nature Genetics (2019). doi:10.1038/s41588-018-0290-x
4.Stuart, T. & Satija, R. Integrative single-cell analysis. Nat. Rev. Genet. 20, 257–272 (2019).
5.Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).
6.Hardoon, D. R., Szedmak, S. & Shawe-Taylor, J. Canonical correlation analysis: An overview with application to learning methods. Neural Computation (2004). doi:10.1162/0899766042321814
7.Haghverdi, L., Lun, A. T. L., Morgan, M. D. & Marioni, J. C. Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nat. Biotechnol. 36, 421–427 (2018).
8.Wen, L. & Tang, F. Boosting the power of single-cell analysis. Nat. Biotechnol. (2018). doi:10.1038/nbt.4131
