CN114360641A - 一种基于变分贝叶斯的基因调控网络结构辨识方法 - Google Patents

Abstract

本发明提供一种基于变分贝叶斯的基因调控网络结构辨识方法，根据不完整和有噪声的基因表达时间序列数据中考虑基因调控网络GRN的结构推断，用含有未知噪声信息的随机非线性状态空间模型，描述了基因表达数据中的动态行为，采用变分贝叶斯VB框架来同时估计参数和基因表达水平，通过生成预测值，可以很容易地处理缺失的观测值；考虑到GRN的稀疏性，利用极端梯度增强树对平滑后的基因数据进行建模，并通过树模型中的重要性得分来识别基因间的调控相互作用。该方法能在观测值缺失的情况下，有效地恢复GRN的调控相互作用，并优于现有的GRN识别方法。

Description

一种基于变分贝叶斯的基因调控网络结构辨识方法

技术领域

本发明涉及基因调控网络识别领域，特别是一种基于变分贝叶斯的基因调控网络结构辨识方法。

背景技术

活细胞中基因调控网络GRN的鉴定是了解基因相互作用进行生物学研究的重要问题，是病变基因鉴定、药物开发、代谢调控等的基础。在过去的几十年里，高通量基因表达测量技术的产生提供了大量的生物数据，并使从基因表达数据中推断GRN成为可能。对于GRN的推理方法有很多，这些方法可以松散地分为两种类型。第一种类型是利用静态表达式数据构建模型，并利用聚类、互信息和相关分析进行GRN推理。一般来说，静态基因表达数据是稳态下的基因表达水平。研究基因调控相互作用的一种更精确的方法是引入一些环境扰动，并测量基因表达时间序列数据。在这种情况下，GRN是由动态数据建立起来的。与基于静态数据的GRN相比，基于动态数据的模型更加精确，这一问题受到了越来越多的关注。

注基于状态空间模型的基因调控网络GRN，其可以明确地描述基因的表达过程。由于基因表达过程是一个高度非线性的过程，当用线性模型来描述GRN时，存在一定的局限性。因此，在最近的研究中，非线性模型一直是主要使用的模型。除了非线性特征外，随机行为是GRN的另一个固有特性。在状态空间模型中，随机行为可以很容易地描述为噪声。此外，数据不完整是几乎所有基于数据的建模问题所遇到的常见现象。由于GRN推理的观测值数相对较少，因此在数据缺失时，推理精度将显著降低。

发明内容

本发明的目的就是提供一种基于变分贝叶斯的基因调控网络结构辨识方法，本发明通过对输入数据的处理，识别并预测基因表达水平。

本发明的目的是通过这样的技术方案实现的，它包括有以下步骤：

1)数据采集：从DREAM4平台的数据库中采集典型基因调控网络GRN的基因表达数据，获得采集数据；

2)构建模型：采用含有未知噪声的随机非线性状态空间模型来描述步骤1)中的采集数据中基因序列的表达过程；

3)模型参数估计：采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据；

4)得到GRN结构：利用极端梯度提升XGBoost方法建立基因调控网络GRN的决策树模型，辨别基因间的相互作用关系，得到最终的基因调控网络GRN结构。

进一步，步骤2)中构建的随机非线性状态空间模型方法如下：

式(1)中，x_t,i为第i基因在t时刻真实的基因表达值，y_t,i为第i基因在t时刻的测量表达值，其中i∈[1,n]，n为基因的个数，C_i为第i个基因的衰减率，G_ij为第j个基因对第i个基因的调控作用，其中j∈[1,n]，则G_ij表示为g_i＝[g_i1,g_i2,...,g_i(i-1),0,g_i+1,...,g_in]，即第i个基因受到除自身外的其他所有基因的调控，v_i为过程噪声，w_i为测量噪声，f(x(t))是***模型中的非线性方程，表示为：

则n基因中任一基因的随机非线性状态空间模型为：

进一步，步骤3)中采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据，具体步骤如下：

