人类常见疾病, 如神经退行性疾病、恶性肿瘤疾病等, 究其原因都是由于基因表达异常所导致的结果.而基因并不是孤立存在的, 基因间促进和抑制的调控关系形成了基因调控网络^[1].通过构建相对精准的基因调控网络来研究基因间的表达关系, 对人体机制和疾病治疗的探索有着重要意义.而如何构建一个更加准确有效的网络成为基因调控网络研究的首要难题^[2].目前, 基因调控网络建模的研究方法主要包括:关联模型^[3]、布尔网络模型^[4]、微分方程模型^[5]和贝叶斯网络模型^[6].近年来, 关联模型和贝叶斯网络模型广泛应用于基因调控网络的构建.在关联模型方面, 文献[7]基于自适应分块策略来估计互信息并构建基因调控网络.在贝叶斯网络模型方面, Frolova和Wilczynski^[8]提出了一种可以支持动态贝叶斯网络的BNFinder2软件来构建基因调控网络.

关联模型可以支持构建大规模的基因调控网络, 但构建出的网络不能描述调控方向^[9]; 而贝叶斯网络模型既可以描述调控关系又可以描述调控方向, 但由于其计算复杂度很高, 限制了网络构建的规模和效率^[10].基于此, 本文将两种模型的优点相结合, 提出了基于父节点筛选的贝叶斯网络模型方法来构建基因调控网络.首先, 对于每个节点, 利用皮尔逊相关系数计算节点间的关联程度, 筛选出与该节点具有强相关性的若干节点; 其次, 将筛选出的节点作为贝叶斯网络模型中结构学习时的候选父节点以缩小搜索空间, 并基于此构建基因调控网络.最后, 通过实验验证了该方法在准确率和计算效率上均优于传统的贝叶斯网络模型, 并且该方法可以支持大规模基因调控网络的构建.

1 基于父节点筛选的贝叶斯网络建模 1.1 总体框架

贝叶斯网络(Bayesian network, BN)^[10]借助有向无环图来刻画基因之间的依赖关系.对于任一变量X_i, 通常可以找到一个与X_i都不独立的最小子集Parent(X_i)⊆{X₁, X₂, …, X_i－1}, 使得P(X_i|X₁, X₂, …, X_i－1)=P(X_i|Parent(X_i)).因此, 当网络变量元组〈X₁, X_i, …, X_n〉赋予具体数据值〈x₁, x_i, …, x_n〉时, 贝叶斯网络的联合概率分布为

(1)

建立贝叶斯网络通常需要2个步骤:结构学习和参数学习.结构学习的方法是针对每个节点, 遍历并通过评分函数评价所有可能的结构, 进而找出最好的结构作为该节点的父节点集; 该方法的寻优策略可以保证所构建出的网络结构的精确性, 然而, 其复杂度也极高.因此, 利用传统贝叶斯网络模型的结构学习构建基因调控网络, 不仅效率非常低且网络规模也有限.本文提出一种基于父节点筛选的贝叶斯网络建模方法, 该方法可以缩小结构学习的搜索空间, 同时, 充分利用结构学习搜索策略的优势保证所构建的基因调控网络的精确性.网络总体框架如图 1所示.

由图 1可以看到, 基于父节点筛选的贝叶斯网络建模方法在贝叶斯网络结构学习前, 先利用皮尔逊相关系数方法按照节点间的关联程度, 针对每个节点筛选出若干个候选父节点, 然后, 将这些父节点作为结构学习时的搜索空间, 但不改变结构学习的搜索策略, 从而在结构学习过程中准确、高效地得出父节点集; 经过结构学习后即可得到贝叶斯网络模型所构建出的基因调控网络.最后, 通过参数学习得到节点间边的权重, 进而得出既有调控方向又有概率值的基因调控网络.

1.2 父节点的筛选方法

父节点筛选方法采用皮尔逊相关系数对节点间的关系进行描述.皮尔逊相关系数衡量了两个变量X_i和Y_i的相似度, 可由式(2)定义:

(2)

式中X和Y分别表示两个变量X_i和Y_i的平均值.

根据式(2)得到两个节点的相关程度后, 将具有一定相关程度的父节点作为候选父节点.基于父节点筛选的贝叶斯网络建模过程如算法1所示.

算法 1  基于父节点筛选的贝叶斯网络建模方法

输入:基因表达数据集合{X_i}_i=1ⁿ, 父节点筛选比例α, 父节点集个数m.

