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

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

【技术实现步骤摘要】
一种基于变分贝叶斯的基因调控网络结构辨识方法


[0001]本专利技术涉及基因调控网络识别领域，特别是一种基于变分贝叶斯的基因调控网络结构辨识方法。

技术介绍

[0002]活细胞中基因调控网络GRN的鉴定是了解基因相互作用进行生物学研究的重要问题，是病变基因鉴定、药物开发、代谢调控等的基础。在过去的几十年里，高通量基因表达测量技术的产生提供了大量的生物数据，并使从基因表达数据中推断GRN成为可能。对于GRN的推理方法有很多，这些方法可以松散地分为两种类型。第一种类型是利用静态表达式数据构建模型，并利用聚类、互信息和相关分析进行GRN推理。一般来说，静态基因表达数据是稳态下的基因表达水平。研究基因调控相互作用的一种更精确的方法是引入一些环境扰动，并测量基因表达时间序列数据。在这种情况下，GRN是由动态数据建立起来的。与基于静态数据的GRN相比，基于动态数据的模型更加精确，这一问题受到了越来越多的关注。
[0003]注基于状态空间模型的基因调控网络GRN，其可以明确地描述基因的表达过程。由于基因表达过程是一个高度非线性的过程，当用线性模型来描述GRN时，存在一定的局限性。因此，在最近的研究中，非线性模型一直是主要使用的模型。除了非线性特征外，随机行为是GRN的另一个固有特性。在状态空间模型中，随机行为可以很容易地描述为噪声。此外，数据不完整是几乎所有基于数据的建模问题所遇到的常见现象。由于GRN推理的观测值数相对较少，因此在数据缺失时，推理精度将显著降低。

技术实现思路

[0004]本专利技术的目的就是提供一种基于变分贝叶斯的基因调控网络结构辨识方法，本专利技术通过对输入数据的处理，识别并预测基因表达水平。
[0005]本专利技术的目的是通过这样的技术方案实现的，它包括有以下步骤：
[0006]1)数据采集：从DREAM4平台的数据库中采集典型基因调控网络GRN的基因表达数据，获得采集数据；
[0007]2)构建模型：采用含有未知噪声的随机非线性状态空间模型来描述步骤1)中的采集数据中基因序列的表达过程；
[0008]3)模型参数估计：采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据；
[0009]4)得到GRN结构：利用极端梯度提升XGBoost方法建立基因调控网络GRN的决策树模型，辨别基因间的相互作用关系，得到最终的基因调控网络GRN结构。
[0010]进一步，步骤2)中构建的随机非线性状态空间模型方法如下：
[0011][0012]式(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))是系统模型中的非线性方程，表示为：
[0013][0014]则n基因中任一基因的随机非线性状态空间模型为：
[0015][0016]进一步，步骤3)中采用变分贝叶斯的方法对状态空间模型中的参数进行估计，并输出基因序列表达数据，具体步骤如下：
[0017]3‑
1)定义非线性状态空间模型中参数：Ξ＝[C,G,y
mis
,x
1:n
,S,R,α]T
，其中：S为过程噪声v的精度，R为测量噪声w的精度，y
mis
＝{y
m1
,y
m2
,...,y
mα
}表示在时刻{m1,m2,...,m
α
}丢失的基因表达测量值，y
obs
＝{y
a1
,y
a2
,...,y
aα
}表示为在时刻{a1,a2,...,a
α
}基因表达的测量值，则对于参数的估计可以表示为在给定观测数据下评估隐变量的后验分P(Ξ|y
obs
)；
[0018]3‑
2)通过变分贝叶斯的方法对随机非线性状态空间模型中的参数C、G、y
mis
、x
1:n
、S、R、α进行估计。
[0019]进一步，步骤3
‑
1)中在给定观测数据下评估隐变量的后验分布P(Ξ|y
obs
)的具体步骤如下：
[0020]基因表达的测量值y
obs
的对数似然函数表示为：
[0021]lnp(y
obs
)＝∫q(Ξ)lnp(y
obs
)dΞ
ꢀꢀ
(4)
[0022]式(4)中，q(Ξ)为任意概率密度函数；
[0023]将对数似然函数写成J+KL的形式：
[0024][0025]式(5)中，KL≥0，当且仅当q(Ξ)＝p(Ξ|y
obs
)时等号成立，J为对数似然函数lnp(y
obs
)的下界，这时计算隐变量的后验分布等价于计算argmaxJ。
[0026]进一步，步骤3
‑
2)中通过变分贝叶斯的方法对随机非线性状态空间模型中的参数的具体步骤为：
[0027]3‑2‑
1)用粒子平滑器更新q(x
1:n
)：
[0028]3‑2‑1‑
1)对于具有状态空间模型的粒子滤波器，其状态的分布近似为：
[0029][0030]式(6)中，δ(
·
)为Dirac delta函数，w
t,j
是对于每一个采样点的归一化权重，
且
[0031]在t＝N时刻，x
N
的后验分布与粒子滤波的结果相同，且
[0032]3‑2‑1‑
2)分别计算w
t,j
,的值：
[0033]令t
‑
1时刻的状态近似用式(7)表示，则t时刻的采样点可以表示为：
[0034][0035]且t时刻的权重相应表示为：
[0036][0037]在丢失数据点处，采用一个简单的预测方法来估计状态，则：
[0038][0039]式(9)中，权重w
t,j
与权重w
t
‑
1,j
相同，的值从(7)中获得；
[0040]令Θ＝[C,G,S,R,α]T
，根据贝叶斯法则：
[0本文档来自技高网...

【技术保护点】

【技术特征摘要】
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α
}表示在时刻{m1,m2,...,m
α
}丢失的基因表达测量值，y
obs
＝{y
a1
,y
a2
,...,y
aα
}表示为在时刻{a1,a2,...,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时刻的权重相应表示为：在...

【专利技术属性】
技术研发人员：刘切，王浩，李俊豪，柴毅，
申请(专利权)人：重庆大学，
类型：发明
国别省市：