基于变分自编码器的单细胞转录因子调控网络构建方法与流程

1.本发明涉及基因转录调控领域，具体涉及基于变分自编码器的单细胞转录因子调控网络构建方法。


背景技术：

2.转录因子调控网络是指由多种转录因子相互作用和调控，调节基因转录的复杂网络。不同的转录因子可以通过与dna结合，改变某个基因的转录速率或转录起始位点进而激活或抑制基因的转录，从而影响该基因的表达水平。转录因子也可以相互作用，形成调控模块或调控通路，实现对基因表达的复杂调控。构建并研究转录因子调控网络可以探究转录调控机制、揭示细胞发育和分化的规律、研究疾病的发生与发展、为药物的研发提供新思路。然而，融合多种单细胞(单一种类的细胞)数据存在诸如批次效应、数据的高维稀疏性、特征信息离散等困难。
3.近年来，生物信息学家们提出了一系列基于单细胞atac-seq和单细胞rna-seq数据的转录因子调控网络构建方法。这些方法的基本思想是利用数理统计、机器学习等手段从海量单细胞数据中提取出转录因子之间的连接关系并应用到下游分析任务中。现有技术提出了深度学习模型deepwalk来学习转录因子之间的语义关系，并维持了转录因子的调控关系在向量空间中和在原始图中的语义一致性。genie3是一种推断转录因子-基因调控网络的机器学习方法。genie3使用调控因子的表达值来预测每个基因的表达值来得到一个随机森林模型。通过对每个基因的预测误差进行排序，从而推断出每个基因最可能受到哪个调控因子的调控作用，进而推断出转录因子调控网络。grnboost是一种基于梯度提升树算法的调控网络预测方法，可以在不需要任何先验知识的情况下，直接从单细胞rna-seq数据中学习调控关系。通过迭代的训练一系列弱分类器来构建一个强分类器，并在每次迭代中使用残差作为新的目标变量来训练下一个弱分类器，在多次迭代后尝试将模型的误差最小化以获得更准确的调控网络。scenic是一种基于共表达模式的转录因子调控网络预测算法，首先通过共表达网络鉴定转录因子的下游目标基因，其次使用基因表达模型来推断每个转录因子的调控模式，用于预测细胞类型并构建调控网络。
4.现有技术存在以下问题：
5.(1)现有方法构建的转录因子调控网络，其中转录因子的不同测序技术存在不一致和数据集不匹配的情况，并采用单一的单细胞atac-seq数据或单细胞rna-seq数据中的高阶特征，难以有效地同时利用单细胞atac-seq数据和单细胞rna-seq数据中的高阶特征；
6.(2)现有方法忽略了转录因子调控网络内部的结构信息，进而无法解释各个转录因子在调控网络内的作用，导致无法全面准确地表示转录因子调控网络。


技术实现要素：

