单细胞测序技术的出现使得人们能够以前所未有的精度观测细胞。然而,单次单细胞转录组测序(scRNA-seq)实验难以捕获所有细胞和基因的信息,单个模态的单细胞数据无法详细阐释细胞状态和系统变化,单细胞数据的整合分析旨在解决这两类问题。整合不同来源的scRNA-seq数据,可以收集完整的细胞类型,为构建细胞图谱提供强大助力;整合多个模态的单细胞数据,可以研究模态间因果关系和基因调控机制。数据整合方法的开发与应用帮助充分挖掘单细胞数据的丰富性和相关性,发现有意义的生物学变化。基于此,本文综述了多源scRNA-seq数据整合和单细胞多模态数据整合的基本原理、方法和应用,并讨论了现有方法的优势和不足,最后对未来的发展前景予以展望。
引用本文: 潘多, 李华梅, 刘宏德, 孙啸. 单细胞数据的整合方法综述. 生物医学工程学杂志, 2021, 38(5): 1010-1017. doi: 10.7507/1001-5515.202104073 复制
引言
细胞是生物体的基本组成单位,每个细胞都是独一无二的。单细胞测序从单个细胞水平上获取细胞信息,解决了传统测序技术基于群体细胞测序而掩盖细胞间异质性的难题,把研究精度从群体细胞精确到单个细胞,有助于深入探索细胞的异质性和生命活动。单细胞转录组测序(single-cell RNA sequencing,scRNA-seq)是当前应用最广泛的单细胞测序技术,它对单个细胞的mRNA进行反转录、扩增和高通量测序,根据转录谱的相似性,可以分辨不同的细胞类型,甚至揭示新的细胞类型[1]。构建细胞图谱是scRNA-seq技术发展的重要应用,是研究基础生物学、衰老、疾病、治疗反应的基本前提[2]。然而,单次scRNA-seq实验难以捕获一个组织或器官中所有的细胞和基因,获得的数据无法提供完整的细胞信息。同时,测序深度和规模、生物噪声和技术噪声都将限制单个scRNA-seq数据集的可挖掘性。因此,整合分析多个来源的scRNA-seq数据能够有效弥补单个数据集的局限,收集更完整的细胞类型和生物信息,助力生成全景细胞图谱。多源scRNA-seq数据整合是指集成多个不同来源(实验、平台、样本、组织等)的scRNA-seq数据集,校正系统差异使细胞基于细胞类型而聚集,生成高质量的大数据集,如图1所示。
scRNA-seq迅速发展的同时,聚焦于其他模态的单细胞技术也在不断地开发和应用中,人们可以获取同一细胞的不同模态数据,如基因组测序、DNA甲基化、染色体构象、染色质可及性、蛋白质丰度等。目前,大多数单细胞技术只是单纯研究细胞某一模态生物分子的变化,提供细胞异质性的局部景观,不足以深入理解细胞状态和功能。因此,整合分析多个模态的单细胞数据能够从不同层次的模态特征解释细胞所处的状态、展现细胞内含的生物逻辑,这将弥补单一模态的片面性,使人们能够从各种生物分子的组成、结构、功能等方面深刻了解和认识细胞。单细胞多模态数据整合是指分析单细胞多个模态组学的数据,发现模态间交互关系,挖掘基因调控机制,增强细胞类型的可解释性,如图1所示。
单细胞数据整合分为多源scRNA-seq数据整合和多模态数据整合,前者侧重某一组织内细胞类型的完整性,后者更侧重调控机制的挖掘,这两者分别从横向和纵向的角度提供细胞的全景图和近景信息,以多视角、多视图来解释细胞的异质性和身份,为生物系统研究提供重要综合性参考。近年来,研究人员开发了多种方法来实现单细胞数据的整合,作为数据分析的核心,发展创新且有效的数据整合方法为充分挖掘单细胞数据、提取有意义的生物学信息提供了强大助力。本文将对多源scRNA-seq数据整合和单细胞多模态数据整合的基本原理、方法和应用进行综述,并根据方法的特色优势和不足提出应用指南,并对未来的发展前景予以展望。
1 多源scRNA-seq数据整合
1.1 多源scRNA-seq数据的整合方法
单细胞转录组数据通常来自多个实验,由于技术性因素(扩增方式、测序平台、实验人员等)、生物性因素(个体、处理、物种等)以及数据自身性质(细胞数量、测序深度)的差异,数据集在合并时会呈现明显的批次效应,表达量具有系统偏差,这将混淆真正的生物逻辑变化。因此,多源scRNA-seq数据整合的目标就是校正系统偏差,使细胞基于真实的生物类型或状态聚集在一起,而不受来源的影响[3]。
scRNA-seq数据通常以基因表达矩阵的形式呈现,矩阵的行为基因、列为细胞,因此scRNA-seq数据整合本质上就是表达矩阵的整合。在整合前,一般会对每个数据集执行相同的预处理,包括质控、归一化、标准化、特征选择,选取一定数量的高变异基因(highly variable genes,HVGs)。然后,通常是取每个数据集HVGs的交集用于整合,使得基因数量保持一致。整合时,为了统一度量细胞间的距离,需要将待整合的数据集放置在同一个空间下。由于在高维的基因表达空间下计算复杂,研究人员通常使用降维策略将所有数据映射到一个共享的低维空间,前提是所有数据集至少共享一种相同细胞类型,并且默认假设来自不同数据集但具有相同细胞类型或生物状态的细胞互相靠近,这些细胞对反映了不同数据集间的对应关系[4]。本文将多源scRNA-seq数据整合方法分为配准和图聚类两种策略。配准类方法的核心思想是在共享空间下,以一个数据集作为基准,另一个数据集作为目标,根据两个数据集相似细胞对间的距离,以及相似细胞对与目标细胞的关联程度,计算目标细胞特异的校正向量,指导目标数据集向基准数据集配准;或者,各个数据集均向公共的基准尺度(例如质心)配准,如图2所示。此类方法结果可能输出一个校正的综合表达矩阵,也可能输出校正后的低维嵌入,即配准后每个细胞在共享低维空间中的坐标。图聚类方法的核心思想是针对单细胞进行图的聚类分析,将不同数据集相似的细胞聚类在一起,如图2所示。在共享空间中,以细胞为节点,连接节点形成边,基于细胞间距离为边分配权值,构建加权联合图。接着,利用社区发现算法对图进行聚类,将来源不同但密切连接的细胞聚类在一起,实现整合,结果通常以图的形式输出。有的方法同时运用了配准和图聚类两种策略,在构建联合图后,还对数据集进行了配准。当输出结果为一个校正的基因表达矩阵时,可以进行广泛的下游分析;校正低维嵌入和联合图可用于可视化、聚类和拟时序分析,由于这两种结果没有改变原始表达值,因此不能直接用于差异表达分析、基因调控网络等基因层面的分析。可以使用批次感知的线性混合效应模型进行差异表达分析,或者对未校正的基因表达矩阵进行单独分析和共同比较[5-6]。此外,除了上述配准类和图聚类的整合方法,还有基于深度学习的整合方法。本文总结了近年来主要的多源scRNA-seq数据整合方法,包括整合策略、应用的主要算法、输出结果、程序语言和参考文献,如表1所示。
配准类方法中,多集典型相关分析(multi-set canonical correlation analysis,MultiCCA)是最早被报道的经典方法之一,支持在单细胞基因组学R工具包Seurat 2.0(Satija实验室,美国)中运行使用[7]。MultiCCA首先利用典型相关分析(canonical correlation analysis,CCA)将两个数据集映射到一个低维空间,接着利用动态时间规整(dynamic time warping,DTW)将两个数据集配准到一个共同的基准尺度,输出对齐的低维嵌入。同期,Haghverdi等[4]提出了以相互最近邻(mutual nearest neighbors,MNN)为核心的MNN校正(MNN-Correct),对于基准数据集的每个细胞,确定在目标数据集中k个最近邻(k-nearest neighbors,KNN),对于目标数据集也执行相同的操作;当一对细胞包含在彼此的k个近邻中,则视为一对MNN,代表不同数据集中具有相同细胞类型或状态的细胞对。MNN-Correct在余弦归一化的表达数据上计算细胞间的欧氏距离,得到一组MNN,MNN对的距离差值可以估算批次效应的大小,根据这些MNN对差异向量的加权平均值计算细胞特异的校正向量,将目标数据集向基准数据集配准。由于MNN-Correct在高维空间中运算,计算量巨大,快速MNN(fastMNN)作为MNN-Correct的新版本,增加了主成分分析(principal component analysis,PCA)降维步骤,运行速度提升。
受MNN的启发,在修拉(Seurat)版本三(Seurat V3)、单细胞合并(scMerge)、全景(Scanorama)中也引入了MNN的思想[6, 8-9]。Seurat V3在CCA定义的低维空间中识别MNN,并取名为锚点。为了确认锚点的可靠性,Seurat V3对锚点进行了过滤和评分:首先,返回高维的表达空间查看锚点的KNN条件是否同样满足,排除表达谱差异很大的错误锚点;然后,根据每对锚点细胞的邻域重叠率为每个锚点打分,在后续计算中降低低分数锚点的权重,以提高结果的鲁棒性和准确率。scMerge利用伽马-高斯混合模型识别一组稳定表达基因(stably expressed genes,SEGs);在PCA空间下识别不同数据集的相互最近集群,其中一定比例的最接近集群质心的细胞集称为伪重复,即不同数据集的同一类细胞;利用三步删除不必要变异模型(remove unwanted variation, 3-step,RUVⅢ)计算SEGs在伪重复中的表达差值,删除来源影响并配准数据集[10]。Scanorama将MNN的识别从两个数据集推广至多个数据集中,对于每个数据集的细胞,KNN的搜寻范围扩展为其余数据集的所有细胞,从而获取所有数据集中转录谱相似的MNN匹配,通过这些匹配将数据集拼接起来形成一张“全景图”。类似于MNN-Correct的校正计算,每个细胞的校正向量由其近邻范围内MNN对差异向量的加权平均得到,然后根据MNN细胞的占比确定每个数据集的合并顺序,依次合并到全景图中。Scanorama使用随机奇异值分解(randomized singular value decomposition,rSVD)进行降维,基于局部敏感哈希算法搜索近似的KNN,这大大加快了运行速度,可灵活应用于大数据集。鸿蒙(Harmony)是一种迭代聚类的整合方法,将所有数据集的细胞向聚类质心进行配准。在PCA空间中,Harmony使用k均值软聚类将相似细胞聚类在一起,并且最大化每个集群内的数据集多样性;在每个集群中,根据集群质心和数据集特异质心,计算细胞特异的校正因子来校正每个细胞;迭代运行以上两个过程,相当于迭代消除数据集来源对细胞聚类分配的影响,直到细胞的低维嵌入收敛到精确值,即完成整合,输出校正的低维嵌入[5]。
批次平衡k近邻(batch balanced k nearest neighbours,BBKNN)、样本网络聚类法(clustering on network of samples,Conos)和基因组实验关系关联推断(linked inference of genomic experimental relationships,LIGER)是图聚类的整合方法[11-13]。BBKNN首先在PCA空间中识别每个细胞在数据集内和不同数据集间的KNN并连接起来,然后按照一致流形近似和投影(uniform manifold approximation and projection,UMAP)的方法将近邻距离转换为连接分数,为每个连接赋予相应的权重,输出加权联合图用于聚类。Conos在低维空间中应用旋转空间和近邻(KNN和MNN)的策略,将数据集两两比较,建立数据集间和数据集内的边,用于构造联合图,其中根据细胞间相似性确定边的权重,并且降低数据集内的边相对于数据集间的边在图中的贡献;然后使用社区发现算法来获得聚类。LIGER应用综合非负矩阵分解(integrative non-negative matrix factorization,iNMF)进行降维,识别数据集特异和共享的两组因子来定义每个细胞的低维嵌入[14]。在因子空间下,计算细胞在数据集内和数据集间的KNN,并连接成对的相似细胞,即具有相似因子载荷模式的细胞,从而构建一个共享因子邻域图,用于社区发现算法聚类。聚类后,LIGER还对每个联合集群的因子载荷进行分位数归一化(quantiles normalization,QN),输出校正的低维嵌入。
近年来,研究人员将深度学习技术引入生物学领域,帮助处理生物数据资料,解决生物信息学问题,在数据整合方面也不断有新的进展。深度学习类方法的核心思想是利用深度网络学习数据结构和配准关系,处理优化问题,以克服批次效应。最大均值差异残差网络(maximum mean discrepancy and residual nets,MMD-ResNet)应用残差神经网络模型,以最小化基准数据集和目标数据集分布的最大均值差为目标,学习两个数据集之间的映射关系,使得两者在分布上相似[15]。深度MNN(deepMNN)也应用了残差神经网络,该方法利用PCA空间中识别的MNN对训练网络,并以最小化MNN对中的细胞间距、使网络的输出与输入相似为目标,指导网络学习从而消除数据集间的批次效应[16]。单细胞基因表达数据模型(scGen)结合了变分自编码器(variational auto-encoder,VAE)模型和潜在空间向量算法,只适用于已知细胞类型的数据整合。该方法的过程是输入标记细胞类型的数据,经过编码器映射到一个潜在空间中,在此空间下计算不同数据集中相同细胞类型的细胞间的差异向量,然后向其中一个数据集进行配准,再经过解码器映射回原本的高维表达空间,输出一个校正的基因表达矩阵[17]。Wang等[18]利用自编码器和生成式对抗网络(generative adversarial network,GAN)构建了对抗式配对风格迁移网络集成多源单细胞数据集(integration of multiple single-cell datasets by adversarial paired-style transfer networks,iMAP)模型,该模型分两个步骤进行数据整合,首先以保留生物变异、去除批次影响为优化目标,在自编码器中重建细胞表达谱,第二步在识别的MNN上训练GAN以正确匹配共享细胞类型的分布。其中,iMAP还引入了随机游走的策略来扩展MNN列表,从而更好地覆盖共享细胞类型的全部分布,有助于GAN的训练。
1.2 多源scRNA-seq数据整合的应用
许多整合方法现已应用于scRNA-seq数据分析中,去除批次效应,鉴定新的细胞类型,绘制单细胞图谱,并在肿瘤、脑科学、病毒等研究领域取得了一定成果。例如,Zhang等[19]利用Harmony整合了一位肝癌患者多个组织的、两种测序技术获得的免疫细胞,综合不同测序技术的优势能够绘制出更高分辨率的肝癌免疫图谱,通过联合分析还识别出了罕见的细胞群,对肝癌的免疫治疗具有指导意义。Liao等[20]利用Harmony和fastMNN整合来自三个不同捐赠者的肾脏组织单细胞转录组数据,鉴定出三个近端小管细胞的亚型和两个集合管细胞的亚型,为肾小管细胞的精确分类和相关疾病研究提供了重要的参考。Trujillo等[21]建立了一个皮层类器官模型用于模拟人脑早期发育,用Seurat V3整合四个培养时间点下所测序的类器官单细胞数据,经细胞注释,类器官发育过程中的细胞类型变化证明该模型具有功能性。这些整合方法在新型冠状病毒(severe acute respiratory syndrome coronavirus 2,SARS-CoV-2)的研究中也提供了有效助力。Qi等[22]收集了13种人体组织的单细胞基因表达数据,通过Harmony整合数据,确定了三个具有与SARS-CoV-2受体血管紧张素转化酶2最相似表达模式的候选基因,其可能是协助SARS-CoV-2入侵人体的共受体编码基因,有助于制定相应的干预策略。Zhang等[23]分别用Seurat V3和LIGER整合了胃和回肠数据集、结肠数据集,综合分析SARS-CoV-2受体血管紧张素转化酶2和跨膜丝氨酸蛋白酶2的共表达模式,揭示了SARS-CoV-2在消化系统中的潜在传播途径。Seurat V3还被用于整合来自健康者和不同严重程度的新冠肺炎患者的细胞数据,绘制出新冠肺炎患者的肺泡灌洗液免疫细胞图谱和免疫反应图谱[24-25]。除此以外,Tran等[26]评估了14种方法在多批次、大数据集、不同测序技术、细胞类型差异大等数据情形下的整合效果,Harmony、LIGER和Seurat V3是总得分前三的方法。
2 单细胞多模态数据整合
2.1 单细胞多模态数据的整合方法
2019年,《自然-方法》杂志将单细胞多模态组学选为“2019年度技术”,揭示了单细胞领域联合多个模态和组学进行测量分析的重要发展趋势[27-28]。单细胞多模态技术在一次实验中测量同一个细胞多个模态数据,表征细胞内不同层面的生物信息:DNA和RNA测序从底层获取基因序列和表达信息,检测基因动态变化,蛋白组和表观组从表层反映基因表达调控机制和分子性状的变化规律。多模态研究为后续的整合分析提供了丰富的数据,然而复杂的数据类型和数据特征也带来了计算上的挑战。例如,转录组数据通常是基因表达矩阵;DNA甲基化提供甲基化信号值矩阵,包含甲基化位点和甲基化水平信息;蛋白组提供蛋白表达矩阵;染色质可及性提供了开放区域的信息,称为峰矩阵。对于这些描述同一细胞不同模态的数据,其整合的目标是揭示不同模态之间的关联性,并且更详细地描述细胞状态和基因调控机制[29]。因此,一个关键问题就是如何将不同模态的信息联系到一起。第一种策略是基于已有假设,将其他模态信息关联到基因层面,进行基因级矩阵转换[28]。例如,对于染色质可及性,可以将峰矩阵转成基因活跃度矩阵,矩阵的行为基因、列为细胞,数值为基因活跃度得分,由跨越基因体和启动子区域(通常上游2~3 kb)的峰读段数相加得到[30];对于DNA甲基化数据,利用基因体的非CG甲基化推断基因表达,生成基因级甲基化数据矩阵[31]。经过这样的转换,其他模态的基因级矩阵与scRNA-seq矩阵的整合等同于两个不同来源scRNA-seq数据的整合,即本文第1.1节讨论的情况。第二种策略是直接从模态数据的本质、特征、关系出发,输出结果是模态信息的综合体,用来揭示模态间的因果关系。根据不同方法模型的处理,结果以不同形式呈现。
传统转录组测序领域早已有许多成熟有效的多组学整合算法和工具,虽然整合原理同样适用于单细胞多模态数据,但还需要适当的调整和实验验证[32-33]。本文总结了适用于单细胞多模态数据整合的主要方法,包括整合的数据类型、主要算法、程序语言和参考文献,如表2所示。
Seurat V3、Conos和LIGER应用了基因级矩阵转换的策略,将单细胞染色质开放区转座酶可及性测序(single-cell assay for transposase-accessible chromatin using sequencing,scATAC-seq)的峰矩阵转换为基因活跃度矩阵,与scRNA-seq矩阵相整合,LIGER还实现了基因级甲基化数据与scRNA-seq数据的整合,原理如本文第1.1节所述[6, 12-13]。
聚类通常是单细胞数据下游分析的第一步,对于scRNA-seq和scATAC-seq都有相应的聚类方法识别细胞类型,而协同两者信息能更好地解释细胞类型。Duren等[34]从聚类的角度,提出耦合scRNA-seq和scATAC-seq两个聚类过程的耦合非负矩阵分解(coupled nonnegative matrix factorization,Coupled NMF)模型,首先向模型中输入两个数据矩阵,然后求解耦合聚类的优化问题,输出两个数据的共同聚类结果以及对应的峰与基因配对。这一模型在反卷积与耦合聚类(de-convolution and coupled-clustering,DC3)模型中得到了升级,加入了群体细胞数据的反卷积化步骤,以获得细胞亚群特异的数据,帮助改善单细胞耦合聚类的结果[35]。多组学因子分析(multi-omics factor analysis,MOFA)是一种因子分析模型,适用于单细胞DNA甲基化与RNA数据的整合分析[36]。向MOFA模型中输入不同模态的数据矩阵,经过矩阵分解,结果输出一组因子,这些因子代表了驱动不同模态数据异质性的因素;下游分析中,可以得到影响因子贡献度的重要基因和甲基化信息,并对所有样本进行可视化、聚类、富集分析。流行对齐表征实验关系(manifold alignment to characterize experimental relationships,MATCHER)应用了流形学习的降维策略,通过高斯过程潜变量模型(Gaussian process latent variable model,GPLVM)将不同模态数据(基因表达、DNA甲基化、染色质可及性)映射到一维流形空间,每种模态数据都可以由一组伪时间值表示,并对伪时间值进行QN以实现统一度量,通过比较和分析该值可以研究多个模态之间的相关性和潜在调控机制[37]。
2.2 单细胞多模态数据整合的应用
有效的计算方法为充分挖掘不同模态的单细胞数据提供了机会,本节讨论单细胞多模态数据整合方法的重要应用,包括揭示细胞异质性、发现模态间交互关系和推断基因调控网络等。Welch等[37]利用MATCHER研究小鼠胚胎干细胞和人类诱导多能干细胞的转录组和表观基因组之间的相关性,结果揭示单细胞基因表达和DNA甲基化、染色质可及性、组蛋白修饰之间具有共同的变异模式,轨迹分析显示了细胞从多能性到分化启动状态的变化。Argelaguet等[36]也分析了小鼠胚胎干细胞的转录组和甲基化数据,MOFA分析结果同样揭示细胞分化过程中转录组和甲基化水平具有协同变化。Argelaguet等[38]应用三重组学测序技术获取小鼠胚胎细胞的单细胞核小体、甲基化和转录组三重模态数据,通过MOFA整合分析揭示了原肠胚形成过程中具有谱系特异的表观遗传模式和标志基因,并在三个模态推断因子中观察到了细胞间的异质性。因此,MOFA从最初的转录组与甲基化二重模态分析拓展到了多重模态分析,应用范围扩大。Lake等[39]通过训练一个梯度提升回归模型,在成人大脑细胞的转录组和表观基因组之间建立映射,发现了驱动细胞异质性的调节元件和转录因子,为研究大脑的复杂过程提供了新的思路。LIGER被用于整合scRNA-seq和DNA甲基化数据,联合定义了小鼠皮质细胞类型,揭示了细胞类型特异性的表观基因组调控机制[13]。研究单个细胞不同模态之间的相关性可帮助研究人员深层次理解基因调控网络,提高细胞类型分类的准确性和可解释性,使全面探索细胞身份和行为成为可能。
3 总结与展望
单细胞数据的整合分析提供了更为全面且深入的见解,帮助从不同层面剖析细胞类型和状态,深入挖掘调控机制。近年来,单细胞数据整合方法的开发与应用方面已经有了显著成果。
多源scRNA-seq数据整合方法中,Seurat系列和MNN-Correct系列的方法假设数据集之间的差异完全源于技术性变化,适用于处理细胞类型相似的数据集。在处理细胞组成差异大的情况下,LIGER相对适用,因为iNMF可以保留数据集间差异,也可以识别相似之处。对于多数方法,基准数据和目标数据的选择会影响整合效果;当整合两个以上的数据集时需要迭代整合过程,整合次序也将影响结果。Scanorama应用全景拼接的原理,避免了输入数据集的顺序对整合效果的影响。Harmony对所有数据执行迭代聚类和校正,对数据集次序不敏感。Harmony和BBKNN运行速度快,均适用于处理大数据集。深度学习类方法在复杂数据中具有可扩展性和高性能,而对于量小的数据集,此类方法的性能较差,不利于网络训练。有些深度学习类方法只能实现有监督或半监督的整合,例如scGen需要输入有细胞类型标签的数据集。deepMNN通过最小化MNN对的细胞间距,促进网络同时消除多个数据集的批次效应,因此能够实现一步整合。此外,深度学习类方法使用图形处理器来加速计算,同样适用于处理大数据集。
单细胞多模态数据整合方法中,Seurat V3、LIGER和MATCHER能够实现基因表达与DNA甲基化、染色质开放程度、组蛋白修饰等多个模态的整合分析,完整地建立起单细胞转录组和表观基因组之间的调控关系。DC3还实现了单细胞与群体细胞数据的联合分析,有效改善细胞的聚类结果。MOFA是相对流行的方法,从因子分析的角度解析数据异质性,使得分析结果也更具解释性。
尽管目前的单细胞数据整合方法具有良好的应用价值,但仍存在以下几个方面的挑战。首先,合适的整合方法是有效整合分析的基础,其选择依赖于对不同整合方法的性能测试,所以形成成熟的评测基准或评测系统至关重要,能够指导数据整合工作,减少方法本身对整合结果的影响。第二个挑战在于整合方法通常基于不同的计算平台开发,方法的推广和使用受限于用户对环境平台的偏好、数据格式和预处理步骤等因素,因此未来需要向跨平台分析、数据转换共享、流程化运行的趋势发展。第三个方面,单细胞数据规模日益扩大,以及适应“人类细胞图谱计划”(Human Cell Atlas,HCA)的需求,如何高效处理大规模数据集是亟待解决的问题,而大数据集恰好与深度学习的数据驱动性质相适应,这可能是未来可以继续推进研究的方向。随着新整合方法不断突破和生物信息学的发展,单细胞数据的整合分析定会提供更全面的细胞视角,帮助解决生物学难题。
利益冲突声明:本文全体作者均声明不存在利益冲突。
引言
细胞是生物体的基本组成单位,每个细胞都是独一无二的。单细胞测序从单个细胞水平上获取细胞信息,解决了传统测序技术基于群体细胞测序而掩盖细胞间异质性的难题,把研究精度从群体细胞精确到单个细胞,有助于深入探索细胞的异质性和生命活动。单细胞转录组测序(single-cell RNA sequencing,scRNA-seq)是当前应用最广泛的单细胞测序技术,它对单个细胞的mRNA进行反转录、扩增和高通量测序,根据转录谱的相似性,可以分辨不同的细胞类型,甚至揭示新的细胞类型[1]。构建细胞图谱是scRNA-seq技术发展的重要应用,是研究基础生物学、衰老、疾病、治疗反应的基本前提[2]。然而,单次scRNA-seq实验难以捕获一个组织或器官中所有的细胞和基因,获得的数据无法提供完整的细胞信息。同时,测序深度和规模、生物噪声和技术噪声都将限制单个scRNA-seq数据集的可挖掘性。因此,整合分析多个来源的scRNA-seq数据能够有效弥补单个数据集的局限,收集更完整的细胞类型和生物信息,助力生成全景细胞图谱。多源scRNA-seq数据整合是指集成多个不同来源(实验、平台、样本、组织等)的scRNA-seq数据集,校正系统差异使细胞基于细胞类型而聚集,生成高质量的大数据集,如图1所示。
scRNA-seq迅速发展的同时,聚焦于其他模态的单细胞技术也在不断地开发和应用中,人们可以获取同一细胞的不同模态数据,如基因组测序、DNA甲基化、染色体构象、染色质可及性、蛋白质丰度等。目前,大多数单细胞技术只是单纯研究细胞某一模态生物分子的变化,提供细胞异质性的局部景观,不足以深入理解细胞状态和功能。因此,整合分析多个模态的单细胞数据能够从不同层次的模态特征解释细胞所处的状态、展现细胞内含的生物逻辑,这将弥补单一模态的片面性,使人们能够从各种生物分子的组成、结构、功能等方面深刻了解和认识细胞。单细胞多模态数据整合是指分析单细胞多个模态组学的数据,发现模态间交互关系,挖掘基因调控机制,增强细胞类型的可解释性,如图1所示。
单细胞数据整合分为多源scRNA-seq数据整合和多模态数据整合,前者侧重某一组织内细胞类型的完整性,后者更侧重调控机制的挖掘,这两者分别从横向和纵向的角度提供细胞的全景图和近景信息,以多视角、多视图来解释细胞的异质性和身份,为生物系统研究提供重要综合性参考。近年来,研究人员开发了多种方法来实现单细胞数据的整合,作为数据分析的核心,发展创新且有效的数据整合方法为充分挖掘单细胞数据、提取有意义的生物学信息提供了强大助力。本文将对多源scRNA-seq数据整合和单细胞多模态数据整合的基本原理、方法和应用进行综述,并根据方法的特色优势和不足提出应用指南,并对未来的发展前景予以展望。
1 多源scRNA-seq数据整合
1.1 多源scRNA-seq数据的整合方法
单细胞转录组数据通常来自多个实验,由于技术性因素(扩增方式、测序平台、实验人员等)、生物性因素(个体、处理、物种等)以及数据自身性质(细胞数量、测序深度)的差异,数据集在合并时会呈现明显的批次效应,表达量具有系统偏差,这将混淆真正的生物逻辑变化。因此,多源scRNA-seq数据整合的目标就是校正系统偏差,使细胞基于真实的生物类型或状态聚集在一起,而不受来源的影响[3]。
scRNA-seq数据通常以基因表达矩阵的形式呈现,矩阵的行为基因、列为细胞,因此scRNA-seq数据整合本质上就是表达矩阵的整合。在整合前,一般会对每个数据集执行相同的预处理,包括质控、归一化、标准化、特征选择,选取一定数量的高变异基因(highly variable genes,HVGs)。然后,通常是取每个数据集HVGs的交集用于整合,使得基因数量保持一致。整合时,为了统一度量细胞间的距离,需要将待整合的数据集放置在同一个空间下。由于在高维的基因表达空间下计算复杂,研究人员通常使用降维策略将所有数据映射到一个共享的低维空间,前提是所有数据集至少共享一种相同细胞类型,并且默认假设来自不同数据集但具有相同细胞类型或生物状态的细胞互相靠近,这些细胞对反映了不同数据集间的对应关系[4]。本文将多源scRNA-seq数据整合方法分为配准和图聚类两种策略。配准类方法的核心思想是在共享空间下,以一个数据集作为基准,另一个数据集作为目标,根据两个数据集相似细胞对间的距离,以及相似细胞对与目标细胞的关联程度,计算目标细胞特异的校正向量,指导目标数据集向基准数据集配准;或者,各个数据集均向公共的基准尺度(例如质心)配准,如图2所示。此类方法结果可能输出一个校正的综合表达矩阵,也可能输出校正后的低维嵌入,即配准后每个细胞在共享低维空间中的坐标。图聚类方法的核心思想是针对单细胞进行图的聚类分析,将不同数据集相似的细胞聚类在一起,如图2所示。在共享空间中,以细胞为节点,连接节点形成边,基于细胞间距离为边分配权值,构建加权联合图。接着,利用社区发现算法对图进行聚类,将来源不同但密切连接的细胞聚类在一起,实现整合,结果通常以图的形式输出。有的方法同时运用了配准和图聚类两种策略,在构建联合图后,还对数据集进行了配准。当输出结果为一个校正的基因表达矩阵时,可以进行广泛的下游分析;校正低维嵌入和联合图可用于可视化、聚类和拟时序分析,由于这两种结果没有改变原始表达值,因此不能直接用于差异表达分析、基因调控网络等基因层面的分析。可以使用批次感知的线性混合效应模型进行差异表达分析,或者对未校正的基因表达矩阵进行单独分析和共同比较[5-6]。此外,除了上述配准类和图聚类的整合方法,还有基于深度学习的整合方法。本文总结了近年来主要的多源scRNA-seq数据整合方法,包括整合策略、应用的主要算法、输出结果、程序语言和参考文献,如表1所示。
配准类方法中,多集典型相关分析(multi-set canonical correlation analysis,MultiCCA)是最早被报道的经典方法之一,支持在单细胞基因组学R工具包Seurat 2.0(Satija实验室,美国)中运行使用[7]。MultiCCA首先利用典型相关分析(canonical correlation analysis,CCA)将两个数据集映射到一个低维空间,接着利用动态时间规整(dynamic time warping,DTW)将两个数据集配准到一个共同的基准尺度,输出对齐的低维嵌入。同期,Haghverdi等[4]提出了以相互最近邻(mutual nearest neighbors,MNN)为核心的MNN校正(MNN-Correct),对于基准数据集的每个细胞,确定在目标数据集中k个最近邻(k-nearest neighbors,KNN),对于目标数据集也执行相同的操作;当一对细胞包含在彼此的k个近邻中,则视为一对MNN,代表不同数据集中具有相同细胞类型或状态的细胞对。MNN-Correct在余弦归一化的表达数据上计算细胞间的欧氏距离,得到一组MNN,MNN对的距离差值可以估算批次效应的大小,根据这些MNN对差异向量的加权平均值计算细胞特异的校正向量,将目标数据集向基准数据集配准。由于MNN-Correct在高维空间中运算,计算量巨大,快速MNN(fastMNN)作为MNN-Correct的新版本,增加了主成分分析(principal component analysis,PCA)降维步骤,运行速度提升。
受MNN的启发,在修拉(Seurat)版本三(Seurat V3)、单细胞合并(scMerge)、全景(Scanorama)中也引入了MNN的思想[6, 8-9]。Seurat V3在CCA定义的低维空间中识别MNN,并取名为锚点。为了确认锚点的可靠性,Seurat V3对锚点进行了过滤和评分:首先,返回高维的表达空间查看锚点的KNN条件是否同样满足,排除表达谱差异很大的错误锚点;然后,根据每对锚点细胞的邻域重叠率为每个锚点打分,在后续计算中降低低分数锚点的权重,以提高结果的鲁棒性和准确率。scMerge利用伽马-高斯混合模型识别一组稳定表达基因(stably expressed genes,SEGs);在PCA空间下识别不同数据集的相互最近集群,其中一定比例的最接近集群质心的细胞集称为伪重复,即不同数据集的同一类细胞;利用三步删除不必要变异模型(remove unwanted variation, 3-step,RUVⅢ)计算SEGs在伪重复中的表达差值,删除来源影响并配准数据集[10]。Scanorama将MNN的识别从两个数据集推广至多个数据集中,对于每个数据集的细胞,KNN的搜寻范围扩展为其余数据集的所有细胞,从而获取所有数据集中转录谱相似的MNN匹配,通过这些匹配将数据集拼接起来形成一张“全景图”。类似于MNN-Correct的校正计算,每个细胞的校正向量由其近邻范围内MNN对差异向量的加权平均得到,然后根据MNN细胞的占比确定每个数据集的合并顺序,依次合并到全景图中。Scanorama使用随机奇异值分解(randomized singular value decomposition,rSVD)进行降维,基于局部敏感哈希算法搜索近似的KNN,这大大加快了运行速度,可灵活应用于大数据集。鸿蒙(Harmony)是一种迭代聚类的整合方法,将所有数据集的细胞向聚类质心进行配准。在PCA空间中,Harmony使用k均值软聚类将相似细胞聚类在一起,并且最大化每个集群内的数据集多样性;在每个集群中,根据集群质心和数据集特异质心,计算细胞特异的校正因子来校正每个细胞;迭代运行以上两个过程,相当于迭代消除数据集来源对细胞聚类分配的影响,直到细胞的低维嵌入收敛到精确值,即完成整合,输出校正的低维嵌入[5]。
批次平衡k近邻(batch balanced k nearest neighbours,BBKNN)、样本网络聚类法(clustering on network of samples,Conos)和基因组实验关系关联推断(linked inference of genomic experimental relationships,LIGER)是图聚类的整合方法[11-13]。BBKNN首先在PCA空间中识别每个细胞在数据集内和不同数据集间的KNN并连接起来,然后按照一致流形近似和投影(uniform manifold approximation and projection,UMAP)的方法将近邻距离转换为连接分数,为每个连接赋予相应的权重,输出加权联合图用于聚类。Conos在低维空间中应用旋转空间和近邻(KNN和MNN)的策略,将数据集两两比较,建立数据集间和数据集内的边,用于构造联合图,其中根据细胞间相似性确定边的权重,并且降低数据集内的边相对于数据集间的边在图中的贡献;然后使用社区发现算法来获得聚类。LIGER应用综合非负矩阵分解(integrative non-negative matrix factorization,iNMF)进行降维,识别数据集特异和共享的两组因子来定义每个细胞的低维嵌入[14]。在因子空间下,计算细胞在数据集内和数据集间的KNN,并连接成对的相似细胞,即具有相似因子载荷模式的细胞,从而构建一个共享因子邻域图,用于社区发现算法聚类。聚类后,LIGER还对每个联合集群的因子载荷进行分位数归一化(quantiles normalization,QN),输出校正的低维嵌入。
近年来,研究人员将深度学习技术引入生物学领域,帮助处理生物数据资料,解决生物信息学问题,在数据整合方面也不断有新的进展。深度学习类方法的核心思想是利用深度网络学习数据结构和配准关系,处理优化问题,以克服批次效应。最大均值差异残差网络(maximum mean discrepancy and residual nets,MMD-ResNet)应用残差神经网络模型,以最小化基准数据集和目标数据集分布的最大均值差为目标,学习两个数据集之间的映射关系,使得两者在分布上相似[15]。深度MNN(deepMNN)也应用了残差神经网络,该方法利用PCA空间中识别的MNN对训练网络,并以最小化MNN对中的细胞间距、使网络的输出与输入相似为目标,指导网络学习从而消除数据集间的批次效应[16]。单细胞基因表达数据模型(scGen)结合了变分自编码器(variational auto-encoder,VAE)模型和潜在空间向量算法,只适用于已知细胞类型的数据整合。该方法的过程是输入标记细胞类型的数据,经过编码器映射到一个潜在空间中,在此空间下计算不同数据集中相同细胞类型的细胞间的差异向量,然后向其中一个数据集进行配准,再经过解码器映射回原本的高维表达空间,输出一个校正的基因表达矩阵[17]。Wang等[18]利用自编码器和生成式对抗网络(generative adversarial network,GAN)构建了对抗式配对风格迁移网络集成多源单细胞数据集(integration of multiple single-cell datasets by adversarial paired-style transfer networks,iMAP)模型,该模型分两个步骤进行数据整合,首先以保留生物变异、去除批次影响为优化目标,在自编码器中重建细胞表达谱,第二步在识别的MNN上训练GAN以正确匹配共享细胞类型的分布。其中,iMAP还引入了随机游走的策略来扩展MNN列表,从而更好地覆盖共享细胞类型的全部分布,有助于GAN的训练。
1.2 多源scRNA-seq数据整合的应用
许多整合方法现已应用于scRNA-seq数据分析中,去除批次效应,鉴定新的细胞类型,绘制单细胞图谱,并在肿瘤、脑科学、病毒等研究领域取得了一定成果。例如,Zhang等[19]利用Harmony整合了一位肝癌患者多个组织的、两种测序技术获得的免疫细胞,综合不同测序技术的优势能够绘制出更高分辨率的肝癌免疫图谱,通过联合分析还识别出了罕见的细胞群,对肝癌的免疫治疗具有指导意义。Liao等[20]利用Harmony和fastMNN整合来自三个不同捐赠者的肾脏组织单细胞转录组数据,鉴定出三个近端小管细胞的亚型和两个集合管细胞的亚型,为肾小管细胞的精确分类和相关疾病研究提供了重要的参考。Trujillo等[21]建立了一个皮层类器官模型用于模拟人脑早期发育,用Seurat V3整合四个培养时间点下所测序的类器官单细胞数据,经细胞注释,类器官发育过程中的细胞类型变化证明该模型具有功能性。这些整合方法在新型冠状病毒(severe acute respiratory syndrome coronavirus 2,SARS-CoV-2)的研究中也提供了有效助力。Qi等[22]收集了13种人体组织的单细胞基因表达数据,通过Harmony整合数据,确定了三个具有与SARS-CoV-2受体血管紧张素转化酶2最相似表达模式的候选基因,其可能是协助SARS-CoV-2入侵人体的共受体编码基因,有助于制定相应的干预策略。Zhang等[23]分别用Seurat V3和LIGER整合了胃和回肠数据集、结肠数据集,综合分析SARS-CoV-2受体血管紧张素转化酶2和跨膜丝氨酸蛋白酶2的共表达模式,揭示了SARS-CoV-2在消化系统中的潜在传播途径。Seurat V3还被用于整合来自健康者和不同严重程度的新冠肺炎患者的细胞数据,绘制出新冠肺炎患者的肺泡灌洗液免疫细胞图谱和免疫反应图谱[24-25]。除此以外,Tran等[26]评估了14种方法在多批次、大数据集、不同测序技术、细胞类型差异大等数据情形下的整合效果,Harmony、LIGER和Seurat V3是总得分前三的方法。
2 单细胞多模态数据整合
2.1 单细胞多模态数据的整合方法
2019年,《自然-方法》杂志将单细胞多模态组学选为“2019年度技术”,揭示了单细胞领域联合多个模态和组学进行测量分析的重要发展趋势[27-28]。单细胞多模态技术在一次实验中测量同一个细胞多个模态数据,表征细胞内不同层面的生物信息:DNA和RNA测序从底层获取基因序列和表达信息,检测基因动态变化,蛋白组和表观组从表层反映基因表达调控机制和分子性状的变化规律。多模态研究为后续的整合分析提供了丰富的数据,然而复杂的数据类型和数据特征也带来了计算上的挑战。例如,转录组数据通常是基因表达矩阵;DNA甲基化提供甲基化信号值矩阵,包含甲基化位点和甲基化水平信息;蛋白组提供蛋白表达矩阵;染色质可及性提供了开放区域的信息,称为峰矩阵。对于这些描述同一细胞不同模态的数据,其整合的目标是揭示不同模态之间的关联性,并且更详细地描述细胞状态和基因调控机制[29]。因此,一个关键问题就是如何将不同模态的信息联系到一起。第一种策略是基于已有假设,将其他模态信息关联到基因层面,进行基因级矩阵转换[28]。例如,对于染色质可及性,可以将峰矩阵转成基因活跃度矩阵,矩阵的行为基因、列为细胞,数值为基因活跃度得分,由跨越基因体和启动子区域(通常上游2~3 kb)的峰读段数相加得到[30];对于DNA甲基化数据,利用基因体的非CG甲基化推断基因表达,生成基因级甲基化数据矩阵[31]。经过这样的转换,其他模态的基因级矩阵与scRNA-seq矩阵的整合等同于两个不同来源scRNA-seq数据的整合,即本文第1.1节讨论的情况。第二种策略是直接从模态数据的本质、特征、关系出发,输出结果是模态信息的综合体,用来揭示模态间的因果关系。根据不同方法模型的处理,结果以不同形式呈现。
传统转录组测序领域早已有许多成熟有效的多组学整合算法和工具,虽然整合原理同样适用于单细胞多模态数据,但还需要适当的调整和实验验证[32-33]。本文总结了适用于单细胞多模态数据整合的主要方法,包括整合的数据类型、主要算法、程序语言和参考文献,如表2所示。
Seurat V3、Conos和LIGER应用了基因级矩阵转换的策略,将单细胞染色质开放区转座酶可及性测序(single-cell assay for transposase-accessible chromatin using sequencing,scATAC-seq)的峰矩阵转换为基因活跃度矩阵,与scRNA-seq矩阵相整合,LIGER还实现了基因级甲基化数据与scRNA-seq数据的整合,原理如本文第1.1节所述[6, 12-13]。
聚类通常是单细胞数据下游分析的第一步,对于scRNA-seq和scATAC-seq都有相应的聚类方法识别细胞类型,而协同两者信息能更好地解释细胞类型。Duren等[34]从聚类的角度,提出耦合scRNA-seq和scATAC-seq两个聚类过程的耦合非负矩阵分解(coupled nonnegative matrix factorization,Coupled NMF)模型,首先向模型中输入两个数据矩阵,然后求解耦合聚类的优化问题,输出两个数据的共同聚类结果以及对应的峰与基因配对。这一模型在反卷积与耦合聚类(de-convolution and coupled-clustering,DC3)模型中得到了升级,加入了群体细胞数据的反卷积化步骤,以获得细胞亚群特异的数据,帮助改善单细胞耦合聚类的结果[35]。多组学因子分析(multi-omics factor analysis,MOFA)是一种因子分析模型,适用于单细胞DNA甲基化与RNA数据的整合分析[36]。向MOFA模型中输入不同模态的数据矩阵,经过矩阵分解,结果输出一组因子,这些因子代表了驱动不同模态数据异质性的因素;下游分析中,可以得到影响因子贡献度的重要基因和甲基化信息,并对所有样本进行可视化、聚类、富集分析。流行对齐表征实验关系(manifold alignment to characterize experimental relationships,MATCHER)应用了流形学习的降维策略,通过高斯过程潜变量模型(Gaussian process latent variable model,GPLVM)将不同模态数据(基因表达、DNA甲基化、染色质可及性)映射到一维流形空间,每种模态数据都可以由一组伪时间值表示,并对伪时间值进行QN以实现统一度量,通过比较和分析该值可以研究多个模态之间的相关性和潜在调控机制[37]。
2.2 单细胞多模态数据整合的应用
有效的计算方法为充分挖掘不同模态的单细胞数据提供了机会,本节讨论单细胞多模态数据整合方法的重要应用,包括揭示细胞异质性、发现模态间交互关系和推断基因调控网络等。Welch等[37]利用MATCHER研究小鼠胚胎干细胞和人类诱导多能干细胞的转录组和表观基因组之间的相关性,结果揭示单细胞基因表达和DNA甲基化、染色质可及性、组蛋白修饰之间具有共同的变异模式,轨迹分析显示了细胞从多能性到分化启动状态的变化。Argelaguet等[36]也分析了小鼠胚胎干细胞的转录组和甲基化数据,MOFA分析结果同样揭示细胞分化过程中转录组和甲基化水平具有协同变化。Argelaguet等[38]应用三重组学测序技术获取小鼠胚胎细胞的单细胞核小体、甲基化和转录组三重模态数据,通过MOFA整合分析揭示了原肠胚形成过程中具有谱系特异的表观遗传模式和标志基因,并在三个模态推断因子中观察到了细胞间的异质性。因此,MOFA从最初的转录组与甲基化二重模态分析拓展到了多重模态分析,应用范围扩大。Lake等[39]通过训练一个梯度提升回归模型,在成人大脑细胞的转录组和表观基因组之间建立映射,发现了驱动细胞异质性的调节元件和转录因子,为研究大脑的复杂过程提供了新的思路。LIGER被用于整合scRNA-seq和DNA甲基化数据,联合定义了小鼠皮质细胞类型,揭示了细胞类型特异性的表观基因组调控机制[13]。研究单个细胞不同模态之间的相关性可帮助研究人员深层次理解基因调控网络,提高细胞类型分类的准确性和可解释性,使全面探索细胞身份和行为成为可能。
3 总结与展望
单细胞数据的整合分析提供了更为全面且深入的见解,帮助从不同层面剖析细胞类型和状态,深入挖掘调控机制。近年来,单细胞数据整合方法的开发与应用方面已经有了显著成果。
多源scRNA-seq数据整合方法中,Seurat系列和MNN-Correct系列的方法假设数据集之间的差异完全源于技术性变化,适用于处理细胞类型相似的数据集。在处理细胞组成差异大的情况下,LIGER相对适用,因为iNMF可以保留数据集间差异,也可以识别相似之处。对于多数方法,基准数据和目标数据的选择会影响整合效果;当整合两个以上的数据集时需要迭代整合过程,整合次序也将影响结果。Scanorama应用全景拼接的原理,避免了输入数据集的顺序对整合效果的影响。Harmony对所有数据执行迭代聚类和校正,对数据集次序不敏感。Harmony和BBKNN运行速度快,均适用于处理大数据集。深度学习类方法在复杂数据中具有可扩展性和高性能,而对于量小的数据集,此类方法的性能较差,不利于网络训练。有些深度学习类方法只能实现有监督或半监督的整合,例如scGen需要输入有细胞类型标签的数据集。deepMNN通过最小化MNN对的细胞间距,促进网络同时消除多个数据集的批次效应,因此能够实现一步整合。此外,深度学习类方法使用图形处理器来加速计算,同样适用于处理大数据集。
单细胞多模态数据整合方法中,Seurat V3、LIGER和MATCHER能够实现基因表达与DNA甲基化、染色质开放程度、组蛋白修饰等多个模态的整合分析,完整地建立起单细胞转录组和表观基因组之间的调控关系。DC3还实现了单细胞与群体细胞数据的联合分析,有效改善细胞的聚类结果。MOFA是相对流行的方法,从因子分析的角度解析数据异质性,使得分析结果也更具解释性。
尽管目前的单细胞数据整合方法具有良好的应用价值,但仍存在以下几个方面的挑战。首先,合适的整合方法是有效整合分析的基础,其选择依赖于对不同整合方法的性能测试,所以形成成熟的评测基准或评测系统至关重要,能够指导数据整合工作,减少方法本身对整合结果的影响。第二个挑战在于整合方法通常基于不同的计算平台开发,方法的推广和使用受限于用户对环境平台的偏好、数据格式和预处理步骤等因素,因此未来需要向跨平台分析、数据转换共享、流程化运行的趋势发展。第三个方面,单细胞数据规模日益扩大,以及适应“人类细胞图谱计划”(Human Cell Atlas,HCA)的需求,如何高效处理大规模数据集是亟待解决的问题,而大数据集恰好与深度学习的数据驱动性质相适应,这可能是未来可以继续推进研究的方向。随着新整合方法不断突破和生物信息学的发展,单细胞数据的整合分析定会提供更全面的细胞视角,帮助解决生物学难题。
利益冲突声明:本文全体作者均声明不存在利益冲突。