一种基于高阶动态贝叶斯的基因调控网络构建方法

1.本发明涉及医学信息学技术领域，具体涉及一种基于高阶动态贝叶斯的基因调控网络构建方法。
技术背景
2.基因在生物遗传进化发展历程中担任着不可或缺的角色，是人们认识和探索生命奥秘必不可少的重要因素。针对基因之间遗传表达关系的研究，基因调控网络提供了一个重要且有效的途径。基因调控网络是由基因与基因之间相互制约的关系构成的，各个基因之间相互影响相互制约，形成一个复杂的网络。通过基因表达数据去寻找不同基因在不同时间延迟(不同调控基因与目标基因表达时刻的时间差)下复杂的调控关系来准确地构建基因调控网络，将有利于更系统地分析出生物遗传进化过程中各个基因的功能以及它们在遗传进化过程中起到的作用。
3.由于基因调控网络十分复杂，不仅表现在所需构建的网络规模较为庞大，而且基因之间的调控关系还会受时间因素的制约，利用传统的重建算法对基因调控网络进行重建时会存在不能准确地描述时间因素对基因之间调控关系的影响和在网络结构学习过程中时间复杂度过高的问题。如何去构建能够更准确描述基因之间调控关系并且构建效率更高的基因调控网络成为一种挑战。
4.申请公布号为cn111223523a的中国专利申请《基于多时滞因果熵的基因调控网络构建方法及系统》将输入的时间序列基因表达数据划分成不同时滞下的时间窗口；对t个时间片的基因表达数据，分别构建tτ个时间窗口下的基因表达矩阵，对每一对基因，计算t时间窗口下的目标基因与之前tτ个时间窗口下基因的多时滞转移熵，得到一个多时滞下的基因相关性矩阵，矩阵的元素代表基因之间存在边的概率，并对矩阵通过k-means将边聚类成两类，过滤掉低概率的边簇，对剩下的每条边计算在条件基因下的多时滞条件转移熵，对最大因果熵小于阈值的间接调控的边过滤掉，获得最终的网络结构。该专利申请提供的技术方案在计算t时间窗口下的目标基因与之前tτ个时间窗口下基因之间的相关关系时，是将前tτ个时间窗口下所有基因都计算了一遍，使得网络结构学习的搜索空间呈指数级增长，虽然为目标基因寻找了不同时间窗口下的相关基因，提高了网络构建的准确率，但是网络构建的时间复杂度非常高、导致网络构建时间长。
5.申请公布号为cn110675912a的中国专利申请《一种基于结构预测的基因调控网络构建方法》首先计算系数矩阵，通过计算基因之间的pearson系数，互信息及最大互信息来确定基因之间的相关性，作为筛选潜在父节点集的依据；然后进行结构预测，将获得的基因之间的系数矩阵作为判定基因潜在父节点集的依据，为每个基因选取潜在父节点集；最后进行基因调控网络的结构学习和参数学习。该专利申请提供的技术方案在进行结构学习之前首先为目标基因筛选父节点集，这样做可以大大缩小搜索空间，但是该方法所构建的基因调控网络是基于非时序数据的，不能表明时间因素在基因相互调控关系中的作用，使得所构建的基因调控网络准确率较低。


技术实现要素：

