整合单细胞转录组与空间转录组数据的分析方法及系统

摘要

本发明提出了整合单细胞转录组与空间转录组数据的分析方法及系统，包括：获取空间转录组数据，对空间转录组数据进行预处理；通过聚类方法将空间转录组数据中的采样点划分为若干个类别，每一个类别表示组织切片图像中一个组织区域；根据去污染后的空间基因表达数据，和空间聚类步骤的结果，筛选出在空间中具有差异性表达的基因；根据输入的单细胞数据进行分析，基于筛选出在空间中具有差异性表达的基因得出每个细胞亚群的表达模式；根据单细胞数据集中的细胞亚群表达模式，将所有spot的基因表达作为输入，获得细胞类型在各个组织区域中的分布。

说明书

技术领域

本发明属于数据分析技术领域，尤其涉及整合单细胞转录组与空间转录组数据的分析方法及系统。

背景技术

本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。

近年来，随着空间测序技术的兴起，对于从RNA水平对基因表达的转录组学研究正在从一维的单细胞转录组学转向二维的空间转录组学。由于空间转录组测序技术自身分辨率的局限性，目前还不能直接通过空间转录组学得到细胞亚群在组织切片中的分布。因此，越来越多的研究通过整合单细胞转录组与空间转录组数据，以给出细胞在组织切片中的分布，进而研究细胞的空间异质性、不同区域的细胞富集和交互等，为更多的生物学发现提供可能。

发明人发现，随着空间转录组测序技术的不断发展和完善，下游的整合分析技术也不断涌现。但是，这些分析工具往往只针对特定的测序技术，或特定的分析任务，无法满足当下空间转录组的分析需求。

发明内容

为克服上述现有技术的不足，本发明提供了整合单细胞转录组与空间转录组数据的分析方法，通过提供空间转录组测序数据和单细胞测序数据，能够直接得到单细胞数据中各种细胞类型在空间中的分布，能够满足大部分的空间转录组和单细胞转录组的联合分析需求。

为实现上述目的，本发明的一个或多个实施例提供了如下技术方案：

第一方面，公开了整合单细胞转录组与空间转录组数据的分析方法，包括：

获取空间转录组数据，对空间转录组数据进行预处理；

预处理之后去除空间转录组数据中游离的mRNA和处理基因表达的缺失值，得到去污染后的空间表达矩阵；

利用去污染后的空间表达矩阵、采样点在组织切片中的位置信息和组织学影像，通过聚类方法将空间转录组数据中的采样点划分为若干个类别，每一个类别表示组织切片图像中一个组织区域；

根据去污染后的空间基因表达数据，和空间聚类步骤的结果，筛选出在空间中具有差异性表达的基因；

根据输入的单细胞数据进行分析，得出每个细胞亚群在单细胞高表达基因下的表达模式；

根据单细胞数据集中在单细胞高表达基因下的细胞亚群表达模式，将所有spot的基因表达作为输入，获得细胞类型在各个组织区域中的分布。

需要说明的是，单细胞高表达基因是基于单细胞数据集得出的，空间变异基因是基于空间数据集得出的。

作为进一步的技术方案，还包括可视化步骤：在组织切片图像上，标注出采样点的所在位置，然后用不同颜色标记每个采样点的类别，每一种类别表示一个组织区域，以表示组织区域分割结果。

作为进一步的技术方案，获取的空间转录组数据，包括：

每个采样点的基因表达组成的表达矩阵、不同分辨率下的组织切片图像、每个采样点在组织切片图像中的具体位置和原始高分辨率图像与低分辨率图像之间的比例因子；

其中，组织切片，指的是载有生物组织的染色芯片，其上的采样点呈点阵排列，一部分采样点位于组织上，一部分采样点位于组织外。

作为进一步的技术方案，对空间转录组数据进行预处理，具体包括：

筛选出位于组织上的采样点；

过滤在所有采样点中表达值过低的基因和线粒体基因；

其中，线粒体基因，指线粒体中携带的遗传信息，在基因信息中以“MT-”为前缀进行标记。

作为进一步的技术方案，去除空间转录组数据中游离的mRNA，具体包括：

读取组织外采样点的表达值进行建模，推断出背景噪声的分布，再将组织区域内的spot表达值分解为实际的表达值和背景噪声两部分，消去组织内的背景噪声，其中，背景噪声指组织外采样点的基因表达模式。

