1.一种推断空间转录组内细胞亚群表达模式的方法，其特征在于，包括：
对scRNA-seq数据集进行质量控制和预处理，得到细胞亚群表达矩阵；
构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布；
在训练好的隐变量分布中进行采样，生成细胞亚群的表达模式；
基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有空间域的表达模式进行解卷积，得到细胞亚群在空间域中分布的最大后验估计，具体的：采用10X Visium测序技术得到的FFPE_Kidney空间转录组数据Y，通过空间聚类方法，为组织切片划分区域，将和Y作为解卷积方法的输入，输出每个区域内中每个细胞亚群所占的比例；
所述对scRNA-seq数据集进行质量控制和预处理，包括：过滤基因含量过低的细胞以及未在细胞中表达的基因和线粒体基因，并筛选出高表达的基因；
构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布方法如下：
对于一个预处理后的单细胞转录组基因表达矩阵，首先输入一层由全连接层组成的编码器，输出为和；再从高斯分布中进行采样得到隐变量Z，最后通过一层由全连接层组成的解码器生成最终的参考数据；
变分神经网络的公式表示如下：

其中，和分别表示编码器和解码器的隐藏层；和表示隐空间高斯分布的参数；表示隐变量；表示细胞亚群重建后的表达矩阵；
设定激活函数、损失函数和重参数化方法；所述损失函数表达式为：

其中，用于表示和的占比。
2.根据权利要求1所述的推断空间转录组内细胞亚群表达模式的方法，其特征在于，对隐变量z进行重参数化的表达式为：

其中，。
3.一种推断空间转录组内细胞亚群表达模式的系统，其特征在于，包括：
质量控制和预处理模块，被配置为：对scRNA-seq数据集进行质量控制和预处理，得到细胞亚群表达矩阵；
隐变量分布学习模块，被配置为：构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布；
表达模式生成模块，被配置为：在训练好的隐变量分布中进行采样，生成细胞亚群的表达模式；
解卷积模块，被配置为：基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有空间域的表达模式进行解卷积，得到细胞亚群在空间域中分布的最大后验估计，具体的：采用10X Visium测序技术得到的FFPE_Kidney空间转录组数据Y，通过空间聚类方法，为组织切片划分区域，将和Y作为解卷积方法的输入，输出每个区域内中每个细胞亚群所占的比例；
所述对scRNA-seq数据集进行质量控制和预处理，包括：过滤基因含量过低的细胞以及未在细胞中表达的基因和线粒体基因，并筛选出高表达的基因；
构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布方法如下：
对于一个预处理后的单细胞转录组基因表达矩阵，首先输入一层由全连接层组成的编码器，输出为和；再从高斯分布中进行采样得到隐变量Z，最后通过一层由全连接层组成的解码器生成最终的参考数据；
变分神经网络的公式表示如下：

其中，和分别表示编码器和解码器的隐藏层；和表示隐空间高斯分布的参数；表示隐变量；表示细胞亚群重建后的表达矩阵；
设定激活函数、损失函数和重参数化方法；所述损失函数表达式为：

其中，用于表示和的占比。
4.一种计算装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现上述权利要求1-2任一所述的方法的步骤。
5.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时执行上述权利要求1-2任一所述的方法的步骤。

