CN103646159B - 一种基于约束性布尔网络的最大评分预测方法 - Google Patents

Abstract

本发明涉及一种基于约束性布尔网络的最大评分预测方法，包括有以下步骤：第一步、计算相对互信息确定候选预测基因集，第二步、采用最大评分预测方法从第一步中得到的预测基因集中选择作为目标基因的父基因，本发明具有本发明适用于依据小样本数据预测多变量之间的关系，本发明对噪声的鲁棒性更强，更适用于真实环境下的生物数据的预测，发明预测的网络结构更加准确、详细，体现在预测调控关系的正确数目及调控关系的方向性和正负调控关系。

Description

一种基于约束性布尔网络的最大评分预测方法

技术领域

本发明涉及一种基因调控网络预测的方法，特别涉及一种利用约束性布尔网络特性设计基因调控网络的最大评分预测方法。

背景技术

系统生物学研究的一个重要目标是描述调控特定细胞行为和过程的分子机制。描述基因调控网络的模型有很多，例如：贝叶斯网络和动态贝叶斯网络提供了一种能阐明基因之间的依赖关系的模型；而布尔网络和概率布尔网络提供了一种通过网络稳态行为研究系统功能的方法；微分方程则是一种连续性模型，它能描述基因之间详细的生化关系。这些模型都被统一的用来研究生物现象(细胞周期)和疾病(癌症)。所以，揭示这些网络的结构是一个非常关键性的问题。

用来阐明潜在的转录调控网络结构的数据有多种类型，包括利用DNA芯片技术或其他高通量技术的基因组转录谱；用染色质免疫沉淀技术得到的转录因子结合位点数据和其他来源的分子间相互作用的数据。从已知的观察数据推理基因调控网络是一个逆向工程。可是这个过程是一个多对一的病态问题，因为不仅仅只有一个结构满足给定数据，尤其是当面对小样本数据时，系统中的变量数(基因)远远大于样本数。

现有的基因调控网络预测方法有很多，例如：Reveal方法、Cod、Best-fit、MDL等等，但是它们都有各自方法的不足，如Reveal可适性差，在小样本数据中很小的噪声会产生比较大的误差；Cod、Best-fit推荐的符合条件的模型空间太大，往往存在选择预测模型上的偏差；MDL则参数较多，难已确定。此外，这些方法预测出来的基因调控网络中基因之间只有简单的调控与被调控关系。而约束性布尔网络由于其定义的特殊性，它能描述基因之间更加详细的正负调控关系，而且它的函数空间比布尔网络小很多(当K＝2时，有个可能的约束性布尔函数，当K＝3，它仅仅只有个可能的约束性布尔函数，K为基因入度)，在预测推理时具有比较多的优势。而现有的基于约束性布尔网络预测方法——三条规则，它对数据质量要求很高，也只适用于小样本数据的推理，而真实环境下的生物数据包含较多噪声，所以它一般仅被用作预测的前期处理。

发明内容

本发明的目的在于提供一种基于约束性布尔网络可适用于真实环境下的基因调控网络预测方法。

本发明的目的是通过如下方案来实现的：一种基于约束性布尔网络的最大评分预测方法，包括有以下步骤：

第一步、计算相对互信息确定候选预测基因集，其包括以下步骤：

一、根据以下公式计算两两变量之间的相对互信息

确定一个相对互信息矩阵，

式中表示预测基因x_j从时刻1到t的表达值；

表示目标基因x_i从时刻2到t+1的表达值；

是预测基因x_j与目标基因x_i延迟一个单位时间的条件互信息；

H(x)表示变量x的熵；

θ_ij表示预测基因x_j与目标基因x_i的相对互信息；

二、计算相对互信息矩阵中每一行的平均值，将之作为阈值δ_i，根据以下公式计算：

${\mathrm{\δ}}_i=\frac{1}{n}\underset{j=1}{\overset{n}{\Σ}}{\mathrm{\θ}}_{ij}$

式中n表示相对互信息矩阵中每一行的个数，

三、定义一个M矩阵，确定目标基因各自的候选预测基因集，根据以下方程计算：