3-1)定义非线性状态空间模型中参数：Ξ＝[C,G,y_mis,x_1:n,S,R,α]^T，其中：S为过程噪声v的精度，R为测量噪声w的精度，y_mis＝{y_m1,y_m2,...,y_mα}表示在时刻{m₁,m₂,...,m_α}丢失的基因表达测量值，y_obs＝{y_a1,y_a2,...,y_aα}表示为在时刻{a₁,a₂,...,a_α}基因表达的测量值，则对于参数的估计可以表示为在给定观测数据下评估隐变量的后验分P(Ξ|y_obs)；

3-2)通过变分贝叶斯的方法对随机非线性状态空间模型中的参数C、G、y_mis、x_1:n、S、R、α进行估计。

进一步，步骤3-1)中在给定观测数据下评估隐变量的后验分布P(Ξ|y_obs)的具体步骤如下：

基因表达的测量值y_obs的对数似然函数表示为：

lnp(y_obs)＝∫q(Ξ)lnp(y_obs)dΞ (4)

式(4)中，q(Ξ)为任意概率密度函数；

将对数似然函数写成J+KL的形式：

式(5)中，KL≥0，当且仅当q(Ξ)＝p(Ξ|y_obs)时等号成立，J为对数似然函数lnp(y_obs)的下界，这时计算隐变量的后验分布等价于计算argmaxJ。

进一步，步骤3-2)中通过变分贝叶斯的方法对随机非线性状态空间模型中的参数的具体步骤为：

3-2-1)用粒子平滑器更新q(x_1:n)：

3-2-1-1)对于具有状态空间模型的粒子滤波器，其状态的分布近似为：

式(6)中，δ(·)为Dirac delta函数，w_t,j是对于每一个采样点

的归一化权重，且

在t＝N时刻，x_N的后验分布与粒子滤波的结果相同，且

3-2-1-2)分别计算w_t,j,

的值：

令t-1时刻的状态近似用式(7)表示，则t时刻的采样点可以表示为：

且t时刻的权重相应表示为：

在丢失数据点处，采用一个简单的预测方法来估计状态，则：

式(9)中，权重w_t,j与权重w_t-1,j相同，

的值从(7)中获得；

令Θ＝[C,G,S,R,α]^T，根据贝叶斯法则：

式(10)中，p(x_t+1|y_1:t,Θ)可以近似地表示为：

假设在t+1时刻，x_t+1的后验分布近似为：

则，x_t的后验分布表示为：

式(13)中，

3-2-1-3)令t＝t-1重复步骤3-2-1-2)，直到t＝1，完成x_1:n的更新；3-2-2)更新q(G)：

基于式(1)，第i个基因的表达过程表示为：

式(15)中：

过程噪声w_i与测量噪声v_i服从零均值的高斯分布：

式(17)中，r_i与s_i表示为正态分布的精度。

假设在

中存在不确定性，且

的先验分布为高斯-伽马分布，r_i与s_i的先验分布为伽马分布，则参数的联合先验分布表示为:

式(18)中，α^-1是

中每个参数的协方差，I_n×n是n×n的单位矩阵，G(·)为伽马分布，a₀,b₀是两个常值超参数；

则

服从高斯分布：

式(19)中：

式(20)中，<·>代表了对于

的期望值，

3-2-3)更新q(r_i)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(r_i)服从一个伽马分布：

式(21)中，<·>代表了对于

的期望值；

3-2-4)更新q(s_i)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(s_i)服从一个伽马分布：

式(22)中，N_α为可获得的观测值数，<·>代表了对于

的期望值；

3-2-5)更新q(α)：

基于式(15)和(18)，使用变分贝叶斯的方法，q(α)服从一个伽马分布：

3-2-6)采用步骤3-2-1)至步骤3-2-5)所述的方法更新模型参数Θ，用模型参数Θ的期望值更新网络模型，直到迭代次数最大。

进一步，步骤4)中得到最终的基因调控网络GRN结构，具体步骤如下：

4-1)基于式(15)所述的基因调控网络GRN的模型，为识别第i个基因的基因相互表达作用，用平滑的基因表达值去建立一个新的基因调控网络GRN模型：

式(24)中，XGBoost方法为GRN推断提供了一个增强型决策树模型，则m棵树的目标函数为：

式(25)中，L是损失函数，T_m-1与T_m分别为m-1颗和m颗树的输出，b_t,c_t为

处L的一阶和二阶梯度，M为叶的数量，w_j为相应的权重，γ和λ为两常数；

采用采用贪婪算法优化目标函数，则将一片叶子分为两片叶的标准为：