输出:基因调控网络矩阵G.

算法描述:

1)   父节点筛选个数λ=α×n;

2)   for i=1 to n do

3)     {Y_j}_j=1^n－1={X₁, …, X_i－1, X_i+1, …, X_n};

4)     for j=1 to n-1 do

5)       计算X_i和Y_j的平均值X, Y;

6)       根据式(2)计算皮尔逊相关系数r_j;

7)       将[r_j, Y_j]添加到父节点集合P_father;

8)     对P_father∈{Pa_k}_k=1^n－1根据r_j由大到小排序;

9)     for k=1 to λ do

10)       将Pa_k中的Y_j存入X_i的候选父节点集合P_i中;

11)       计算X_i父节点集为空的BDE分数emptyscore;

12)      设置最优父节点集optimal_set =[X_i, []];

13)      设置最优父节点分数optimal_score=

           emptyscore;

14)      for p=1 to m do

15)       计算可能的父节点组合数

16)       for q=1 to pc do

17)         在候选父节点集合P_i中计算父节

              点集组合Pa_q的BDE分数score;

18)         if score>optimal_score

19)           optimal_set = [X_i, [Pa_q]];

20)           optimal_score = score

21)       将[optimal_set, optimal_score]保存到G_i;

22)       G_i添加到G;

23)   输出G.

算法分两部分:候选父节点筛选过程(第1~10行)和贝叶斯调控网络结构学习(第11~23行).在父节点筛选过程中, 先根据输入的父节点筛选比例计算父节点筛选个数λ(第1行); 其次, 取基因表达数据集合{X_i}_i=1ⁿ中的一个节点X_i, 将该节点以外其余所有节点设置为{Y_j}_j=1^n－1(第2~3行); 再次, 计算节点X_i与其余所有节点的皮尔逊相关系数, 并将相关系数及其表达数据存储到父节点集合P_father中(第4~7行); 最后, 根据皮尔逊相关系数的取值对P_father中的行从大到小进行排序, 并取前λ行中的数据Y_j作为X_i节点所对应的候选父节点集合(第8~10行); 在贝叶斯网络结构学习过程中, 根据X_i节点对应的候选父节点及其基因表达数据, 先计算父节点集为空时的BDE分数, 并且将该结构及其对应的BDE分数作为X_i节点的父节点集的最优组合和最优分数初始值(第11~13行); 然后, 根据输入的父节点集个数和父节点筛选个数λ, 计算从1个父节点作为父节点集到m个父节点作为父节点集的所有组合数, 即为可能的父节点组合数(第14~15行); 接着, 依次在父节点筛选时所筛选的集合P_i中计算所有可能父节点的BDE分数并与最优分数进行比较, 并根据计算结果更新最优组合和最优系数，直至得到X_i节点最终的最优组合和最优分数并保存到G_i中(第16~22行); 最后, 重复上述所有步骤, 计算出所有节点的结构及其BDE分数并输出基因调控网络矩阵G.

2 实验结果及分析 2.1 实验设置

实验选用的数据为GeneNetWeaver上获取的大肠杆菌基因调控网络.选取了其中20个基因的子网络, 对父节点筛选的贝叶斯网络模型(PS-BN)和传统的贝叶斯网络模型(BN)的性能进行评价; 还分别选取了包含100, 200, 300, 400和500个基因的子网络对PS-BN的大规模基因调控网络构建效率进行评价.

在对PS-BN和BN进行性能评价时, 选用了准确率、精确率、召回率、F值4项评估指标, 另外还对两种方法的运行时间进行了比较, 验证了PS-BN的效率.

准确率描述了构建基因调控网络时判断为有边或无边的总体的准确率, 可由式(3)计算:

(3)

精确率描述了判断为有边的结果中, 预测正确的概率, 可由式(4)计算:

(4)

召回率描述的是在所有标签为有边的金标准中, 实验结果也判断为有边的概率, 可由式(5)计算:

(5)

F值是一个综合评价指标, 表示的是精确率和召回率的调和平均评估指标, 可由式(6)计算:

(6)

上述评估指标中，

TP表示调控网络真实为有边, 构建结果也为有边;

TN表示真实为无边, 结果也为无边;

FP表示真实为无边, 结果却为有边;

FN表示真实为有边, 结果却为无边.

2.2 实验结果