技术领域
[0001]本发明属于生物信息学空间转录组测序数据分析技术领域，尤其是涉及一种推断空间转录组内细胞亚群表达模式的方法及系统。
背景技术
[0002]本部分的陈述仅仅是提供了与本公开相关的背景技术信息，不必然构成在先技术。
[0003]空间转录组学是一门生命科学与计算机科学的交叉学科。该领域的突破性进展已经为疾病和生物过程的研究带来了新的发现。但是，由于目前测序技术的局限性：空间转录组学技术能够测量转录物产生的位置，却无法分别哪个细胞产生了转录物。而单细胞技术(scRNA-seq)可以获得每个细胞的转录物，尽管会丢失空间信息。
[0004]一些分析工具整合了单细胞数据和空间转录组数据，并提出了解卷积的方法，即将每一个采样点(spot或bead)看成多种细胞类型的混合。它以单细胞中细胞亚群的表达模式为依据构建模型，以空间转录组每个spot的实验数据作为输入，产生的输出为在给定spot的基因表达分布下，细胞亚群在空间中分布的最大后验估计。
[0005]发明人发现，目前解卷积法对于细胞亚群表达模式的要求非常高，原始的scRNA-seq数据规模大，噪声多，这会导致解卷积法的运行速度慢、效果一般。而直接在数据中进行下采样则会丢失大量有价值的信息。
[0006]因此，开发一种获得细胞亚群表达模式的方法以解决以上问题是非常有必要的。
发明内容
[0007]本发明的目的在于提供一种推断空间转录组内细胞亚群表达模式的方法及系统，使得空间转录组中解卷积法所需的单细胞参考数据在降低维度的同时，保留大量有关信息，从而提高解卷积方法的运行速度和准确性。
[0008]为实现上述目的，本发明的一个或多个实施例提供了如下技术方案：
[0009]第一方面，本发明的一种推断空间转录组内细胞亚群表达模式的方法，包括：
[0010]对scRNA-seq数据集进行质量控制和预处理，得到细胞亚群表达矩阵；
[0011]对细胞亚群表达矩阵进行标准化和归一化；
[0012]构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布；
[0013]在训练好的隐变量分布中进行采样，生成细胞亚群的表达模式；
[0014]基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有空间域的表达模式进行解卷积，得到细胞亚群在空间域中分布的最大后验估计。
[0015]优选的，所述的对scRNA-seq数据集进行质量控制和预处理，包括：过滤基因含量过低的细胞以及未在细胞中表达的基因和线粒体基因，并筛选出高表达的基因。
[0016]优选的，对细胞亚群表达矩阵的标准化和归一化的方法如下：
[0017]Xi＝log(Xi+1)，i∈C
[0018]
[0019]其中Xi表示每个细胞亚群的表达矩阵，标准化采用的是log标准化方法，归一化采用的是min-max归一化方法；得到的表达矩阵X′i的取值范围为[0，1]。
[0020]优选的，构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布方法如下：
[0021]对于一个预处理后的单细胞转录组基因表达矩阵Xi，首先输入一层由全连接层组成的编码器，输出为μ和σ；再从高斯分布Norm(μ，σ2)中进行采样得到隐变量Z，最后通过一层由全连接层组成的解码器生成最终的参考数据；
[0022]神经网络的公式表示如下：
[0023]E＝ReLU(XiWE)
[0024]μ＝ReLU(XiWμ)
[0025]σ＝ReLU(XiWσ)
[0026]Z＝Sample[Norm(μ，σ2)]
[0027]D＝ReLU(ZWD)
[0028]
[0029]其中，E和D分别表示编码器和解码器的隐藏层；μ和σ表示隐空间高斯分布的参数；Z表示隐变量；X′i表示细胞亚群i重建后的表达矩阵。
[0030]优选的，在上述方法中还包括步骤：设定激活函数、损失函数和重参数化方法。
[0031]优选的，损失函数表达式为：
[0032]
[0033]其中，α用于表示||Xi-X′i||2和的占比。
[0034]优选的，对隐变量z进行重参数化的表达式为：
[0035]Z＝Sample[Norm(μ，σ2)]＝μ+εσ
[0036]其中，ε～Norm(0，1)。
[0037]第二方面，本发明的一种推断空间转录组内细胞亚群表达模式的系统，包括：
[0038]质量控制和预处理模块，被配置为：对scRNA-seq数据集进行质量控制和预处理，得到细胞亚群表达矩阵；
[0039]标准化和归一化模块，被配置为：对细胞亚群表达矩阵进行标准化和归一化；
[0040]隐变量分布学习模块，被配置为：构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布；
[0041]表达模式生成模块，被配置为：在训练好的隐变量分布中进行采样，生成细胞亚群的表达模式；
[0042]解卷积模块，被配置为：基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有空间域的表达模式进行解卷积，得到细胞亚群在空间域中分布的最大后验估计。
[0043]以上一个或多个技术方案存在以下有益效果：
[0044]本发明利用变分自编码器能够精确的获取scRNA-seq数据集中各个细胞亚群的表达模式，使得空间转录组中解卷积方法能够在给定spot的基因表达分布下，精确地得到细胞亚群在空间中分布的最大后验估计。
[0045]本发明使得空间转录组中解卷积法所需的单细胞参考数据在降低维度的同时，保留大量有关信息，从而提高解卷积方法的运行速度和准确性，使细胞在组织切片内的分布更加精确。
[0046]当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有优点。
附图说明
[0047]为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0048]图1为本发明变分自编码器结构图；
具体实施方式
[0049]下面结合附图与实施例对本公开做进一步说明。
[0050]应该指出，以下详细说明都是例示性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。
[0051]需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
[0052]此后参照附图描述本发明的具体实施例；然而，应当理解，所发明的实施例仅仅是本发明的实例，其可采用多种方式实施。熟知和/或重复的功能和结构并未详细描述以避免不必要或多余的细节使得本发明模糊不清。因此，本文所发明的具体的结构性和功能性细节并非意在限定，而是仅仅作为权利要求的基础和代表性基础用于教导本领域技术人员以实质上任意合适的详细结构多样地使用本发明。
[0053]在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。
[0054]实施例一
[0055]本发明实施例一提供了一种推断空间转录组内细胞亚群表达模式的方法，该方法可以应用在空间转录组学、单细胞转录组学等领域中，能够结合变分自编码器从而精确地获取细胞亚群的表达模式，进而利用解卷积方法，给出细胞亚群在空间中分布的最大后验估计，该方法包括以下步骤：
[0056]步骤1：对scRNA-seq数据集进行质量控制，本实施例中选用了Tabula-muris数据集中年龄为18个月的小鼠的肾脏细胞数据，是一个由3138个细胞和20138个基因组成的表达矩阵，记为X。对其进行质量控制，过滤基因含量过低的细胞以及未在细胞中表达的基因，以及筛选出具有高表达的基因。经过预处理，表达矩阵X分别由2771个细胞和3000个高变基因组成。
[0057]步骤2：对细胞亚群表达矩阵进行标准化和归一化。Tabula-muris给出了每个细胞所属的细胞亚群C，本实施例中剔除了细胞数小于25的细胞亚群，并对细胞亚群表达矩阵Xi(i∈C)进行log标准化和min-max归一化，步骤公式所示：
[0058]Xi＝log(Xi+1)，i∈C
[0059]
[0060]步骤3：构建变分自编码器(VAE)，以学习scRNA-seq数据集中细胞亚群表达矩阵的隐变量分布。本实施例中，变分自编码器属于神经网络的一种，它通过节点之间的连接来实现细胞表达模式的学习，并以高斯分布的方式描述对隐变量的观察，最后通过隐变量重建细胞亚群表达模式。在本实施例中，对于一个单细胞转录组的基因表达矩阵Xi，首先经过一层由全连接层组成的编码器(Encoder)，其输出为μ和σ，再从高斯分布Norm(μ，σ2)中进行采样得到隐变量Z，最后通过一层由全连接层组成的解码器(Decoder)生成最终的参考数据。
[0061]E＝ReLU(XiWE)
[0062]μ＝ReLU(XiWμ)
[0063]σ＝ReLU(XiWσ)
[0064]Z＝Sample[Norm(μ，σ2)]
[0065]D＝ReLU(ZWD)
[0066]
[0067]其中，E和D分别表示编码器和解码器的隐藏层，在本实施例中其维度为400；WE和WD分别表示全连接层的权重参数；μ和σ表示隐空间高斯分布的参数；Z表示隐变量，在本实施例中其维度为20；X′i表示细胞亚群i重建后的表达矩阵。
[0068]进一步的，由于标准化后的表达矩阵输入值均在0到1之间，故隐藏层采用ReLU激活函数，输出层采用sigmoid函数。VAE的损失函数可表示为：
[0069]Loss＝Ez～q(z|x)[logp(x|z)]+KL(N(μ，σ2)||N(O，I))其中第一项也称为reconstruction loss，即重建损失，本文的模型采用L2 loss，即：
[0070]||Xi-X′i||2
[0071]第二项为KL loss，用于反映重建后的表达模式和原始细胞亚群表达模式的拟合程度，在VAE中可表示为：
[0072]
[0073]因此最终的损失函数表示为：
[0074]
[0075]其中α用于表示重建损失和KL loss的占比，在本实施例中设置为2。在反向传播时，由于采样操作是不可导的，因此我们需要对隐变量z进行重参数化(Reparameterization)。
[0076]由于Z～N(μ，σ2)，可令：
[0077]Z＝Sample[Norm(μ，σ2)]＝μ+εσ
[0078]其中ε～Norm(0，1)。通过这个技巧，梯度可以直接通过μ和σ反向传播。
[0079]步骤4：在训练好的细胞亚群隐变量分布中进行采样，生成细胞亚群的表达模式，具体为：对于每一个细胞数大于25的细胞亚群，作为变分自编码器的输入在本实施例中，设置最大的迭代次数为1000，学习率为10-3，当KL loss值小于10-5时，训练停止。对于输出的结果，进行降采样至维度为25，得到标准的参考细胞亚群。
[0080]步骤5：基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有spot的表达模式进行解卷积，得到细胞亚群在空间中分布的最大后验估计，具体为：在本实施例中，采用了10X Visium测序技术得到的FFPE_Kidney空间转录组数据Y，其在组织上的spot有3124个，具有19465个基因，与S4步骤得到的细胞亚群中，有交集的基因有2675个。通过空间聚类方法，为组织切片划分区域，将X′和Y作为解卷积方法的输入，输出的是每个区域内X′中每个细胞亚群所占的比例。
[0081]需要注意的是，空间聚类方法可以采用Seurat、BayesSpace和SpaGCN等方法，解卷积方法可以采用SPOTlight、spacexr以及stereoScope方法，这些方法都是公知的，均在本专利的保护范围内。
[0082]实施例二
[0083]本实施例的目的是提供一种计算装置,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述实施例一中的方法的步骤。
[0084]实施例三
[0085]本实施例的目的是提供一种计算机可读存储介质。
[0086]一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时执行上述实施例一中的方法的步骤。
[0087]实施例四
[0088]本实施例的目的是提供一种推断空间转录组内细胞亚群表达模式的系统，包括：
[0089]质量控制和预处理模块，被配置为：对scRNA-seq数据集进行质量控制和预处理；
[0090]标准化和归一化模块，被配置为：对细胞亚群表达矩阵进行标准化和归一化；
[0091]隐变量分布学习模块，被配置为：构建变分神经网络，以学习scRNA-seq数据集中每个细胞亚群表达矩阵的隐变量分布；
[0092]表达模式生成模块，被配置为：在训练好的隐变量分布中进行采样，生成细胞亚群的表达模式；
[0093]解卷积模块，被配置为：基于所述的细胞亚群的表达模式，对空间转录组组织切片中所有空间域的表达模式进行解卷积，得到细胞亚群在空间域中分布的最大后验估计。
[0094]以上实施例二、三和四的装置中涉及的各步骤与方法实施例一相对应，具体实施方式可参见实施例一的相关说明部分。术语“计算机可读存储介质”应该理解为包括一个或多个指令集的单个介质或多个介质；还应当被理解为包括任何介质，所述任何介质能够存储、编码或承载用于由处理器执行的指令集并使处理器执行本公开中的任一方法。
[0095]本领域技术人员应该明白，上述本公开的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本公开不限制于任何特定的硬件和软件的结合。
[0096]以上所述仅为本公开的优选实施例而已，并不用于限制本公开，对于本领域的技术人员来说，本公开可以有各种更改和变化。凡在本公开的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。
[0097]上述虽然结合附图对本公开的具体实施方式进行了描述，但并非对本公开保护范围的限制，所属领域技术人员应该明白，在本公开的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本公开的保护范围以内。