6.针对现有技术存在的不足，本发明提供一种基于高阶动态贝叶斯的基因调控网络构建方法，旨在解决构建基因调控网络过程中出现的准确率低、网络构建时间长的问题。
7.为达到此目的，本发明的技术方案是：
8.一种基于高阶动态贝叶斯的基因调控网络构建方法，该方法包括如下步骤：
9.步骤1：根据时序基因表达数据，为每个目标基因在不同时间延迟下筛选潜在的调控因子形成每个目标基因的潜在调控因子集；
10.步骤2：将每个目标基因的潜在调控因子集作为在利用高阶动态贝叶斯模型为每个目标基因寻找最优父节点集的搜索空间，构建局部最优的转移网络；
11.步骤3：将各个目标基因对应的局部最优的转移网络进行合并形成全局最优的转移网络，并将其输出作为最终的基因调控网络。
12.进一步地，根据所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述步骤1包括如下步骤：
13.步骤1.1：将时序基因表达数据处理成矩阵形式，其中所述矩阵的每一列表示一个基因；所述矩阵的每一行表示基因在一个时间片的样本；
14.步骤1.2：根据为目标基因所要寻找的不同时间延迟下的潜在调控因子，将目标基因和与目标基因具有不同时间延迟的其他基因在时序基因表达数据的选取上依据数据对齐原则进行数据对齐；
15.步骤1.3：根据数据对齐后的时序基因表达数据，结合互信息方法和预测最小描述长度方法为目标基因进行第一次筛选潜在调控因子；
16.步骤1.4：同样根据数据对齐后的时序基因表达数据，结合皮尔逊相关系数方法和基于假设检验的断点检测方法为目标基因进行第二次筛选潜在调控因子；
17.步骤1.5：将步骤1.3和1.4分别获得的目标基因的潜在调控因子集取并集，得到目标基因最终的潜在调控因子集。
18.进一步地，根据所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述数据对齐，就是在计算目标基因i和与其有r个时间延迟的基因j之间的关联程度时，从时序基因表达数据中选取所使用的数据时，将目标基因i所对应的表达数据与要计算的与目标基因i具有r个时间延迟的基因j所对应的表达数据进行匹配。
19.进一步地，根据所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述为目标基因进行第一次筛选潜在调控因子的方法为：首先根据对齐好的数据计算目标基因i和与目标基因i具有不同时间延迟的其他基因j之间的互信息值，形成一个互信息矩阵；然后分别以该矩阵中的每一个互信息值分别作为阈值进行筛选，得到一个模型矩阵；再然后针对这个模型矩阵利用预测最小描述长度方法计算数据长度，选取具有最小长度的模型矩阵所对应的互信息值作为阈值对互信息矩阵进行筛选，将矩阵中大于阈值的互信息值保留，小于阈值的互信息值设为0；最后选择与目标基因之间互信息值大于0的基因作为目标基因的潜在调控因子。
20.进一步地，根据所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述为目标基因进行第二次筛选潜在调控因子的方法为：首先根据对齐好的数据计算基因i和与目标基因i具有不同时间延迟的其他基因j之间的皮尔逊相关系数值，形成一个皮尔逊相关系
数矩阵；然后，假设计算的是d阶转移网络，此时分别对比目标基因i与1至(d-1)时间延迟下的基因j的皮尔逊相关系数值，选最大值作为目标基因i与基因j的皮尔逊相关系数值，并对每个基因遍历该过程；再然后将目标基因i与其他所有基因的皮尔逊系数值进行排序，基于假设检验原理进行断点检测，所寻找到的断点就是阈值；又然后利用阈值对皮尔逊相关系数矩阵进行筛选，将矩阵中大于阈值的皮尔逊相关系数值保留，小于阈值的皮尔逊相关系数值设为0；最后选择与目标基因之间皮尔逊相关系数值大于0的基因作为目标基因的潜在调控因子。
21.进一步地，根据所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述步骤2包括如下步骤：
22.步骤2.1：依次遍历各个目标基因的潜在调控因子集，将各个目标基因的潜在调控因子集作为接下来进行结构学习的搜索空间；
23.步骤2.2：根据各个目标基因的搜索空间，利用高阶动态贝叶斯模型为每个目标基因在搜索空间中选择不同的父节点组合并进行评分，且从中选出得分最高的父节点组合作为目标基因的父节点集，得到局部最优的转移网络。
24.进一步地，根据上述任一项所述的基于高阶动态贝叶斯的基因调控网络构建方法，所述时序基因表达数据是从公开数据库dream challenge中下载的dream4时序基因表达数据。
25.总体而言，通过本发明所构思的以上技术方案较现有技术具有以下有益效果：本方法在进行基因调控网络构建时通过使用时序基因表达数据，在不同时间延迟下为目标基因寻找父节点，提高了网络构建的准确率；同时为了缩小网络构建的时间，使用基于互信息、皮尔逊相关系数、预测最小描述长度和断点检测的方法，在结构学习前先进行了目标基因的潜在调控因子筛选，缩小了结构学习时的搜索空间，大大减少了网络结构学习的时间，在一定程度上减少了基因调控网络的构建时间，提升了计算性能，可以更加快速，准确地构建大规模基因调控网络。
附图说明
26.图1为本实施方式基于高阶动态贝叶斯的基因调控网络构建方法的流程示意图；
27.图2为本发明实施例的目标基因筛选潜在调控因子集的过程示意图；
28.图3为本发明实施例的目标基因和与目标基因具有不同时间延迟的其他基因在时序基因表达数据的选取上进行数据对齐的过程示意图；
29.图4为本发明实施例的基于5个基因构建2阶动态贝叶斯基因调控网络的转移网络的过程示意图。
具体实施方式
30.为了使本发明的目的、技术方案及优势更加清晰，下面结合附图和具体实施例对本发明做进一步详细说明。此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
31.本发明的核心思路：通过利用时序基因表达数据去构建基于时间变量的高阶动态贝叶斯基因调控网络，在不同的时间延迟下为目标基因寻找父节点，以提高网络结构构建
的准确率。由于高阶网络结构学习的搜索空间非常大，导致结构学习的时间复杂度特别高，为了加快学习速度，提出了在结构学习之前进行潜在调控因子筛选的办法，首先筛选出与目标基因关联度较高的基因作为目标基因的潜在调控因子集，删除与目标基因关联程度较低的基因，以达到缩小搜索空间的目的，可以加快后面进行网络结构学习的速度。
32.图1是本实施方式基于高阶动态贝叶斯的基因调控网络构建方法的流程示意图，如图1所示，所述基于高阶动态贝叶斯的基因调控网络构建方法包括如下步骤：
33.步骤1：根据时序基因表达数据，为每个目标基因在不同时间延迟下筛选潜在的调控因子形成每个目标基因的潜在调控因子集。
34.在本实施例中是从公开数据库dream challenge中下载的dream4时序基因表达数据。图2示出了基因数量为5，在2个时间延迟下寻找潜在调控因子的过程。图2中所包含的5个基因既是目标基因，同时也是彼此在进行网络结构学习搜索父节点时的待选基因。
35.本发明为目标基因寻找潜在调控因子的整个过程包括两种方法，目的是为了从不同的角度对基因之间的关联程度进行分析，以便可以为目标基因保留下来更多可能存在调控关系的其他基因。在方法一中，根据时序基因表达数据，使用互信息(mutual information,mi)去寻找目标基因和与目标基因具有不同时间延迟的其他基因之间的关联程度，但是互信息方法只能计算出基因之间关联程度的大小，不能确定到底是哪些基因应该被保留下来，哪些基因应该被删除，因此在这里使用预测最小描述长度(predictive minimum description length,pmdl)原理去进行阈值的确定，使得与目标基因关联程度大于阈值的基因被保留下来作为目标基因的潜在调控因子。在方法二中，利用时序基因表达数据使用皮尔逊相关系数(pearson correlation coefficient,pccs)去寻找目标基因和与目标基因具有不同时间延迟的其他基因之间的关联程度，同样用皮尔逊相关系数方法也是只能计算出基因之间关联程度的大小，不能确定到底是哪些基因应该被保留下来，哪些基因应该被删除，在这里使用基于假设检验的断点检测方法去进行阈值的确定，使得与目标基因关联程度大于阈值的基因被保留下来作为目标基因的潜在调控因子。最后，将上述两种方法分别为目标基因寻找到的潜在调控因子取并集，作为目标基因的潜在调控因子集。之所以使用两种方法去为目标基因进行潜在调控因子筛选并在最后将两种方法为目标基因寻找到的潜在调控因子取并集，是因为互信息和皮尔逊关联程度两种方法是从不同的角度对基因之间的关联程度进行衡量，该过程既可以揭示目标基因和与其具有不同时间延迟的其他基因间是否存在关系，也可以根据时间延迟来确定基因间的调控方向，这样做的目的就是为了减少对无关联或关联程度较小基因之间的计算，同时也为下面的构建高阶基因调控网络时为目标基因选择父节点缩小搜索空间，达到减小网络构建时间的目的。具体包括如下步骤：
36.步骤1.1：使用标准离散化方法对时序基因表达数据进行预处理，消除由于实验技术所导致的基因芯片的数据缺失或者数据异常，使各个样本和平行实验的数据处于相同的水平，并将预处理后的基因表达数据处理成矩阵形式，获得时序基因表达数据矩阵a，a的每一列表示一个基因，a的每一行表示基因在一个时间片的样本；
37.步骤1.2：根据为目标基因所要寻找的不同时间延迟下的潜在调控因子，在计算基因之间的相关关系之前，首先将目标基因和与目标基因具有不同时间延迟的其他基因在时序基因表达数据的选取上依据数据对齐原则进行数据对齐。
38.在本实施方式中，根据为目标基因所要寻找的不同时间延迟下的潜在调控因子，在计算基因之间的相关关系之前，首先将目标基因i和与目标基因i具有不同时间延迟下的其他基因j所对应矩阵a中的列向量ai与列向量aj依据数据对齐原则进行数据对齐。
39.所谓数据对齐，就是在计算目标基因i和与其有r个时间延迟的基因j之间的关联程度时，从时序基因表达数据中选取所使用的数据时，将目标基因i所对应的表达数据与要计算的与目标基因i具有r个时间延迟的基因j所对应的表达数据进行匹配。
40.具体来说，假设有一个长度为n的时间序列，要寻找d阶转移网络，则目标基因i所选取的数据是时间序列中序号d+1到n的数据，而与目标基因i具有r时间延迟的基因j选取的数据是时间序列中序号d-r+1到n-r的数据。如图3所示的对齐过程，假定x2是目标基因，有一个长度为n的时间序列，要构建2阶转移网络，因此要在2个时间延迟下寻找潜在调控因子，x1与x2具有1个时间延迟，在计算x1和x2之间的关联程度时，x1所选取的数据是时间序列中序号2到n-1的数据，x2选取的数据是时间序列中序号3到n的数据，同理x3表示的是与x2有2个时间延迟的基因节点，所选取的数据是时间序列中1到n-2的数据。
41.步骤1.3：根据数据对齐后的时序基因表达数据，结合互信息方法和预测最小描述长度方法为目标基因进行第一次筛选潜在调控因子；
42.具体是，根据对齐好的数据计算目标基因i与不同时间延迟下基因j之间的互信息值，形成一个互信息矩阵，然后分别以该矩阵中的每一个互信息值分别作为阈值进行筛选，得到一个模型矩阵，针对这个模型矩阵利用预测最小描述长度方法计算数据长度，选取具有最小长度的模型矩阵所对应的互信息值作为阈值对互信息矩阵进行筛选，将矩阵中大于阈值的互信息值保留，小于阈值的互信息值设为0，最后选择与目标基因之间互信息值大于0的基因作为目标基因的潜在调控因子。
43.步骤1.4：同样根据数据对齐后的时序基因表达数据，结合皮尔逊相关系数方法和基于假设检验的断点检测方法为目标基因进行第二次筛选潜在调控因子；
44.具体是，根据对齐好的数据计算基因i与不同时间延迟下基因j之间的皮尔逊相关系数值，形成一个皮尔逊相关系数矩阵，然后，假设计算的是d阶转移网络，此时分别对比目标基因i与1至(d-1)时间延迟下基因j的皮尔逊相关系数值，选最大值作为目标基因i与基因j的皮尔逊相关系数值，并对每个基因遍历该过程，然后将目标基因i与其他所有基因的皮尔逊系数值进行排序，基于假设检验原理进行断点检测，所寻找到的断点就是阈值，最后利用阈值对皮尔逊相关系数矩阵进行筛选，将矩阵中大于阈值的皮尔逊相关系数值保留，小于阈值的皮尔逊相关系数值设为0，最后选择与目标基因之间皮尔逊相关系数值大于0的基因作为目标基因的潜在调控因子。
45.步骤1.5：为了更多地为目标基因保留下来可能与其具有调控关系的基因节点，将步骤1.3和1.4计算所得到的目标基因的潜在调控因子集取并集，得到目标基因最终的潜在调控因子集。
46.步骤2：将每个目标基因的潜在调控因子集作为在利用高阶动态贝叶斯模型为每个目标基因寻找最优父节点集的搜索空间，构建局部最优的转移网络。
47.步骤2.1：依次遍历所述步骤1中获得的各个目标基因的潜在调控因子集，将各个目标基因的潜在调控因子集作为接下来进行结构学习的搜索空间。
48.步骤2.2：根据步骤2.1中各个目标基因搜索空间，利用高阶动态贝叶斯模型为每
个目标基因在搜索空间中选择不同的父节点组合并进行评分，且选出得分最高的父节点组合作为目标基因的父节点集，得到局部最优的转移网络。
49.步骤3：将各个目标基因对应的局部最优的转移网络进行合并形成全局最优的转移网络，并将其输出作为最终的基因调控网络。
50.虽然高阶动态贝叶斯基因调控网络的构建包括先验网络的构建和转移网络的构建，但是先验网络寻找时间片内的调控关系，先验网络可以设置为空；转移网络寻找时间片内和时间片间的调控关系，在生物背景下，时间片间的边对应于时延相互作用，时间片内的边对应于瞬时相互作用，只有当总时间延迟与采样间隙相比很小时，才被认为是瞬时交互，所以本发明只考虑时间片间的调控关系。
51.图4示出了基于5个基因构建2阶动态贝叶斯基因调控网络的转移网络的过程。首先基于5个基因分别得到由具有不同时间延迟的潜在调控因子所形成的潜在调控因子集后，分别对每个基因以该基因的潜在调控因子集作为搜索空间，进行高阶动态贝叶斯网络结构学习，最终将结构学习所得到的局部最优的转移网络进行合并形成全局最优的转移网络。
52.最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解；其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；因而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