式(26)中，G_L和G_R为***后左右结点的b_t之和,H_R和H_L为***后左右结点的c_t之和；

则采用第m颗树的重要性得分去评价给定基因与其他基因之间的调控关系可以被计算为：

式(27)中，

表示为在m棵树中第i表基因与第j个基因之间的调控链的重要性得分，

表示为m棵树的节点数，如果在第k个结点，选择了第j个基因，则o(k,j)＝1，否则，o(k,j)＝0；

则，第j个基因的重要性总体得分表示为：

式(28)中，N_tree为集合树的数量；

4-2)基于步骤4-1)所得的重要性总体得分，对第i个基因与其他基因之间的调控关系进行排序，并通过选择w_i,j的阈值来确定其中一个基因是否对给定第i个基因有影响，得到最终的基因调控网络GRN结构。

由于采用了上述技术方案，本发明具有如下的优点：

1、本申请通过构建一个具有未知噪声信息的非线性随机状态空间模型能够对基因调控网络中的高度非线性和随机性进行识别，适用范围广。

2、本申请提出采用变分贝叶斯的方法进行非线性的GRN推理，该方法提供了状态估计与噪声估计，从而提高了模型识别精度。

3、本申请考虑了存在数据缺失时的GRN推理，提出采用预测数据来代替数据缺失，从而有效地证明实验的准确性。

4、本申请采用VB方法来处理不确定性，并将XGBoost方法用于探索基因相互作用的稀疏性，并在人工基因表达数据和真实数据上对结果进行验证。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书和权利要求书来实现和获得。

附图说明

本发明的附图说明如下：

图1为本发明的流程示意图。

图2为本发明的采用VB方法估计的基因表达水平。

图3为本发明的在数据缺失下采用VB方法估计的基因表达水平。

图4为本发明的估计GAL基因网络。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

一种基于变分贝叶斯的基因调控网络结构辨识方法，具体步骤如下：

1)数据采集：从DREAM4平台的数据库中采集典型基因调控网络GRN的基因表达数据，获得采集数据；

2)构建模型：采用含有未知噪声的随机非线性状态空间模型来描述步骤1)中的采集数据中基因序列的表达过程，方法如下：

式(29)中，x_t,i为第i基因在t时刻真实的基因表达值，y_t,i为第i基因在t时刻的测量表达值，其中i∈[1,n]，n为基因的个数，C_i为第i个基因的衰减率，G_ij为第j个基因对第i个基因的调控作用，其中j∈[1,n]，则G_ij表示为g_i＝[g_i1,g_i2,...,g_i(i-1),0,g_i+1,...,g_in]，即第i个基因受到除自身外的其他所有基因的调控，v_i为过程噪声，w_i为测量噪声，f(x(t))是***模型中的非线性方程，表示为：

则n基因中任一基因的随机非线性状态空间模型为：

3)模型参数估计：采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据；具体步骤如下：

3-1)定义非线性状态空间模型中参数：Ξ＝[C,G,y_mis,x_1:n,S,R,α]^T，其中：S为过程噪声v的精度，R为测量噪声w的精度，y_mis＝{y_m1,y_m2,...,y_mα}表示在时刻{m₁,m₂,...,m_α}丢失的基因表达测量值，y_obs＝{y_a1,y_a2,...,y_aα}表示为在时刻{a₁,a₂,...,a_α}基因表达的测量值，则对于参数的估计可以表示为在给定观测数据下评估隐变量的后验分P(Ξ|y_obs)；具体方法如下：

基因表达的测量值y_obs的对数似然函数表示为：

lnp(y_obs)＝∫q(Ξ)lnp(y_obs)dΞ (32)

式(32)中，q(Ξ)为任意概率密度函数；

将对数似然函数写成J+KL的形式：

式(33)中，KL≥0，当且仅当q(Ξ)＝p(Ξ|y_obs)时等号成立，J为对数似然函数lnp(y_obs)的下界，这时计算隐变量的后验分布等价于计算argmaxJ。

3-2)通过变分贝叶斯的方法对随机非线性状态空间模型中的参数C、G、y_mis、x_1:n、S、R、α进行估计，具体步骤如下：

3-2-1)用粒子平滑器更新q(x_1:n)：