7.针对现有技术中的上述不足，本发明提供了基于变分自编码器的单细胞转录因子调控网络构建方法，能提取并对齐单细胞atac-seq和单细胞rna-seq数据中的高阶特征，结
构化转录因子的调控关系，进而更全面准确地构建转录因子调控网络。
8.为了达到上述发明目的，本发明采用的技术方案为：
9.基于变分自编码器的单细胞转录因子调控网络构建方法，包括以下步骤：
10.s1、根据scatac-seq数据和scrna-seq数据，确定输入矩阵和真实标签矩阵；
11.s2、构建转录因子调控网络模型，并根据步骤s1中的输入矩阵和转录因子调控网络模型确定预测标签矩阵；
12.s3、根据步骤s1中的真实标签矩阵和步骤s2中的预测标签矩阵，确定转录因子调控网络模型损失；
13.s4、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对转录因子调控网络模型的内部参数进行更新；
14.s5、采用步骤s4中更新内部参数后的转录因子调控网络模型构建转录因子调控网络。
15.进一步地，步骤s1包括以下分步骤：
16.s11、利用生信数据分析工具对scatac-seq数据和scrna-seq数据进行数据筛选；
17.s12、利用转录因子对齐分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，得到第一输入矩阵和第二输入矩阵；
18.s13、利用基序扫描工具扫描分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，确定真实标签矩阵。
19.进一步地，步骤s2包括以下分步骤：
20.s21、构建包括第一变分自编码器、第二变分自编码器和解码器的转录因子调控网络模型；
21.s22、训练第一变分自编码器，并输入分步骤s12中的第一输入矩阵至训练后的第一变分自编码器，得到第一特征信息，表示为：
[0022][0023]
其中：为特征信息，此处为第一特征信息，tanh为双曲正切激活函数，w3为第一变分自编码器的第三全连接层神经元的权重，w2为第一变分自编码器的第二全连接层神经元的权重，dropout为dropout层，该层以50％的概率随机失活全连接层中的神经元，w1为第一变分自编码器的第一全连接层神经元的权重，xi为输入矩阵，此处为第一输入矩阵，bias为第一变分自编码器的第一全连接层神经元的偏置系数；
[0024]
s23、训练第二变分自编码器，并输入分步骤s12中的第二输入矩阵至训练后的第二变分自编码器，得到第二特征信息，表示为：
[0025][0026]
其中：为特征信息，此处为第二特征信息，w6为第二变分自编码器的第三全连接层神经元的权重，w5为第二变分自编码器的第二全连接层神经元的权重，w3为第二变分自
编码器的第一全连接层神经元的权重，xi为输入矩阵，此处为第二输入矩阵，biasr为第二变分自编码器的第一全连接层神经元的偏置系数；
[0027]
s24、训练解码器，并通过训练后的解码器得到预测标签矩阵，表示为：
[0028][0029]
其中：y
ip
为预测标签矩阵，a为包含转录因子调控关系的邻接矩阵，w8为解码器的第二层全连接层神经元的权重，w7为解码器的第一层全连接层神经元的权重，xi为输入矩阵，此处包括第一输入矩阵和第二输入矩阵，为特征信息，此处包括第一特征信息和第二特征信息。
[0030]
进一步地，在分步骤s21中，第一变分自编码器和第二变分自编码器相互独立但结构相同，均包括一个有偏置系数的全连接层、一个dropout层和两个无偏置系数的全连接层；解码器包括双层全连接层模块，并通过点乘的方式嵌入包含转录因子调控关系的邻接矩阵。
[0031]
进一步地，在分步骤s24中，分别通过真实数据和伪数据依次对解码器进行训练，伪数据表示为：
[0032][0033]
其中：为伪数据，θ
d*
为解码器的内部更新参数，为特征信息，此处包括第一特征信息和第二特征信息。
[0034]
进一步地，步骤s3包括以下分步骤：
[0035]
s31、计算转录因子调控网络模型的交叉熵损失函数；
[0036]
s32、计算转录因子调控网络模型的kl散度；
[0037]
s33、计算步骤s1中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数；
[0038]
s34、根据分步骤s31中转录因子调控网络模型的交叉熵损失函数、分步骤s32中转录因子调控网络模型的kl散度和分步骤s33中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数，计算转录因子调控网络模型损失。
[0039]
进一步地，在分步骤s34中，计算转录因子调控网络模型损失，表示为：
[0040][0041]
其中：l为转录因子调控网络模型损失，为转录因子调控网络模型的交叉熵损失函数，xi为输入矩阵，y
ip
为预测标签矩阵，β为，kl(μr,σr)||(μa,σa)为转录因子调控网络模型的kl散度，μr为scrna-seq数据特征分布的均值，σr位scrna-seq数据特征分布的标准差，μa为scatac-seq数据特征分布的均值，σa为scatac-seq数据特征分布的标准差，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵。
[0042]
进一步地，步骤s4包括以下分步骤：
[0043]
s41、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对编码器
的内部参数进行更新；
[0044]
s42、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对解码器器的内部参数进行更新。
[0045]
进一步地，步骤s41包括以下分步骤：
[0046]
s411、计算编码器的内部更新参数，表示为：
[0047][0048]
其中：θ
e*
为编码器的内部更新参数，为当转录因子调控网络模型的交叉熵损失最小时所取得的编码器内部参数，v
l
为第一数据集，vi为第一数据集中的元素，为转录因子调控网络模型的交叉熵损失函数，为特征信息，y
ip
为预测标签矩阵；
[0049]
s412、根据步骤s3中的转录因子调控网络模型损失计算编码器的反向传播梯度；
[0050]
s413、判断分步骤s412中的反向传播梯度是否下降；若是则将编码器的内部参数更新为分步骤s411中的计算值，否则对编码器的内部参数进行更新并返回分步骤s411；
[0051]
对编码器的内部参数进行更新，表示为：
[0052][0053]
其中：θ
′e为更新后的编码器的内部参数，为编码器的学习速率，为编码器的反向传播梯度。
[0054]
进一步地，步骤s42包括以下分步骤：
[0055]
s421、计算解码器的内部更新参数，表示为：
[0056][0057]
其中：θ
d*
为解码器的内部更新参数，为当真实标签矩阵和预测标签矩阵的均方误差最小时所取得的解码器内部参数，v
p
为第二数据集，vi为数据集中的元素，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵，y
ip
为预测标签矩阵；
[0058]
s422、根据步骤s3中的转录因子调控网络模型损失计算解码器的反向传播梯度；
[0059]
s423、判断分步骤s422中的反向传播梯度是否下降；若是则将解码器的内部参数更新为分步骤s421中的计算值，否则对解码器的内部参数进行更新并返回分步骤s421；
[0060]
对解码器的内部参数进行更新，表示为：
[0061][0062]
其中：θ
′d为更新后的解码器的内部参数，为解码器的学习速率，为解码器的反向传播梯度。
[0063]
本发明具有以下有益效果：
[0064]
(1)本发明通过构建的第一变分自编码器、第二变分自编码器和解码器的转录因子调控网络模型，能对齐并提取单细胞atac-seq和单细胞rna-seq数据中的高阶特征；
[0065]
(2)本发明通过点乘的方式将包含转录因子调控关系的邻接矩阵嵌入解码器，能结构化转录因子的调控关系；
[0066]
(3)本发明采用双级优化策略对转录因子调控网络模型的内部参数进行更新，进而能更全面准确地构建转录因子调控网络。
附图说明
[0067]
图1为基于变分自编码器的单细胞转录因子调控网络构建方法流程示意图。
具体实施方式
[0068]
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
[0069]
基于变分自编码器的单细胞转录因子调控网络构建方法，包括步骤s1-s5：
[0070]
s1、根据scatac-seq数据和scrna-seq数据，确定输入矩阵和真实标签矩阵。
[0071]
在本发明的一个可选实施例中，本发明根据测序得到的10x格式的scatac-seq数据和scrna-seq数据，使用生信分析工具进行初步的数据筛选和质量控制，确定输入矩阵和真实标签矩阵。
[0072]
步骤s1包括以下分步骤：
[0073]
s11、利用生信数据分析工具对scatac-seq数据和scrna-seq数据进行数据筛选。
[0074]
具体地，本发明利用生信数据分析工具seurat去除检出量低于测序样本数量10％的scatac-seq数据和scrna-seq数据。
[0075]
s12、利用转录因子对齐分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，得到第一输入矩阵和第二输入矩阵。
[0076]
具体地，本发明利用转录因子对齐分步骤s11中数据筛选后的scatac-seq数据，得到第一输入矩阵；本发明利用转录因子对齐分步骤s11中数据筛选后的scrna-seq数据，得到第二输入矩阵。
[0077]
s13、利用基序扫描工具扫描分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，确定真实标签矩阵。
[0078]
具体地，本发明使用基序扫描工具fimo，设定p-value为0.000001扫描数分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，确定真实标签矩阵。
[0079]
s2、构建转录因子调控网络模型，并根据步骤s1中的输入矩阵和转录因子调控网络模型确定预测标签矩阵。
[0080]
在本发明的一个可选实施例中，本发明通过构建转录因子调控网络模型，并根据步骤s1中的输入矩阵提取scatac-seq数据和scrna-seq数据的特征信息，进而确定预测标签矩阵。
[0081]
步骤s2包括以下分步骤：
[0082]
s21、构建包括第一变分自编码器、第二变分自编码器和解码器的转录因子调控网络模型。
[0083]
第一变分自编码器和第二变分自编码器相互独立但结构相同，均包括一个有偏置系数的全连接层、一个dropout层和两个无偏置系数的全连接层；解码器包括双层全连接层模块，并通过点乘的方式嵌入包含转录因子调控关系的邻接矩阵。
[0084]
s22、训练第一变分自编码器，并输入分步骤s12中的第一输入矩阵至训练后的第一变分自编码器，得到第一特征信息，表示为：
[0085][0086]
其中：为特征信息，此处为第一特征信息，tanh为双曲正切激活函数，w3为第一变分自编码器的第三全连接层神经元的权重，w2为第一变分自编码器的第二全连接层神经元的权重，dropout为dropout层，该层以50％的概率随机失活全连接层中的神经元，w1为第一变分自编码器的第一全连接层神经元的权重，xi为输入矩阵，此处为第一输入矩阵，bias为第一变分自编码器的第一全连接层神经元的偏置系数。
[0087]
具体地，scatac-seq数据处理成了第一输入矩阵。当第一变分自编码器开始训练时，本发明按照样本随机划分第一输入矩阵得到第一变分自编码器的批次输入矩阵。每一个批次的输入矩阵形状都是64个样本
×
转录因子数量的矩阵。当第一变分自编码器训练完成时，本发明输入分步骤s12中的第一输入矩阵至训练后的第一变分自编码器，得到第一特征信息。
[0088]
s23、训练第二变分自编码器，并输入分步骤s12中的第二输入矩阵至训练后的第二变分自编码器，得到第二特征信息，表示为：
[0089][0090]
其中：为特征信息，此处为第二特征信息，w6为第二变分自编码器的第三全连接层神经元的权重，w5为第二变分自编码器的第二全连接层神经元的权重，w3为第二变分自编码器的第一全连接层神经元的权重，xi为输入矩阵，此处为第二输入矩阵，biasr为第二变分自编码器的第一全连接层神经元的偏置系数。
[0091]
具体地，scrna-seq数据处理成了与第一输入矩阵相同形状的第二输入矩阵。当第二变分自编码器开始训练时，本发明按照样本随机划分第二输入矩阵得到第二变分自编码器的批次输入矩阵。每一个批次的输入矩阵形状都是64个样本
×
转录因子数量的矩阵。当第二变分自编码器训练完成时，本发明输入分步骤s12中的第二输入矩阵至训练后的第二变分自编码器，得到第二特征信息。
[0092]
s24、训练解码器，并通过训练后的解码器得到预测标签矩阵，表示为：
[0093][0094]
其中：y
ip
为预测标签矩阵，a为包含转录因子调控关系的邻接矩阵，w8为解码器的第二层全连接层神经元的权重，w7为解码器的第一层全连接层神经元的权重，xi为输入矩
阵，此处包括第一输入矩阵和第二输入矩阵，为特征信息，此处包括第一特征信息和第二特征信息。
[0095]
本发明分别通过真实数据和伪数据依次对解码器进行训练，伪数据表示为：
[0096][0097]
其中：为伪数据，θ
d*
为解码器的内部更新参数，为特征信息，此处包括第一特征信息和第二特征信息。
[0098]
解码器包括双层全连接层模块。相比用两个独立的模块分别进行计算，双层全连接层模块相当于共享了不同解码器的参数。本发明分别使用真实数据和伪数据训练双层全连接层模块，以匹配不同数据的特征分布，并进行特征转换。
[0099]
s3、根据步骤s1中的真实标签矩阵和步骤s2中的预测标签矩阵，确定转录因子调控网络模型损失。
[0100]
在本发明的一个可选实施例中，本发明计算转录因子调控网络模型的交叉熵损失函数、转录因子调控网络模型的kl散度，并根据步骤s1中的真实标签矩阵和步骤s2中的预测标签矩阵，确定对应的均方误差函数，进而确定转录因子调控网络模型损失。
[0101]
步骤s3包括以下分步骤：
[0102]
s31、计算转录因子调控网络模型的交叉熵损失函数。
[0103]
s32、计算转录因子调控网络模型的kl散度。
[0104]
s33、计算步骤s1中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数。
[0105]
s34、根据分步骤s31中转录因子调控网络模型的交叉熵损失函数、分步骤s32中转录因子调控网络模型的kl散度和分步骤s33中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数，计算转录因子调控网络模型损失。
[0106]
本发明计算转录因子调控网络模型损失，表示为：
[0107][0108]
其中：l为转录因子调控网络模型损失，为转录因子调控网络模型的交叉熵损失函数，xi为输入矩阵，y
ip
为预测标签矩阵，β为，kl(μr,σr)||(μa,σa)为转录因子调控网络模型的kl散度，μr为scrna-seq数据特征分布的均值，σr位scrna-seq数据特征分布的标准差，μa为scatac-seq数据特征分布的均值，σa为scatac-seq数据特征分布的标准差，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵。
[0109]
s4、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对转录因子调控网络模型的内部参数进行更新。
[0110]
在本发明的一个可选实施例中，本发明根据步骤s3中的转录因子调控网络模型损失计算梯度，并采用双级优化策略对转录因子调控网络模型的内部参数进行更新。本发明对转录因子调控网络模型的内部参数进行更新的过程具体为：转录因子调控网络模型逐层计算梯度，通过反向传播将梯度返回转录因子调控网络模型的每一层；转录因子调控网络模型通过反向传播算法，可以自动计算模型内部参数对于损失函数的梯度，并利用梯度优
化算法进行参数更新，从而使得模型逐步学习和调整以使得转录因子调控网络模型损失最小化。
[0111]
步骤s4包括以下分步骤：
[0112]
s41、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对编码器的内部参数进行更新。
[0113]
步骤s41包括以下分步骤：
[0114]
s411、计算编码器的内部更新参数，表示为：
[0115][0116]
其中：θ
e*
为编码器的内部更新参数，为当转录因子调控网络模型的交叉熵损失最小时所取得的编码器内部参数，v
l
为第一数据集，vi为第一数据集中的元素，为转录因子调控网络模型的交叉熵损失函数，为特征信息，y
ip
为预测标签矩阵。
[0117]
具体地，本发明采用双级优化策略会计算编码器的内部更新参数两次。
[0118]
s412、根据步骤s3中的转录因子调控网络模型损失计算编码器的反向传播梯度。
[0119]
具体地，本发明采用采用双级优化策略，通过分步骤s411中计算的两次编码器的内部更新参数，并根据步骤s3中的转录因子调控网络模型损失计算两次编码器的反向传播梯度。
[0120]
s413、判断分步骤s412中的反向传播梯度是否下降；若是则将编码器的内部参数更新为分步骤s411中的计算值，否则对编码器的内部参数进行更新并返回分步骤s411；
[0121]
对编码器的内部参数进行更新，表示为：
[0122][0123]
其中：θ
′e为更新后的编码器的内部参数，为编码器的学习速率，为编码器的反向传播梯度。
[0124]
具体地，本发明通过判断两次编码器的反向传播梯度是否下降，即第二次编码器的反向传播梯度是否小于第一次编码器的反向传播梯度，对编码器的内部参数进行更新；若是则将编码器的内部参数更新为分步骤s411中第二次编码器的内部更新参数的计算值，否则采用上述更新公式对编码器的内部参数进行更新并返回分步骤s411。
[0125]
s42、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对解码器器的内部参数进行更新。
[0126]
步骤s42包括以下分步骤：
[0127]
s421、计算解码器的内部更新参数，表示为：
[0128][0129]
其中：θ
d*
为解码器的内部更新参数，为当真实标签矩阵和预测标签矩阵的均方误差最小时所取得的解码器内部参数，v
p
为第二数据集，vj为第二数据集中的元素，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵，yip
为预测标签矩阵。
[0130]
具体地，本发明采用双级优化策略会计算解码器的内部更新参数两次。
[0131]
s422、根据步骤s3中的转录因子调控网络模型损失计算解码器的反向传播梯度。
[0132]
具体地，本发明采用采用双级优化策略，通过分步骤s421中计算的两次解码器的内部更新参数，并根据步骤s3中的转录因子调控网络模型损失计算两次解码器的反向传播梯度。
[0133]
s423、判断分步骤s422中的反向传播梯度是否下降；若是则将解码器的内部参数更新为分步骤s421中的计算值，否则对解码器的内部参数进行更新并返回分步骤s421；
[0134]
对解码器的内部参数进行更新，表示为：
[0135][0136]
其中：θ
′d为更新后的解码器的内部参数，为解码器的学习速率，为解码器的反向传播梯度。
[0137]
具体地，本发明通过判断两次解码器的反向传播梯度是否下降，即第二次解码器的反向传播梯度是否小于第一次解码器的反向传播梯度，对解码器的内部参数进行更新；若是则将解码器的内部参数更新为分步骤s421中第二次解码器的内部更新参数的计算值，否则采用上述更新公式对编码器的内部参数进行更新并返回分步骤s421。
[0138]
s5、采用步骤s4中更新内部参数后的转录因子调控网络模型构建转录因子调控网络。
[0139]
在本发明的一个可选实施例中，本发明通过输入处理scatac-seq数据后得到的第一输入矩阵和处理scrna-seq数据后得到的第二输入矩阵至步骤s4中更新内部参数后的转录因子调控网络模型，可以确定转录因子调控关系，进而构建转录因子调控网络。
[0140]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0141]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0142]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0143]
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内
容不应理解为对本发明的限制。
[0144]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