作为进一步的技术方案，处理空间转录组基因表达的缺失值，具体包括：

分别在每个采样点在组织切片图像中的具体位置信息层面，以及空间表达矩阵层面层面上分别寻找相似的采样点，然后进行KNN平滑，以反映出组织中的RNA丰度。

作为进一步的技术方案，通过聚类方法将空间转录组数据中的采样点划分为若干个组织区域，具体包括：

对去污染后的空间表达矩阵中的基因进行降维，之后对降维后的每一个维度采用多元t分布模型进行建模和聚类，最后进行参数的更新；

或者，整合空间位置信息和组织学影像信息，计算去污染后的空间表达矩阵中每个spot之间的距离，并采用图卷积神经网络，整合距离信息与基因表达量，最后基于图卷积网络的输出，采用无监督深度嵌入聚类分析方法将空间转录组数据中的采样点划分为若干个组织区域；

或者，通过对每个采样点对应染色图区域和邻域信息对基因的表达进行均一化，再进行无监督聚类。

作为进一步的技术方案，根据输入的单细胞数据进行分析，得出每个细胞亚群的表达模式，具体包括：

对单细胞转录组数据集进行质量控制；

过滤基因含量过低的细胞以及未在细胞中表达的基因；

筛选出具有高表达的基因；

利用数据集给出的细胞类型标注，对细胞亚群进行随机下采样，每个细胞亚群中提取出少量细胞作为该细胞亚群在高表达基因下的表达模式，最后得到单细胞参考表达矩阵。

第二方面，公开了整合单细胞转录组与空间转录组数据的分析系统，包括：

数据获取模块，其用于获取空间转录组数据，判断测序数据类型；

去污染模块，其用于去除背景中游离的mRNA和处理基因表达的缺失值，提高空间中每个spot中mRNA的纯度和丰度；

空间聚类模块，其利用spot的基因表达、spot在组织切片中的位置信息和组织学影像，通过聚类方法将组织区域划分为若干个类别；

空间差异基因识别模块，其根据基因表达的数据流形，和空间聚类的先验知识，筛选出在空间中具有差异性表达的基因；

单细胞分析模块，其根据用户输入的单细胞数据进行分析，得出每个细胞亚群的表达模式；

解卷积模块，其根据单细胞数据集中的细胞亚群表达模式，将所有spot的基因表达作为输入，获得细胞类型在各个组织区域中的分布；

可视化模块，其用于将组织区域中的空间差异基因和细胞亚群分布进行可视化。

以上一个或多个技术方案存在以下有益效果：

本发明整合现有的分析工具，兼容测序技术，提出了一种标准化的分析工具，满足了空间转录组领域的分析需求。

本发明增加了空间去污染，有利于提高空间下游分析方法的准确性；只需提供空间转录组数据和单细胞转录组数据，即可自动化地实现联合分析工作流(pipline)；在每一个单元中，提供了R语言和python中多种高性能的分析工具，用户可以通过配置文件的修改，在每个处理步骤中可以实现对工具的自由配置和对工具软件参数的自定义设置。

本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1为本发明实施例所述的联合分析工作流程图。

具体实施方式

应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。

在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

mRNA：信使RNA，即信使核糖核酸，是由DNA的一条链作为模板转录而来的、携带遗传信息能指导蛋白质合成的一类单链核糖核酸。

spot：一个mRNA的采样点，通常包含3至10个细胞的mRNA。

FFPE：用福尔马林和石蜡包埋技术制备的组织切片。

FSV：每个基因的变异中空间变异解释的变异比例(Fraction of varianceexplained by spatial variation，FSV)。

H&E染色切片：苏木精和伊红染色切片。

实施例一

本实施例公开了整合单细胞转录组与空间转录组数据的分析方法，包括：

S1，数据获取：获取空间转录组数据，并进行预处理。

上述空间转录组数据，指由10XVisium测序技术得到的空间转录组测序数据，具体包括：

每个采样点(Barcodes)的基因表达(Features)组成的表达矩阵(spMatrix)、不同分辨率下的组织切片图像(Images)、每个采样点在组织切片图像中的具体位置(coordinates)和原始高分辨率图像与低分辨率图像之间的比例因子(scaleFactors)。其中，组织切片，指的是载有生物组织的H&E染色芯片，其上的采样点呈点阵排列，一部分采样点位于组织上，一部分采样点位于组织外。