将父节点筛选比例分别设置为总节点个数的10 %, 20 %, 30 %, 40 %和50 %, BN和PS-BN两种方法的评估指标与运行时间的实验结果对比如表 1所示.可以看出, 对于PS-BN方法, 父节点筛选比例高的评估指标优于筛选比例较低的; 当筛选比例大于或等于30 %时, 各项评估指标值相同.与BN方法相比, PS-BN方法筛选比例大于或等于20 %时的实验结果均优于BN方法.在运行时间方面, 各个比例的PS-BN方法的运行时间均远远少于BN方法.因此, PS-BN方法在保证评估指标优于BN方法的前提下, 计算效率得到了显著提高.

在准确率方面, BN方法与筛选比例为10 %的PS-BN方法的准确率相同.因为, 在经过父节点筛选后, 贝叶斯网络搜索空间中的绝大部分冗余信息都被过滤掉, 因而假阳边非常少; 但由于筛选比例过小而导致部分真阳边也被过滤掉.因此, 准确率是在牺牲真阳边数量的前提下得以保证的.筛选比例为20 %和30 %时, 准确率逐渐升高, 且均高于BN方法.随着筛选比例的增高, 搜索空间中对应的真阳边越来越多, 使得所构建的基因调控网络真阳边也增加，而大部分冗余信息仍会被过滤掉；因此在保证真阳边数量的前提下, 假阳边的比例少于BN方法, 从而提升其准确率.而继续增加父节点筛选比例后，由于有效节点已经被筛选进搜索空间, 而增加少部分冗余信息并没有使假阳边增多，因此, 其准确率与30 %时相比不变.

精确率和召回率的评估考虑的均是真阳边的预测情况.当PS-BN方法的父节点筛选比例设置为10 %时, 精确率和召回率均低于BN方法, 说明该比例所得到的真阳边数量少于BN方法.当筛选比例增加至20 %及以上时, 精确率和召回率得到提升并高于BN方法, 也就是说, 筛选比例增加, 真阳边数量也随之增多, 并且超过了BN方法所得到的真阳边.当筛选比例高于30 %时, 真阳边数量不再增加, 并且网络结构没有变化, 因此, 精确率和召回率保持不变.

评估指标F值是平衡精确率和召回率的一项综合指标, 因此该指标与精确率和召回率的相关程度较高, F值的变化趋势与精确率和召回率的变化趋势一致.筛选比例为10 %的PS-BN方法在精确率和召回率上均低于BN方法, 因此其F值也较低.当筛选比例为20 %和30 %时, 两项评估指标均高于BN方法, 因而其F值也所有提高.当筛选比例为40 %和50 %时, 两项评估指标均与30 %相同, 故其F值也与30 %相同.

在运行时间方面, PS-BN方法均优于BN方法.因为父节点经过筛选后, 贝叶斯网络的搜索空间大幅缩小, 在搜索空间中进行父节点集遍历搜索的结构学习时, 搜索的节点数量和各节点的组合数目都大量减少, 因此, PS-BN方法的运行时间大大缩短.随着筛选比例的增加, 贝叶斯网络模型的搜索空间逐渐增大, 不但导致需要搜索的节点数量增加, 还导致需要进行评分计算的组合数目大量增多, 使结构学习耗费的时间呈指数增长.因此, 随着筛选比例的增加, 运行时间也呈指数增长.

由表 1可知, 在评估指标值较高的筛选比例中, 比例为30 %时运行时间最短.因此, 在对大规模基因调控网络的构建进行效率评估时, 选用的筛选比例为30 %.大规模基因调控网络构建运行时间如图 2所示，其中带标记的实线为运行时间, 虚线表示运行时间增加的趋势线.由图 2可知, 随着基因数量的增加, 运行时间呈幂增长趋势.这是由于基因数量增加, 同样的筛选比例所筛选出的基因数量也随之增加; 但由于筛选比例较小, 基因数量增加有限, 因此, 在该搜索空间进行结构学习时, 运行时间呈幂增长趋势.

3 结语

为解决贝叶斯网络模型在构建基因调控网络时效率低且构建规模有限的问题, 本文提出了基于父节点筛选的贝叶斯网络(PS-BN)建模方法.PS-BN方法充分利用传统贝叶斯网络建模方法在结构学习中的搜索策略, 同时, 通过筛选父节点缩小了搜索空间且去除了部分冗余信息, 从而在极大提高效率的同时, 准确率等4项评估指标也均有所提升.实验证明了上述结论.