$M_{ij}=\left\{\begin{array}{cc}0& {\mathrm{\&theta;}}_{ij}<{\mathrm{\&delta;}}_i\\ 1& {\mathrm{\&theta;}}_{ij}\&GreaterEqual;{\mathrm{\&delta;}}_i\end{array}\right.$

M_ij中i行第j列为1，则选择第j个基因作为i基因的候选预测基因集中的一个；第二步、采用最大评分预测方法从第一步中得到的预测基因集中选择作为目标基因的父基因，其包括以下步骤：

①根据预测基因时间点t的值x_j(t)与目标基因时间点t和时间点t+1的值x_i(t)、x_i(t+1)，确定调控关系表R(m-1)×p(x_i)，其中，p(x_i)表示目标基因x_i的预测基因的个数，m表示样本点/时间点个数；

调控关系表R(m-1)×p(x_i)每一行的元素将由x_j(t)、x_i(t)和x_i(t+1)按照调控关系a_ij的约束关系进行求解，具体过程如下:当x_i(t)＝0，x_i(t+1)＝0，对a_ij的约束为当x_i(t)＝0，x_i(t+1)＝1，对a_ij的约束为当x_i(t)＝1，x_i(t+1)＝0，对a_ij的约束为当x_i(t)＝1，x_i(t+1)＝1，对a_ij的约束为

并根据调控关系a_ij进行编码，其中：a_ij＝-1表示预测基因对目标基因具有负调控；a_ij＝1表示预测基因对目标基因具有正调控；a_ij＝2表示未完全确定，表示预测基因对目标基因的调控关系不确定但存在一定的约束；a_ij＝0表示预测基因对目标基因调控关系不确定且不存在任何约束；

②定义预测基因x_j对目标基因x_i的调控关系aⁱ _j的确定性为：

$d_{ij}=\left|N_{ij}^{-1}-N_{ij}^1\right|,$

式中 分别表示a_ij＝-1，a_ij＝1的个数，并通过和计算；

③按照步骤②得到的确定性d_ij从大到小依次确定预测基因x_j对目标基因x_i最可能的调控关系a_ij的值，如果那么a_ij＝-1，否则a_ij＝1，将已确定的a_ij代入所述的调控关系表再确定下一个a_ik，k≠j；

④统计错误个数：

采用计算关系冲突错误；

根据步骤①中获得的编码进行统计计算1的个数无输入错误：

⑤定义评分准则：

$s_i=\underset{j}{\&Sigma;}(d_{ij}-{\mathrm{\&epsiv;}}_{ij}^{-1,1}-{\mathrm{\&epsiv;}}_i^{null})/P\left(x_i\right)$

式中s_i表示预测结果，从预测基因集合中选择使得s_i最大的预测基因组合作为目标基因x_i的父基因，且基因入度K≤3。

本发明有以下优点：

1.本发明适用于依据小样本数据预测多变量之间的关系。

2.本发明对噪声的鲁棒性更强，更适用于真实环境下的生物数据的预测。

3.本发明预测的网络结构更加准确、详细，体现在预测调控关系的正确数目及调控关系的方向性和正负调控关系。

附图说明

附图1为真实的芽殖酵母细胞的关键性基因组成的基因调控网络结构图。

附图2为用三条规则预测得到的基因调控网络结构图。

附图3为用最大评分预测方法预测得到的基因调控网络结构图。

附图4为目标基因的状态变化对预测基因权值的约束关系表。

附图5为调控关系的编码方式。

附图6为无输入状态下错误判断表。

附图7为芽殖酵母细胞基因调控网络的一个时序数据。

附图8为用三条规则和最大评分预测方法在不同噪声数据中的预测结果。

具体实施方式

本发明一种基于约束性布尔网络的最大评分预测方法，包括有以下步骤：

第一步、计算相对互信息确定候选预测基因集，

①定义一个M矩阵，根据以下公式计算两两变量之间的相对互信息

${\mathrm{\&theta;}}_{ij}=\frac{MI(x_j^{t+1},x_i^{t+1})}{min\left(H\right(x_j^t),H(x_i^{t+1}\left)\right)}$

式中表示预测基因x_j从时刻1到t的表达值；

表示目标基因x_i从时刻2到t+1的表达值；

是预测基因x_j与目标基因x_i延迟一个单位时间的条件互信息；

H(x)表示变量x的熵；

θ_ij表示预测基因x_j与目标基因x_i的相对互信息。

②计算每一行的平均值，将之作为阈值δ_i，根据以下公式计算：

${\mathrm{\&delta;}}_i=\underset{j=1}{\overset{n}{\&Sigma;}}{\mathrm{\&theta;}}_{ij}/\underset{j=1}{\overset{n}{\&Sigma;}}j$

式中n表示相对互信息矩阵中每一行的个数。

③确定目标基因各自的候选预测基因集，根据以下方程计算：

$M_{ij}=\left\{\begin{array}{cc}0& {\mathrm{\&theta;}}_{ij}<{\mathrm{\&delta;}}_i\\ 1& {\mathrm{\&theta;}}_{ij}\&GreaterEqual;{\mathrm{\&delta;}}_i\end{array}\right.$

M_ij中i行第j列为1，则选择第j个基因作为i基因的候选预测基因集中的一个。

第二步、采用最大评分预测方法从第一步中得到的预测基因集中选择作为目标基因的父基因，

①构建预测基因与目标基因每两相邻时间点之间确定的调控关系表R_(m-1)×n(m表示样本点个数，n表示预测基因个数)，将调控关系表R_(m-1)×n中每一行的元素由每两相邻时间点之间按照调控关系a_ij的约束关系进行求解，并根据调控关系a_ij进行编码；步骤①中所述调控关系a_ij的约束关系包括 和详见附图4，步骤①中所述的调控关系a_ij包括预测基因对目标基因具有负调控、预测基因对目标基因具有正调控、预测基因对目标基因的调控关系不确定但存在一定的约束以及预测基因对目标基因的调控关系不确定，且不存在任何约束，详见附图5。

调控关系表R_(m-1)×n中每一行的元素由每两相邻时间点之间根据附图4中相应情况所列不等式解得并按附图5进行编码赋值。

附图5中-1和1表示调控关系完全确定，其中-1表示预测基因对目标基因具有负调控，1表示预测基因对目标基因具有正调控；2表示未完全确定，表示预测基因对目标基因的调控关系不确定但存在一定的约束，例如a_i1+a_i2+a_i3≥0，那么当a_i1为1时，a_i2和a_i3既能为-1也能为1，但不能同时为-1；0表示预测基因对目标基因的调控关系不确定，且不存在任何约束，例如预测基因状态为0。

②定义预测基因x_j对目标基因x_i的调控关系a_ij的确定性为：

$d_{ij}=\left|N_{ij}^{-1}-N_{ij}^1\right|$

d_ij越大，确定性越大。式中 分别表示a_ij＝-1，a_ij＝1的个数，并通过和计算。

③按照确定性d_ij从大到小依次确定预测基因x_j对目标基因x_i最可能的调控关系a_ij的值，如果那么a_ij＝-1，否则a_ij＝1，应当注意将已确定的a_ij代入前面的各个不等式更新R表，再确定下一个a_ik(k≠j)。

④统计错误个数：错误分为两部分，一部分为关系冲突错误，另一部分为无输入错误。

关系冲突错误：

无输入错误：根据附图6进行统计计算1的个数。其中当无输入且目标基因对自身的调控关系为负调控时，则为self-degradation，否则无输入时的其他情况为noself-degradation。

⑤定义评分准则：

$s_i=\underset{j}{\&Sigma;}(d_{ij}-{\mathrm{\&epsiv;}}_{ij}^{-1,1}-{\mathrm{\&epsiv;}}_i^{null})/\left|P\left(x_i\right)\right|$

s_i表示预测结果，当s_i越大，说明预测结果越好。从预测基因集合中选择使得s_i最大的预测基因组合(K≤3)作为目标基因x_i的父基因。

本发明中根据相对互信息的计算排除掉与目标基因极可能不相关的基因，缩小搜索空间，提高计算效率，根据约束性布尔网络的定义，从目标基因状态的变化特点出发反推预测基因与目标基因之间可能的调控关系，结合数理统计知识，确定调控关系、确定性大小及误差个数，综合确定性及误差个数定义评分机制，选择最大评分预测基因组合作为最终的父基因，它能发现更多基因之间的调控关系，同时消除了噪声的连带效应，对噪声的适应性更强，更适合存在噪声的真实环境中的生物数据预测推理。

本发明中对目标基因的预测是相互独立的，同时它可以限定预测基因的个数，使得在小样本数据中能对多变量进行预测，这恰好符合真实生物数据的特点。

下面通过示例以进一步说明本发明，该示例不对本发明构成任何限制：示例一：有一芽殖酵母细胞的基因调控网络如附图1。附图7为它的一个时序表达数据。分别用三条规则和最大评分预测方法预测图的基因调控网络，结果如附图2和附图3所示。

附图2和附图3中三条规则的图中粗实线表示预测正确的调控关系，细虚线表示预测出的未完全确定且不包含正确的调控关系，粗虚线则表示预测出的未完全确定且包含正确的调控关系。这里的未完全确定指的是两基因之间的调控关系可能不存在或者是正负调控关系中的一种。最大评分预测图中的粗实线表示预测正确的调控关系，细虚线表示预测错误的调控关系。

示例二：对示例一中的时序数据分别添加0.5％和1％的随机噪声200次。用得到的噪声数据分别用三条规则和最大评分预测方法预测基因调控网络。然后用召回率和精确度来评价结果的优劣，结果如附图8。

由示例一可见，最大评分预测方法预测的芽殖酵母细胞的基因调控网络中正确的调控关系数比三条规则多，而且它的网络连通度比三条规则好，更符合真实网络的结构，并且结果中没有不确定性的成分存在，更具预测指导意义。由示例二可见，从平均召回率来看，在两种噪声的情况下，最大评分预测都比三条规则大很多。这说明最大评分预测比三条规则能够预测到更多正确的调控关系。从平均精确度来看，在两种噪声的情况下，最大评分预测都要小于三条规则，但是相差不大，尤其是当噪声等于1％时，两者仅相差0.0092。综合两者来看，最大评分预测明显要优于三条规则，因为平均召回率的提高幅度明显比平均精确度的下降幅度大。

由此可见最大评分预测方法不仅对小样本数据的多变量关系预测适用，并且在多噪声的情况下，能够获得比较好的预测效果，更适合用真实环境下的生物数据来预测基因调控网络。

虽然已结合附图对本发明实施例做了详细的描述，但是本领域的技术人员应该明白以上的实施例仅是本发明的实施方式，而并非对本发明的保护范围的限定，任何基于本发明的发明所做的改进都理所当然在本发明保护范围内。

Claims (1)

1.一种基于约束性布尔网络的最大评分预测方法，其特征在于：包括有以下步骤：

第一步、计算相对互信息确定候选预测基因集，其包括以下步骤：

一、根据以下公式计算两两变量之间的相对互信息

确定一个相对互信息矩阵，

式中表示预测基因x_j从时刻1到t的表达值；

表示目标基因x_i从时刻2到t+1的表达值；

是预测基因x_j与目标基因x_i延迟一个单位时间的条件互信息；

H(x)表示变量x的熵；

θ_ij表示预测基因x_j与目标基因x_i的相对互信息；

二、计算相对互信息矩阵中每一行的平均值，将之作为阈值δ_i，根据以下公式计算：

${\mathrm{\&delta;}}_i=\frac{1}{n}\underset{j=1}{\overset{n}{\&Sigma;}}{\mathrm{\&theta;}}_{ij}$

式中n表示相对互信息矩阵中每一行的个数，

三、定义一个M矩阵，确定目标基因各自的候选预测基因集，根据以下方程计算：

$M_{ij}=\left\{\begin{array}{cc}0& {\mathrm{\&theta;}}_{ij}<{\mathrm{\&delta;}}_i\\ 1& {\mathrm{\&theta;}}_{ij}\&GreaterEqual;{\mathrm{\&delta;}}_i\end{array}\right.$

M_ij中i行第j列为1，则选择第j个基因作为i基因的候选预测基因集中的一个；

第二步、采用最大评分预测方法从第一步中得到的预测基因集中选择作为目标基因的父基因，其包括以下步骤：

①根据预测基因时间点t的值x_j(t)与目标基因时间点t和时间点t+1的值x_i(t)、x_i(t+1)，确定调控关系表R(m-1)×p(x_i)，其中，p(x_i)表示目标基因x_i的预测基因的个数，m表示样本点/时间点个数；

调控关系表R(m-1)×p(x_i)每一行的元素将由x_j(t)、x_i(t)和x_i(t+1)按照调控关系a_ij的约束关系进行求解，具体过程如下:当x_i(t)＝0，x_i(t+1)＝0，对a_ij的约束为当x_i(t)＝0，x_i(t+1)＝1，对a_ij的约束为当x_i(t)＝1，x_i(t+1)＝0，对a_ij的约束为当x_i(t)＝1，x_i(t+1)＝1，对a_ij的约束为

并根据调控关系a_ij进行编码，其中：a_ij＝-1表示预测基因对目标基因具有负调控；a_ij＝1表示预测基因对目标基因具有正调控；a_ij＝2表示未完全确定，表示预测基因对目标基因的调控关系不确定但存在一定的约束；a_ij＝0表示预测基因对目标基因调控关系不确定且不存在任何约束；

②定义预测基因x_j对目标基因x_i的调控关系a_ij的确定性为：

$d_{ij}=|N_{ij}^{-1}-N_{ij}^1|,$

式中分别表示a_ij＝-1，a_ij＝1的个数，并通过和计算；

③按照步骤②得到的确定性d_ij从大到小依次确定预测基因x_j对目标基因x_i最可能的调控关系a_ij的值，如果那么a_ij＝-1，否则a_ij＝1，将已确定的a_ij代入所述的调控关系表再确定下一个a_ik，k≠j；

④统计错误个数：

采用计算关系冲突错误；

根据步骤①中获得的编码进行统计计算1的个数无输入错误：

⑤定义评分准则：

$s_i=\underset{j}{\&Sigma;}(d_{ij}-{\mathrm{\&epsiv;}}_{ij}^{-1,1}-{\mathrm{\&epsiv;}}_i^{null})/P\left(x_i\right)$

式中s_i表示预测结果，从预测基因集合中选择使得s_i最大的预测基因组合作为目标基因x_i的父基因，且基因入度K≤3。