3-2-1-1)对于具有状态空间模型的粒子滤波器，其状态的分布近似为：

式(34)中，δ(·)为Dirac delta函数，w_t,j是对于每一个采样点

的归一化权重，且

在t＝N时刻，x_N的后验分布与粒子滤波的结果相同，且

3-2-1-2)分别计算w_t,j,

的值：

令t-1时刻的状态近似用式(7)表示，则t时刻的采样点可以表示为：

且t时刻的权重相应表示为：

在丢失数据点处，采用一个简单的预测方法来估计状态，则：

式(37)中，权重w_t,j与权重w_t-1,j相同，

的值从(7)中获得；

令Θ＝[C,G,S,R,α]^T，根据贝叶斯法则：

式(38)中，p(x_t+1|y_1:t,Θ)可以近似地表示为：

假设在t+1时刻，x_t+1的后验分布近似为：

则，x_t的后验分布表示为：

式(41)中，

3-2-1-3)令t＝t-1重复步骤3-2-1-2)，直到t＝1，完成x_1:n的更新；

3-2-2)更新q(G)：

基于式(29)，第i个基因的表达过程表示为：

式(43)中：

过程噪声w_i与测量噪声v_i服从零均值的高斯分布：

式(45)中，r_i与s_i表示为正态分布的精度。

假设在

中存在不确定性，且

的先验分布为高斯-伽马分布，r_i与s_i的先验分布为伽马分布，则参数的联合先验分布表示为:

式(46)中，α^-1是

中每个参数的协方差，I_n×n是n×n的单位矩阵，G(·)为伽马分布，a₀,b₀是两个常值超参数；

则

服从高斯分布：

式(19)中：

式(48)中，<·>代表了对于

的期望值，

3-2-3)更新q(r_i)：

基于式(43)和(46)，采用变分贝叶斯的方法，q(r_i)服从一个伽马分布：

式(49)中，<·>代表了对于

的期望值；

3-2-4)更新q(s_i)：

基于式(43)和(46)，采用变分贝叶斯的方法，q(s_i)服从一个伽马分布：

式(50)中，N_α为可获得的观测值数，<·>代表了对于

的期望值；

3-2-5)更新q(α)：

基于式(43)和(46)，使用变分贝叶斯的方法，q(α)服从一个伽马分布：

3-2-6)采用步骤3-2-1)至步骤3-2-5)所述的方法更新模型参数Θ，用模型参数Θ的期望值更新网络模型，直到迭代次数最大。

4)得到GRN结构：利用极端梯度提升XGBoost方法建立基因调控网络GRN的决策树模型，辨别基因间的相互作用关系，得到最终的基因调控网络GRN结构，具体步骤如下：

4-1)基于式(43)所述的基因调控网络GRN的模型，为识别第i个基因的基因相互表达作用，用平滑的基因表达值去建立一个新的基因调控网络GRN模型：

式(52)中，XGBoost方法为GRN推断提供了一个增强型决策树模型，则m棵树的目标函数为：

式(53)中，L是损失函数，T_m-1与T_m分别为m-1颗和m颗树的输出，b_t,c_t为

处L的一阶和二阶梯度，M为叶的数量，w_j为相应的权重，γ和λ为两常数；

采用采用贪婪算法优化目标函数，则将一片叶子分为两片叶的标准为：

式(54)中，G_L和G_R为***后左右结点的b_t之和,H_R和H_L为***后左右结点的c_t之和；

则采用第m颗树的重要性得分去评价给定基因与其他基因之间的调控关系可以被计算为：

式(55)中，

表示为在m棵树中第i表基因与第j个基因之间的调控链的重要性得分，

表示为m棵树的节点数，如果在第k个结点，选择了第j个基因，则o(k,j)＝1，否则，o(k,j)＝0；

则，第j个基因的重要性总体得分表示为：

式(56)中，N_tree为集合树的数量；

4-2)基于步骤4-1)所得的重要性总体得分，对第i个基因与其他基因之间的调控关系进行排序，并通过选择w_i,j的阈值来确定其中一个基因是否对给定第i个基因有影响，得到最终的基因调控网络GRN结构。

5)数据模拟仿真：在InSilico_Size10提供的基因网络与酿酒酵母GAL基因网络中验证模型方法的准确性，其步骤如下：