在本实施例子中，获取空间转录组数据，具体包括：

利用计算机指令，读取硬盘中的空间转录组数据表达矩阵、采样点信息(Barcodes)、基因信息(Features)、不同分辨率下的组织切片图像(Images)、每个采样点在组织切片图像中的具体位置(coordinates)和原始高分辨率图像与低分辨率图像之间的比例因子(scaleFactors)，用于后续的分析。

在本实施例子中，对空间转录组数据进行预处理，具体包括：

筛选出位于组织上的采样点；

过滤在所有采样点中表达值过低的基因和线粒体基因。线粒体基因，指线粒体中携带的遗传信息，在基因信息中以“MT-”为前缀进行标记。

本实施例中，选用10XVisium测序技术的FFPE小鼠肾脏组织切片作为空间转录组数据集，其二维空间位置和组织学影像如图所示。对样本的表达矩阵进行质量控制，得到在组织上的spot有3124个，具有19465种基因类型。

S2，去污染：去除空间转录组数据中游离的mRNA和处理空间转录组基因表达的缺失值，提高空间中每个spot中mRNA的纯度和丰度。

在本实施例子中，去除空间转录组数据中游离的mRNA，具体包括：

采用SpotClean方法，读取组织外采样点的表达值进行建模，推断出背景噪声的分布，再将组织区域内的spot表达值分解为实际的表达值和背景噪声两部分，从而消去组织内的背景噪声，更好地反映出组织切片中的RNA纯度。其中，背景噪声，指组织外采样点的基因表达模式。

在本实施例子中，处理空间转录组基因表达的缺失值，具体包括：

采用SPCS方法，分别在每个采样点在组织切片图像中的具体位置信息层面，以及空间表达矩阵层面层面上分别寻找相似的采样点，然后进行KNN平滑，以更好地反映出组织中的RNA丰度。

经过上述两步的处理后，得到一个去污染后的空间表达矩阵。

S3，空间聚类：利用去污染后的空间表达矩阵、采样点在组织切片中的位置信息和组织学影像，通过聚类方法将空间转录组数据中的采样点划分为若干个类别，每一个类别表示组织切片图像中一个组织区域。

若干个类别可以作为空间差异基因识别的先导步骤，具体地，可以在每个类别中筛选标记基因左右空间差异基因；也是探究细胞亚群在组织切片中分布的基础。

利用去污染后的空间表达矩阵、采样点在组织切片中的位置信息和组织学影像，通过聚类方法将空间转录组数据中的采样点划分为若干个组织区域，具体包括：

采用BayesSpace方法，首先对去污染后的空间表达矩阵中的基因进行降维，之后对降维后的每一个维度采用多元t分布模型进行建模和聚类，最后用Metropolis-Hastings算法进行参数的更新，其中该算法采用了Potts模型整合的二维空间信息作为其先验分布；

或者，采用SpaGCN方法，整合空间位置信息和组织学影像信息(H&E染色图)，计算了去污染后的空间表达矩阵中每个spot之间的距离，并采用图卷积神经网络，整合距离信息与基因表达量，最后基于图卷积网络的输出，采用一种无监督深度嵌入聚类分析方法将空间转录组数据中的采样点划分为若干个组织区域；

或者，采用stlearn方法，通过对每个采样点对应H&E染色图区域和邻域信息对基因的表达进行均一化，再进行无监督聚类。

本公开的空间聚类步骤提供以上三种方法，用户可以选择采用一种或多种方法进行空间聚类，以充分利用二维空间位置信息和组织学影像，得到组织区域分割结果，结合去污染步骤种的去噪和增强效果，达到了比单独使用上述方法和传统的工作流更好的性能。

S4，空间差异基因识别：根据去污染后的空间基因表达数据，和空间聚类步骤的结果，筛选出在空间中具有差异性表达的基因。具有空间异质性的基因是组织区域的重要标识，对空间中细胞间的交互研究具有重要意义。

需要说明的是，空间差异表达基因是组织区域的一个重要“标志”，是本方案的输出结果之一，在细胞间交互研究等下游分析中使用。由于空间差异基因与单细胞高表达基因的交集较小，本技术方案在S6步骤中使用了所有的空间基因。