技术特征：
1.基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，包括以下步骤：s1、根据scatac-seq数据和scrna-seq数据，确定输入矩阵和真实标签矩阵；s2、构建转录因子调控网络模型，并根据步骤s1中的输入矩阵和转录因子调控网络模型确定预测标签矩阵；s3、根据步骤s1中的真实标签矩阵和步骤s2中的预测标签矩阵，确定转录因子调控网络模型损失；s4、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对转录因子调控网络模型的内部参数进行更新；s5、采用步骤s4中更新内部参数后的转录因子调控网络模型构建转录因子调控网络。2.根据权利要求1所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s1包括以下分步骤：s11、利用生信数据分析工具对scatac-seq数据和scrna-seq数据进行数据筛选；s12、利用转录因子对齐分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，得到第一输入矩阵和第二输入矩阵；s13、利用基序扫描工具扫描分步骤s11中数据筛选后的scatac-seq数据和scrna-seq数据，确定真实标签矩阵。3.根据权利要求2所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s2包括以下分步骤：s21、构建包括第一变分自编码器、第二变分自编码器和解码器的转录因子调控网络模型；s22、训练第一变分自编码器，并输入分步骤s12中的第一输入矩阵至训练后的第一变分自编码器，得到第一特征信息，表示为：其中：为特征信息，此处为第一特征信息，tanh为双曲正切激活函数，w3为第一变分自编码器的第三全连接层神经元的权重，w2为第一变分自编码器的第二全连接层神经元的权重，dropout为dropout层，该层以50％的概率随机失活全连接层中的神经元，w1为第一变分自编码器的第一全连接层神经元的权重，x
i
为输入矩阵，此处为第一输入矩阵，bias为第一变分自编码器的第一全连接层神经元的偏置系数；s23、训练第二变分自编码器，并输入分步骤s12中的第二输入矩阵至训练后的第二变分自编码器，得到第二特征信息，表示为：其中：为特征信息，此处为第二特征信息，w6为第二变分自编码器的第三全连接层神经元的权重，w5为第二变分自编码器的第二全连接层神经元的权重，w3为第二变分自编码器
的第一全连接层神经元的权重，x
i
为输入矩阵，此处为第二输入矩阵，bias
r
为第二变分自编码器的第一全连接层神经元的偏置系数；s24、训练解码器，并通过训练后的解码器得到预测标签矩阵，表示为：其中：y
ip
为预测标签矩阵，a为包含转录因子调控关系的邻接矩阵，w8为解码器的第二层全连接层神经元的权重，w7为解码器的第一层全连接层神经元的权重，x
i
为输入矩阵，此处包括第一输入矩阵和第二输入矩阵，为特征信息，此处包括第一特征信息和第二特征信息。4.根据权利要求3所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，在分步骤s21中，第一变分自编码器和第二变分自编码器相互独立但结构相同，均包括一个有偏置系数的全连接层、一个dropout层和两个无偏置系数的全连接层；解码器包括双层全连接层模块，并通过点乘的方式嵌入包含转录因子调控关系的邻接矩阵。5.根据权利要求1所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，在分步骤s24中，分别通过真实数据和伪数据依次对解码器进行训练，伪数据表示为：其中：为伪数据，γ
d*
为解码器的内部更新参数，为特征信息，此处包括第一特征信息和第二特征信息。6.根据权利要求1所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s3包括以下分步骤：s31、计算转录因子调控网络模型的交叉熵损失函数；s32、计算转录因子调控网络模型的kl散度；s33、计算步骤s1中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数；s34、根据分步骤s31中转录因子调控网络模型的交叉熵损失函数、分步骤s32中转录因子调控网络模型的kl散度和分步骤s33中真实标签矩阵和步骤s2中预测标签矩阵的均方误差函数，计算转录因子调控网络模型损失。7.根据权利要求6所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，在分步骤s34中，计算转录因子调控网络模型损失，表示为：其中：l为转录因子调控网络模型损失，为转录因子调控网络模型的交叉熵损失函数，x
i
为输入矩阵，y
ip
为预测标签矩阵，β为，kl(μ
r
,σ
r
)||(μ
a
,σ
a
)为转录因子调控网络模型的kl散度，μ
r
为scrna-seq数据特征分布的均值，σ
r
位scrna-seq数据特征分布的标准差，μ
a
为scatac-seq数据特征分布的均值，σ
a
为scatac-seq数据特征分布的标准差，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵。
8.根据权利要求1所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s4包括以下分步骤：s41、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对编码器的内部参数进行更新；s42、根据步骤s3中的转录因子调控网络模型损失，并采用双级优化策略对解码器器的内部参数进行更新。9.根据权利要求8所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s41包括以下分步骤：s411、计算编码器的内部更新参数，表示为：其中：θ
e*
为编码器的内部更新参数，为当转录因子调控网络模型的交叉熵损失最小时所取得的编码器内部参数，v
l
为第一数据集，v
i
为第一数据集中的元素，为转录因子调控网络模型的交叉熵损失函数，为特征信息，y
ip
为预测标签矩阵；s412、根据步骤s3中的转录因子调控网络模型损失计算编码器的反向传播梯度；s413、判断分步骤s412中的反向传播梯度是否下降；若是则将编码器的内部参数更新为分步骤s411中的计算值，否则对编码器的内部参数进行更新并返回分步骤s411；对编码器的内部参数进行更新，表示为：其中：θ
′
e
为更新后的编码器的内部参数，为编码器的学习速率，为编码器的反向传播梯度。10.根据权利要求8所述的基于变分自编码器的单细胞转录因子调控网络构建方法，其特征在于，步骤s42包括以下分步骤：s421、计算解码器的内部更新参数，表示为：其中：θ
d*
为解码器的内部更新参数，为当真实标签矩阵和预测标签矩阵的均方误差最小时所取得的解码器内部参数，v
p
为第二数据集，v
j
为第二数据集中的元素，mse(y,y
ip
)为真实标签矩阵和预测标签矩阵的均方误差函数，y为真实标签矩阵，y
ip
为预测标签矩阵；s422、根据步骤s3中的转录因子调控网络模型损失计算解码器的反向传播梯度；s423、判断分步骤s422中的反向传播梯度是否下降；若是则将解码器的内部参数更新为分步骤s421中的计算值，否则对解码器的内部参数进行更新并返回分步骤s421；对解码器的内部参数进行更新，表示为：
其中：θ
′
d
为更新后的解码器的内部参数，为解码器的学习速率，为解码器的反向传播梯度。

技术总结
本发明公开了基于变分自编码器的单细胞转录因子调控网络构建方法，涉及基因转录调控领域，该方法包括：根据scATAC-seq数据和scRNA-seq数据，确定输入矩阵和真实标签矩阵；构建转录因子调控网络模型，并根据输入矩阵和转录因子调控网络模型确定预测标签矩阵；根据真实标签矩阵和预测标签矩阵，确定转录因子调控网络模型损失；采用双级优化策略对转录因子调控网络模型的内部参数进行更新；采用更新内部参数后的转录因子调控网络模型构建转录因子调控网络。本发明能提取并对齐单细胞ATAC-seq和单细胞RNA-seq数据中的高阶特征，结构化转录因子的调控关系，进而更全面准确地构建转录因子调控网络。录因子调控网络。录因子调控网络。


技术研发人员：张永清 牛颢 龙树全 邹权 龚美琴 王茂丞 何宇辰 丁春利 杨显华
受保护的技术使用者：四川省计算机研究院
技术研发日：2023.07.12
技术公布日：2023/10/15