5-1)首先计算假阳性率(FRR)和真阳性率(TRR)，并将其定义为：

式(57)中，TP为真阳性数，TN为真阴性数，FP为假阳性数，FN为假阴性数；基于这四个指标，Precision和Recall定义为：

基于这些指标，计算指标AUROC(FRR-TRR下方的面积)与AUPR(Precision-Recall下方的面积)，如表1于表3所示为本申请所提出的方法与现有方法的性能比较。

表1现有方法与本申请所提出的方法在InSideco_Size10数据集上的性能比较

表2在数据缺失存在的情况下BiXGBoost与本申请所提出方法在InSideco_Size10数据集上的性能比较

表3在数据缺失存在的情况下现有方法与本申请所提出方法在InSideco_Size10数据集上的性能比较

5-2)在InSilico_Size10提供的基因网络与酿酒酵母GAL基因网络中验证模型方法的准确性；如表4所示为本申请在GAL基因网络上现有方法与本申请所提出方法的性能比较。

表4在GAL基因网络上现有方法与本申请所提出方法的性能比较

采用本发明提出的变分贝叶斯方法的评估曲线和测量曲线，如2所示，通过对比可以看出，VB方法可以很好拟合测量曲线，且误差较小。

如图3所示，在数据缺失存在的情况下采用本发明提出的变分贝叶斯方法评估曲线依然可以较好地拟合测量曲线，该方法克服了现有的经典GRN推理的局限性，提供了一种有效的基因表达序列推断方法，表明了本发明所提出的在数据缺失的情况下，采用变分贝叶斯的方法进行GRN识别的有效性。

本领域内的技术人员应明白，本申请的实施例可提供为方法、***、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(***)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。

Claims (6)

1.一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，具体步骤如下：

1)数据采集：从DREAM4平台的数据库中采集典型基因调控网络GRN的基因表达数据，获得采集数据；

2)构建模型：采用含有未知噪声的随机非线性状态空间模型来描述步骤1)中的采集数据中基因序列的表达过程；

3)模型参数估计：采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据；

4)得到GRN结构：利用极端梯度提升XGBoost方法建立基因调控网络GRN的决策树模型，辨别基因间的相互作用关系，得到最终的基因调控网络GRN结构。

2.如权利要求1所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤2)中构建的随机非线性状态空间模型方法如下：

式(1)中，x_t,i为第i基因在t时刻真实的基因表达值，y_t,i为第i基因在t时刻的测量表达值，其中i∈[1,n]，n为基因的个数，C_i为第i个基因的衰减率，G_ij为第j个基因对第i个基因的调控作用，其中j∈[1,n]，则G_ij表示为g_i＝[g_i1,g_i2,...,g_i(i-1),0,g_i+1,...,g_in]，即第i个基因受到除自身外的其他所有基因的调控，v_i为过程噪声，w_i为测量噪声，f(x(t))是***模型中的非线性方程，表示为：

则n基因中任一基因的随机非线性状态空间模型为：

3.如权利要求2所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤3)中采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据，具体步骤如下：

3-1)定义非线性状态空间模型中参数：Ξ＝[C,G,y_mis,x_1:n,S,R,α]^T，其中：S为过程噪声v的精度，R为测量噪声w的精度，y_mis＝{y_m1,y_m2,...,y_mα}表示在时刻{m₁,m₂,...,m_α}丢失的基因表达测量值，y_obs＝{y_a1,y_a2,...,y_aα}表示为在时刻{a₁,a₂,...,a_α}基因表达的测量值，则对于参数的估计可以表示为在给定观测数据下评估隐变量的后验分P(Ξ|y_obs)；

3-2)通过变分贝叶斯的方法对随机非线性状态空间模型中的参数C、G、y_mis、x_1:n、S、R、α进行估计。

4.如权利要求3所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤3-1)中在给定观测数据下评估隐变量的后验分布P(Ξ|y_obs)的具体步骤如下：

基因表达的测量值y_obs的对数似然函数表示为：

lnp(y_obs)＝∫q(Ξ)lnp(y_obs)dΞ (4)

式(4)中，q(Ξ)为任意概率密度函数；

将对数似然函数写成J+KL的形式：