技术特征：
1.一种基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，该方法包括如下步骤：步骤1：根据时序基因表达数据，为每个目标基因在不同时间延迟下筛选潜在的调控因子形成每个目标基因的潜在调控因子集；步骤2：将每个目标基因的潜在调控因子集作为在利用高阶动态贝叶斯模型为每个目标基因寻找最优父节点集的搜索空间，构建局部最优的转移网络；步骤3：将各个目标基因对应的局部最优的转移网络进行合并形成全局最优的转移网络，并将其输出作为最终的基因调控网络。2.根据权利要求1所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述步骤1包括如下步骤：步骤1.1：将时序基因表达数据处理成矩阵形式，其中所述矩阵的每一列表示一个基因；所述矩阵的每一行表示基因在一个时间片的样本；步骤1.2：根据为目标基因所要寻找的不同时间延迟下的潜在调控因子，将目标基因和与目标基因具有不同时间延迟的其他基因在时序基因表达数据的选取上依据数据对齐原则进行数据对齐；步骤1.3：根据数据对齐后的时序基因表达数据，结合互信息方法和预测最小描述长度方法为目标基因进行第一次筛选潜在调控因子；步骤1.4：同样根据数据对齐后的时序基因表达数据，结合皮尔逊相关系数方法和基于假设检验的断点检测方法为目标基因进行第二次筛选潜在调控因子；步骤1.5：将步骤1.3和1.4分别获得的目标基因的潜在调控因子集取并集，得到目标基因最终的潜在调控因子集。3.根据权利要求2所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述数据对齐，就是在计算目标基因i和与其有r个时间延迟的基因j之间的关联程度时，从时序基因表达数据中选取所使用的数据时，将目标基因i所对应的表达数据与要计算的与目标基因i具有r个时间延迟的基因j所对应的表达数据进行匹配。4.根据权利要求2所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述为目标基因进行第一次筛选潜在调控因子的方法为：首先根据对齐好的数据计算目标基因i和与目标基因i具有不同时间延迟的其他基因j之间的互信息值，形成一个互信息矩阵；然后分别以该矩阵中的每一个互信息值分别作为阈值进行筛选，得到一个模型矩阵；再然后针对这个模型矩阵利用预测最小描述长度方法计算数据长度，选取具有最小长度的模型矩阵所对应的互信息值作为阈值对互信息矩阵进行筛选，将矩阵中大于阈值的互信息值保留，小于阈值的互信息值设为0；最后选择与目标基因之间互信息值大于0的基因作为目标基因的潜在调控因子。5.根据权利要求2所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述为目标基因进行第二次筛选潜在调控因子的方法为：首先根据对齐好的数据计算基因i和与目标基因i具有不同时间延迟的其他基因j之间的皮尔逊相关系数值，形成一个皮尔逊相关系数矩阵；然后，假设计算的是d阶转移网络，此时分别对比目标基因i与1至(d-1)时间延迟下的基因j的皮尔逊相关系数值，选最大值作为目标基因i与基因j的皮尔逊相关系数值，并对每个基因遍历该过程；再然后将目标基因i与其他所有基因的皮尔逊系数值进行
排序，基于假设检验原理进行断点检测，所寻找到的断点就是阈值；又然后利用阈值对皮尔逊相关系数矩阵进行筛选，将矩阵中大于阈值的皮尔逊相关系数值保留，小于阈值的皮尔逊相关系数值设为0；最后选择与目标基因之间皮尔逊相关系数值大于0的基因作为目标基因的潜在调控因子。6.根据权利要求1所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述步骤2包括如下步骤：步骤2.1：依次遍历各个目标基因的潜在调控因子集，将各个目标基因的潜在调控因子集作为接下来进行结构学习的搜索空间；步骤2.2：根据各个目标基因的搜索空间，利用高阶动态贝叶斯模型为每个目标基因在搜索空间中选择不同的父节点组合并进行评分，且从中选出得分最高的父节点组合作为目标基因的父节点集，得到局部最优的转移网络。7.根据上述任一项权利要求所述的基于高阶动态贝叶斯的基因调控网络构建方法，其特征在于，所述时序基因表达数据是从公开数据库dream challenge中下载的dream4时序基因表达数据。

技术总结
本发明公开了一种基于高阶动态贝叶斯的基因调控网络构建方法，涉及医学信息学技术领域。通过利用时序基因表达数据去构建基于时间变量的高阶动态贝叶斯基因调控网络，在不同的时间延迟下为目标基因寻找父节点，以提高网络结构构建的准确率。由于高阶网络结构学习的搜索空间非常大，导致结构学习的时间复杂度特别高，为了加快学习速度，提出了在结构学习之前进行潜在调控因子筛选的办法，首先筛选出与目标基因关联度较高的基因作为目标基因的潜在调控因子集，删除与目标基因关联程度较低的基因，以达到缩小搜索空间的目的，可以加快后面进行网络结构学习的速度。进行网络结构学习的速度。进行网络结构学习的速度。


技术研发人员：李婵 王之琼 曲璐渲 王炜祎祺 信俊昌
受保护的技术使用者：东北大学
技术研发日：2021.12.13
技术公布日：2022/3/8