根据去污染后的空间基因表达数据，筛选出在空间中具有差异性表达的基因，具体包括：

采用SpatialDE方法，利用高斯过程回归模型，计算每个基因的变异中空间变异解释的变异比例(Fraction of variance explained by spatial variation，FSV)，并根据FSV的值筛选出空间变异基因。

或者，采用SPARK方法，利用多个空间核的广义线性模型计算空间差异基因。

或者，采用SPARK-X方法，通过非参数模型，提高SPARK的运行速度和有效性。

根据去污染后的空间基因表达数据，和空间聚类步骤的结果，筛选出在空间中具有差异性表达的基因，还包括：

采用SpaGCN方法，利用空间聚类步骤得到的组织区域划分结果，计算每个组织区域中采样点的基因表达，与相邻区域中采样点的基因基因表达的差异，得到每个区域的差异基因，并将其作为空间差异基因。

本公开的空间差异基因识别步骤提供以上四种方法，用户可以选择采用一种或多种方法进行空间差异基因识别。在本实施例中，本公开利用空间聚类步骤提供的聚类信息，筛选出了188个空间差异基因，相比于其他工作流，这些差异基因具有更高的空间异质性。

S5，单细胞分析：根据用户输入的单细胞数据进行分析，得出每个细胞亚群的表达模式。具体包括：

对单细胞转录组数据集进行质量控制；

过滤基因含量过低的细胞以及未在细胞中表达的基因；

筛选出具有高表达的基因；

利用数据集给出的细胞类型标注，对细胞亚群进行随机下采样，每个细胞亚群中提取出少量细胞作为该细胞亚群在单细胞高表达基因下的表达模式，最后得到单细胞参考表达矩阵。由于目前空间转录组尚不能达到单细胞分辨率，因此一个采样点是由多种细胞类型混合而成。单细胞参考表达矩阵可以作为单细胞的表达模式，利用之后的解卷积步骤，可以得到细胞亚群在每个组织区域上的表达分布。

上述单细胞参考表达矩阵可以作为单细胞的表达模式，利用之后的解卷积步骤，可以得到细胞亚群在每个组织区域上的表达分布。

在本实施例中，选用了Tabula-muris数据集中年龄为18个月的小鼠的肾脏细胞数据，是一个由3138个细胞和20138个基因组成的表达矩阵。对其进行质量控制，过滤基因含量过低的细胞以及未在细胞中表达的基因，以及筛选出具有高表达的基因。经过预处理，表达矩阵分别由2771个细胞和3000个高变基因组成。最后利用数据集给出的细胞类型标注，对14个细胞亚群分别进行随机下采样，每个细胞亚群中提取出25个细胞作为该细胞亚群的表达模式，最后得到的参考表达矩阵分别由350个细胞和3000个高变基因组成。

S6，解卷积：根据单细胞数据集中在单细胞高表达基因上的细胞亚群表达模式，将所有spot的基因表达作为输入，根据需要可以选择一种或多种方法获得细胞类型在各个组织区域中的分布。解卷积法将每一个spot看成多种细胞类型的混合，其根据单细胞数据集中的细胞亚群表达模式，将所有spot的基因表达作为输入，获得细胞类型在各个组织区域中的分布。

具体包括：

采用SPOTlight方法，该方法将空间数据与单细胞参考数据进行集成，非负矩阵因子回归方法进行初始化，然后根据每种细胞类型的标记基因和空间基因表达数据，采用非负最小二乘法得到每个采样点的细胞类型比例。

或者，采用spacexr方法，该方法认为空间转录组数据符合泊松分布，在给定提供的单细胞数据后，读取每种细胞类型中所有基因的平均表达，然后采用逐步拟合的方法，找到最符合采样点基因表达的一种或两种细胞类型，然后用最大似然估计法拟合泊松分布的参数，进而推断出采样点中的细胞类型比例。

或者，采用stereoscope方法，该方法认为单细胞参考表达数据和空间转录组数据都符合负二项分布，在给定提供的单细胞数据后，通过查找最大似然估计(MLE)获得细胞类型分布特定参数的值，最后根据观测到的空间数据，利用细胞类型的先验分布，估计每个采样点的细胞类型比例。