式(5)中，KL≥0，当且仅当q(Ξ)＝p(Ξ|y_obs)时等号成立，J为对数似然函数lnp(y_obs)的下界，这时计算隐变量的后验分布等价于计算argmaxJ。

5.如权利要求3所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤3-2)中通过变分贝叶斯的方法对随机非线性状态空间模型中的参数的具体步骤为：

3-2-1)用粒子平滑器更新q(x_1:n)：

3-2-1-1)对于具有状态空间模型的粒子滤波器，其状态的分布近似为：

式(6)中，δ(·)为Dirac delta函数，w_t,j是对于每一个采样点

的归一化权重，且

在t＝N时刻，x_N的后验分布与粒子滤波的结果相同，且

3-2-1-2)分别计算w_t,j,

的值：

令t-1时刻的状态近似用式(7)表示，则t时刻的采样点可以表示为：

且t时刻的权重相应表示为：

在丢失数据点处，采用一个简单的预测方法来估计状态，则：

式(9)中，权重w_t,j与权重w_t-1,j相同，

的值从(7)中获得；

令Θ＝[C,G,S,R,α]^T，根据贝叶斯法则：

式(10)中，p(x_t+1|y_1:t，Θ)可以近似地表示为：

假设在t+1时刻，x_t+1的后验分布近似为：

则，x_t的后验分布表示为：

式(13)中，

3-2-1-3)令t＝t-1重复步骤3-2-1-2)，直到t＝1，完成x_1:n的更新；

3-2-2)更新q(G)：

基于式(1)，第i个基因的表达过程表示为：

式(15)中：

过程噪声w_i与测量噪声v_i服从零均值的高斯分布：

式(17)中，r_i与s_i表示为正态分布的精度。

假设在

中存在不确定性，且

的先验分布为高斯-伽马分布，r_i与s_i的先验分布为伽马分布，则参数的联合先验分布表示为:

式(18)中，α^-1是

中每个参数的协方差，I_n×n是n×n的单位矩阵，G(·)为伽马分布，a₀,b₀是两个常值超参数；

则

服从高斯分布：

式(19)中：

式(20)中，<·>代表了对于

的期望值，

3-2-3)更新q(r_i)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(r_i)服从一个伽马分布：

式(21)中，<·>代表了对于

的期望值；

3-2-4)更新q(s_i)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(s_i)服从一个伽马分布：

式(22)中，N_α为可获得的观测值数，<·>代表了对于

的期望值；

3-2-5)更新q(α)：

基于式(15)和(18)，使用变分贝叶斯的方法，q(α)服从一个伽马分布：

3-2-6)采用步骤3-2-1)至步骤3-2-5)所述的方法更新模型参数Θ，用模型参数Θ的期望值更新网络模型，直到迭代次数最大。

6.如权利要求5所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤4)中得到最终的基因调控网络GRN结构，具体步骤如下：

4-1)基于式(15)所述的基因调控网络GRN的模型，为识别第i个基因的基因相互表达作用，用平滑的基因表达值去建立一个新的基因调控网络GRN模型：

式(24)中，XGBoost方法为GRN推断提供了一个增强型决策树模型，则m棵树的目标函数为：

式(25)中，L是损失函数，T_m-1与T_m分别为m-1颗和m颗树的输出，b_t,c_t为

处L的一阶和二阶梯度，M为叶的数量，w_j为相应的权重，γ和λ为两常数；

采用采用贪婪算法优化目标函数，则将一片叶子分为两片叶的标准为：

式(26)中，G_L和G_R为***后左右结点的b_t之和,H_R和H_L为***后左右结点的c_t之和；

则采用第m颗树的重要性得分去评价给定基因与其他基因之间的调控关系可以被计算为：

式(27)中，

表示为在m棵树中第i表基因与第j个基因之间的调控链的重要性得分，

表示为m棵树的节点数，如果在第k个结点，选择了第j个基因，则o(k,j)＝1，否则，o(k,j)＝0；

则，第j个基因的重要性总体得分表示为：

式(28)中，N_tree为集合树的数量；

4-2)基于步骤4-1)所得的重要性总体得分，对第i个基因与其他基因之间的调控关系进行排序，并通过选择w_i,j的阈值来确定其中一个基因是否对给定第i个基因有影响，得到最终的基因调控网络GRN结构。

Images