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

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

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

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

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

步骤2)中构建的随机非线性状态空间模型方法如下：式(1)中，xt,i为第i基因在t时刻真实的基因表达值，yt,i为第i基因在t时刻的测量表达值，其中i∈[1,n]，n为基因的个数，Ci为第i个基因的衰减率，Gij为第j个基因对第i个基因的调控作用，其中j∈[1,n]，则Gij表示为gi＝[gi1,gi2,...,gi(i‑1),0,gi+1,...,gin]，即第i个基因受到除自身外的其他所有基因的调控，vi为过程噪声，wi为测量噪声，f(x(t))是系统模型中的非线性方程，表示为：则n基因中任一基因的随机非线性状态空间模型为：

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

3‑1)定义非线性状态空间模型中参数：Ξ＝[C,G,ymis,x1:n,S,R,α]，其中：S为过程噪声v的精度，R为测量噪声w的精度，ymis＝{ym1,ym2,...,ymα}表示在时刻{m1,m2,...,mα}丢失的基因表达测量值，yobs＝{ya1,ya2,...,yaα}表示为在时刻{a1,a2,...,aα}基因表达的测量值，则对于参数的估计表示为在给定观测数据下评估隐变量的后验分P(Ξ|yobs)；

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

3.如权利要求2所述的一种基于变分贝叶斯的基因调控网络结构辨识方法，其特征在于，步骤3‑1)中在给定观测数据下评估隐变量的后验分布P(Ξ|yobs)的具体步骤如下：基因表达的测量值yobs的对数似然函数表示为：ln p(yobs)＝∫q(Ξ)lnp(yobs)dΞ               (4)式(4)中，q(Ξ)为任意概率密度函数；

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

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

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

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

3‑2‑1‑1)对于具有状态空间模型的粒子滤波器，其状态的分布近似为：式(6)中，δ(·)为Dirac delta函数，wt,j是对于每一个采样点 的归一化权重，且在t＝N时刻，xN的后验分布与粒子滤波的结果相同，且

3‑2‑1‑2)分别计算wt,j, 的值：令t‑1时刻的状态近似用式(7)表示，则t时刻的采样点表示为：且t时刻的权重相应表示为：

在丢失数据点处，采用一个简单的预测方法来估计状态，则：式(9)中，权重wt,j与权重wt‑1,j相同， 的值从(7)中获得；

T

令Θ＝[C,G,S,R,α]，根据贝叶斯法则：式(10)中，p(xt+1|y1:t,Θ)近似地表示为：假设在t+1时刻，xt+1的后验分布近似为：则，xt的后验分布表示为：

式(13)中，

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

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

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

过程噪声wi与测量噪声vi服从零均值的高斯分布：式(17)中，ri与si表示为正态分布的精度；

假设在 中存在不确定性，且 的先验分布为高斯‑伽马分布，ri与si的先验分布为伽马分布，则参数的联合先验分布表示为:‑1

式(18)中，α 是 中每个参数的协方差，In×n是n×n的单位矩阵，G(·)为伽马分布，a0,b0是两个常值超参数；

则 服从高斯分布：

式(19)中：

式(20)中，＜·＞代表了对于 的期望值，

3‑2‑3)更新q(ri)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(ri)服从一个伽马分布：式(21)中，＜·＞代表了对于 的期望值；

3‑2‑4)更新q(si)：

基于式(15)和(18)，采用变分贝叶斯的方法，q(si)服从一个伽马分布：式(22)中，Nα为可获得的观测值数，＜·＞代表了对于 的期望值；

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

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

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

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

4‑1)基于式(15)所述的基因调控网络GRN的模型，为识别第i个基因的基因相互表达作用，用平滑的基因表达值去建立一个新的基因调控网络GRN模型：式(24)中，XGBoost方法为GRN推断提供了一个增强型决策树模型，则m棵树的目标函数为：式(25)中，L是损失函数，Tm‑1与Tm分别为m‑1颗和m颗树的输出，bt,ct为处L的一阶和二阶梯度，M为叶的数量，wj为相应的权重，γ和λ为两常数；

采用采用贪婪算法优化目标函数，则将一片叶子分为两片叶的标准为：式(26)中，GL和GR为分裂后左右结点的bt之和,HR和HL为分裂后左右结点的ct之和；

则采用第m颗树的重要性得分去评价给定基因与其他基因之间的调控关系被计算为：式(27)中， 表示为在m棵树中第i表基因与第j个基因之间的调控链的重要性得分，表示为m棵树的节点数，如果在第k个结点，选择了第j个基因，则o(k,j)＝1，否则，o(k,j)＝0；

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

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

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