本公开的解卷积步骤提供以上三种方法，用户可以选择采用一种或多种方法进行解卷积，得到每个采样点的细胞类型分布，结合空间聚类步骤中得到的组织区域，可以得到每个组织区域的细胞类型分布，对于探究细胞的空间异质性、细胞在微环境中的交互作用具有重要意义。

S7，可视化：通过计算机指令，将组织区域、组织区域中的空间差异基因表达和细胞亚群分布进行可视化，具体包括：

可视化空间聚类步骤的结果，具体为：组织区域分割结果。本公开通过计算机指令，在组织切片图像上，标注出采样点的所在位置，然后用不同颜色标记每个采样点的类别，每一种类别表示一个组织区域，以表示组织区域分割结果。

可视化空间差异基因识别步骤的结果具体为：空间变异基因在组织切片上的表达。本公开通过计算机指令，在组织切片图像上，标注出采样点的所在位置，然后用不同深浅程度的颜色标注出空间差异基因在采样点上的表达。

可视化解卷积步骤的结果，具体为：细胞在组织区域中的分布、细胞亚群—组织区域热图以及细胞类型在组织区域中的比例。本公开通过计算机指令，在组织切片图像上，标注出采样点的所在位置，然后用不同深浅程度的颜色标注出细胞类型在采样点中所占比例；计算每个细胞类型在每个组织区域的采样点上的平均表达，得到细胞类型分布矩阵，进一步地，得到细胞亚群—组织区域热图；再通过堆叠柱状图的方式，将细胞类型在组织区域中的比例进行可视化。

上述方法可以应用在空间转录组学、单细胞转录组学等领域中，能够结合空间转录组数据和单细胞转录组数据，得到细胞类型在组织区域中的分布。该结果对探究细胞的空间异质性、细胞在微环境中的交互作用具有重要意义，是进行更多空间转录组下游分析的基础。

实施例二

本实施例的目的是提供一种计算机装置,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述方法的步骤。

实施例三

本实施例的目的是提供一种计算机可读存储介质。

一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时执行上述方法的步骤。

实施例四

本实施例的目的是提供整合单细胞转录组与空间转录组数据的联合分析系统，包括：

数据获取模块，其用于获取空间转录组数据(SpatialTranscriptomicsData)，判断测序数据类型；

去污染模块(Decontamination)，其用于去除背景中游离的mRNA和处理基因表达的缺失值，提高空间中每个spot中mRNA的纯度和丰度。

空间聚类模块(Clustering)，其利用spot的基因表达、spot在组织切片中的位置信息和组织学影像，通过聚类方法将组织区域划分为若干个类别。

空间差异基因识别模块(FindingSpatialVariableGenes，SVGs)，其根据基因表达的数据流形，和空间聚类的先验知识，筛选出在空间中具有差异性表达的基因。

单细胞分析模块(SingleCellAnalysis)，其根据用户输入的单细胞数据进行分析，得出每个细胞亚群的表达模式。

解卷积模块(Deconvolution)，其根据单细胞数据集中的细胞亚群表达模式，将所有spot的基因表达作为输入，获得细胞类型在各个组织区域中的分布。

可视化模块(Visualization)，其用于将组织区域中的空间差异基因和细胞亚群分布进行可视化。

该系统相比于之前的系统增加了空间去污染模块，在每个模块中整合了多种分析方法，提高了空间下游分析方法的准确性；用户只需提供空间转录组数据和单细胞转录组数据，即可自动化地实现联合分析工作流(pipline)；用户可以通过配置文件的修改，在每个处理步骤中可以实现对工具的自由配置和对工具软件参数的自定义设置。

该系统增加了去污染模块，有利于提高空间下游分析方法的准确性；只需提供空间转录组数据和单细胞转录组数据，即可自动化地实现联合分析工作流(pipline)；在每一个单元中，提供了R语言和python中多种高性能的分析工具，用户可以通过配置文件的修改，在每个处理步骤中可以实现对工具的自由配置和对工具软件参数的自定义设置。

以上实施例二、三和四的装置中涉及的各步骤与方法实施例一相对应，具体实施方式可参见实施例一的相关说明部分。术语“计算机可读存储介质”应该理解为包括一个或多个指令集的单个介质或多个介质；还应当被理解为包括任何介质，所述任何介质能够存储、编码或承载用于由处理器执行的指令集并使处理器执行本发明中的任一方法。

本领域技术人员应该明白，上述本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。