1.一种基于随机变分贝叶斯的基因调控网络拓扑结构识别方法,其特征在于,具体步骤如下:
1)数据获取:获取待识别基因表达数据集,该基因表达数据集满足正态分布的高斯噪声信号;
2)构建模型:采用基于状态空间模型的动态结构函数DSF建模重构基因调控网络模型;
3)模型参数估计:采用基于随机变分贝叶斯的方法对基因调控网络模型中的参数进行估计;
4)网络拓扑结构识别:采用前向选择方法对基因调控网络的拓扑结构进行识别,用ARD变量来更新网络模型,并使用下界函数J选择模型结构,绘制基因调控网络拓扑图;
步骤2)中重构基因调控网络模型的具体结构为:基于状态空间模型,所述基因调控网络通过动态结构函数DSF进行建模,则第i基因被描述为:式中,yi(t)为第i基因在t时刻的测量表达值,yj(t)为第j基因在t时刻的测量表达值,Qij为第j基因对第i基因的调控作用,n为基因的个数,i,j∈[1,n],m为基因调控网络输入的个数,Pij为第j输入对第i基因的调控作用,uj(t)为基因调控网络的第j个输入,ei(t)为测量噪音,Qij和Pij是Q和P中的元素,Q和P是传递函数矩阵;
传递函数矩阵Q为0,则表示第i基因不受第j基因的调控,传递函数矩阵P为0,则表示第i基因不受第j输入的影响,在基因调控网络拓扑图中表现为该节点无输入;
步骤3)中对基因调控网络模型中的参数进行估计的具体步骤为:
3‑1)将对基因调控网络拓扑结构的识别转换为对基因调控网络参数的估计,表示为q(Θ),对基因调控网络参数的估计等同于计算其下界函数J的最大值;
3‑2)基于随机优化的贝叶斯变分方法计算下界函数J最大化;
3‑3)通过随机变分贝叶斯方法估计基因调控网络的参数;
步骤4)中采用前向选择方法对基因调控网络的拓扑结构进行识别的具体方法如下:使用ARD变量来更新网络模型,并使用下界函数J选择模型结构,即通过下界函数确定传递矩阵Q和P的结构:
4‑1)识别Q的结构:对于全连接的基因调控网络模型,使用步骤3)中的随机变分贝叶斯方法来估计ARD变量α;
4‑1‑1)将ARD变量αj(j=1,2,…,n)按升序排列:
4‑1‑2)令R=1,则基因调控网络模型为:
4‑1‑3)对于基因调控网络模型,使用步骤3)中的随机变分贝叶斯方法计算下界,表示R为J;
R R‑1
4‑1‑4)若J
4‑2)识别P的结构:对于全连接的基因调控网络模型,使用步骤3)中的随机变分贝叶斯方法来估计ARD变量β;
4‑2‑1)将ARD变量βj(j=1,2,…,m)按升序排列:
4‑2‑2)令R=1,则基因调控网络模型为:
4‑2‑3)对于基因调控网络模型,使用步骤3)中的随机变分贝叶斯方法计算下界,表示R为J;
R R‑1
4‑2‑4)若J
2.如权利要求1所述的一种基于随机变分贝叶斯的基因调控网络拓扑结构识别方法,其特征在于,步骤3‑1)中计算其下界函数J的最大值的具体方法为:第i个基因的测量表达值被描述为:式中,hij为第i个基因模型表示中Qij的脉冲响应,lij为第i个基因模型表示中Pij的脉冲响应,Z为足够大的正整数,用于使|hij,Z|≈0,|lij,Z|≈0;
设定噪声服从均值为零的高斯分布,即:‑1
ei(t)~N(0,τ ) (3)式中,τ为测量噪音的高斯分布的精度参数;
设定τ的先验分布为一个伽马分布,即:式中,a0,b0为两个超参数;
hij和lij的先验分布为:
式中, 和 表示TC核,αj和βj为基因调控网络拓扑结构的自动相关确定ARD变量;
设定αj和βj的先验分布为共轭先验的伽马分布:p(αj)=p(βj)=G(τ|a0,b0) (8)核矩阵 和 中的每一个元素有如下定义:式中,γ1,j和γ2,j是核矩阵的常数;
则参数表示为 ARD变量表示为α=[α1,…,T
αn],β=[β1,…,βn],令Yi=[yi(1),…,yi(N)];
基于贝叶斯方法的网络识别是似然函数最大化:式中,<·>表示期望算子,q<·>表示后验分布;
T
令Θ=[h,l,τ,α,β],则基因调控网络第i个基因输出的对数似然函数有:式中, 为KL散度, 为似然函数的下限,用J表示,q(Θ)的计算等同于计算下限J的最大值,即:
3.如权利要求1所述的一种基于随机变分贝叶斯的基因调控网络拓扑结构识别方法,其特征在于,步骤3‑2)中基于随机优化的贝叶斯变分方法计算下界函数J最大化具体方法如下:关于ξ的下界函数为:
式中,ξ∈Θ,Θ‑ξ={Θ/ξ}, 表示在q(Θ)方面的目标期望,const.表示一个常数;给出J在q(ξ)方面有最大的值:其先验分布p(ξ)和后验分布q(ξ)具有共轭性质,属于指数族,即:T
lnq(ξ)=lns(ξ)+λt(ξ)‑ag(λ) (17)式中,λ是自然参数,t(ξ)是充分统计量,s(ξ)和ag(λ)分别是基本度量项和对数规则化器;对于q(ξ)的下界函数J改写为:式中:
J(λ)的自然梯度具有以下简单的形式:式中,自然梯度被定义为标准梯度与q(ξ)的费歇尔信息矩阵的逆的乘积,即:对识别数据进行子采样来得到自然梯度的噪声估计,如下:式中,t′从{1,2,…,N}均匀采样,|t′|是采样数据的长度;
基因调控网络所有的基因节点均为独立设置,则:k
基于随机梯度算法,λ将收敛到最优,它由下式更新:k
λ的收敛条件为:
k k 2
∑kρ=∞,∑k(ρ) <∞ (26)。
4.如权利要求1所述的一种基于随机变分贝叶斯的基因调控网络拓扑结构识别方法,其特征在于,步骤3‑3)中通过随机变分贝叶斯方法估计基因调控网络参数的具体方法为:
3‑2‑1)初始化后验分布q(Θ)、参数γ1,j和γ2,j、学习率ρ、参数k;
3‑2‑2)更新q(h)和q(l):q(hij)的先验分布是一个高斯分布,它的指数族模型为:式中,λ1是自然参数,s1 (hij) =1,充分统计量t1(hij)被定义为其中Vec(·)是一个由矩阵的行向量组成的向量;
下界函数J(λ1)的噪声自然梯度为:式中:
T T
式中,h‑j={h\hij},φj=[yi(t‑1),…,yi(t‑Z)] ,ψj=[uj(t‑1),…,uj(t‑Z)] ,Φ‑j={φ\φj};
自然参数λ1被更新为:
q(hij)是一个高斯分布,其均值和协方差为:式中,iVec被定义为将向量转换为矩阵的算子,它是Vec的逆运算,end表示向量的最后一个索引;
同样地,q(lij)的指数族模型由以下公式给出:下界函数J(λ2)的噪声自然梯度为:式中:
自然参数λ2被更新为:
其均值和协方差为:
3‑2‑3)更新q(αj)和q(βj):q(αj)为伽马分布,即:
T
式中,s3(αj)=1,t3(αj)=[αj,lnαj],下界函数J(λ3)的自然梯度为:q(αj)更新为:
αj的均值计算为:
同样地,q(βj)也为伽马分布,即 更新规则为:式中:
3‑2‑4)更新q(τ):
q(τ)的后验分布为伽马分布,即:q(τ)的下界函数用J(λ5)表示,其噪声自然梯度为:q(τ)更新为:
3‑2‑5)更新γ1,j和γ2,j:在每次最大化下界的迭代中,这两个值通过求解两个简单的优化问题来更新,γ1,j的下界函数为:每次迭代中,γ1,j和γ2,j更新为:k k k‑1
3‑2‑6)计算出下界J,若|J‑J |<∈,其中∈为阈值,则停止;否则跳转到步骤3‑2‑2);
3‑2‑7)所有的参数被估计为其期望值。
5.一种存储介质,其特征在于,所述存储介质存储有若干指令,所述指令适用于处理器进行加载,以执行权利要求1至4任意一项所述的方法。