一种面向大规模基因调控网络的智能推断方法与流程

1.本发明涉及基因调控领域，尤其涉及一种面向大规模基因调控网络的智能推断方法。


背景技术：

2.基因调控网络是由基因之间相互作用形成的动态生化网络，能够从系统的角度认识、理解生物遗传信号的传输、细胞的分裂等活动规律，对于阐明细胞周期、损伤修复和凋亡的生物学机制至关重要。此外，基因调控网可以分析和识别病变细胞的关键调控基因，有助于肿瘤等复杂疾病的诊断、治疗以及靶向药物的研发。因此，推断基因调控网络是系统生物学的一项重要研究内容。随着高通量测序技术的发展，研究人员已经获得了大量的生物基因表达数据，例如大肠杆菌数据集包含的基因数量接近4000个，为基因调控网络的推断奠定了基础。近年来使用人工智能方法和基因表达数据推断基因调控网络成为生物信息学研究的热点。
3.目前已有的面向大规模基因表达数据的基因调控网络推断算法，主要存在两个缺点，一是推断算法的计算效率较低，二是推断算法的准确性有待提高。主要原因是基因调控网络具有高维度、非线性等特点，使得大规模基因调控网络的推断非常困难。


技术实现要素：

4.本发明提供一种名称，以克服计算效率较低、准确性不高等技术问题。
5.为了实现上述目的，本发明的技术方案是：
6.一种面向大规模基因调控网络的智能推断方法，其特征在于，包括：
7.步骤1，在基因表达数据点分布的两维空间中，使用m*n的网格划分数据空间，使落在第(x，y)格子中的数据点的频率作为联合概率分布函数p(x，y)的估计值o(x，y)；
8.步骤2，利用联合概率分布函数p(x，y)的估计值o(x，y)确定基因x和基因y表达值之间的互信息i(x；y)；
9.步骤3，在不同网格中利用互信息i(x；y)计算最大互信息系数mic(x；y)，筛选出调控基因集x中关于靶基因y的关键调控基因集r；
10.步骤4，利用关键调控基因重要性融合算法计算关键调控基因集r中的调控基因对靶基因的重要性，获得重要性矩阵i；
11.步骤5，在重要性矩阵i中保留控制强度超过设定阈值的调控关系，控制强度超过设定阈值的调控关系构成基因调控网络。
12.进一步的，步骤4具体为：
13.步骤4.1，整合基因真实表达数据中的时间序列数据和平稳状态数据及关键调控基因集r，获得用于建模的基因表达数据集；
14.步骤4.2，将基因表达数据集分解为p个子集，每一个子集均有一个靶基因；
15.步骤4.3，在每个子集中，分别运用xgboost、lightgbm和catboost算法建立树模型
学习非线性函数，计算得到关键调控基因对靶基因的重要性值；
16.步骤4.4，将关键调控基因对靶基因的重要性值排序，得到一个有序的关键调控基因的重要性列表f1；
17.步骤4.5，重复4.3和4.4，依据基因编号将p个重要性列表合成为重要性矩阵i。
18.进一步的，步骤3筛选出调控基因集x中关于靶基因y的关键调控基因集r具体为：
19.步骤3.1：给定基因真实表达数据，其中基因真实表达数据共有s个样本和g个基因，依次选取一个基因作为靶基因j，剩下的基因作为候选调控基因；
20.步骤3.2：获取给定靶基因j与候选调控基因之间的最大信息系数mic(x；y)，选取阈值剔除冗余调控基因，得到给定靶基因的关键调控基因列表r
j
；
21.步骤3.3：重复步骤3.1和步骤3.2，得到所有靶基因的关键调控基因集r。
22.进一步的，步骤4.3中在每个子集中，分别运用xgboost、lightgbm和catboost算法建立树模型学习非线性函数，计算得到关键调控基因对靶基因的重要性具体为：
23.步骤4.3.1，利用xgboost算法获得第一重要性列表f
xgb
；
24.步骤4.3.2，利用lightgbm算法获得第二重要性列表f
lg
；
25.步骤4.3.3，利用catboost算法获得第三重要性列表f
cb
；
26.步骤4.3.4，由第一重要性列表f
xgb
、第二重要性列表f
lg
和第三重要性列表f
cb
构成关键调控基因对靶基因的重要性值。
27.进一步的，步骤4.4将关键调控基因对靶基因的重要性值排序，得到一个有序的关键调控基因的重要性列表f1具体公式为：
28.f1＝α1·
f
xgb
α2·
f
lg
α3·
f
cb
29.其中，α1、α2、α3为权重系数。
30.进一步的，步骤1中获取联合概率分布函数p(x，y)的估计值o(x，y)的公式为：
[0031][0032]
其中，g为第(x，y)个网格中的数据点数，g为总的数据点数。
[0033]
进一步的，步骤2中确定基因x和基因y表达值之间的互信息i(x；y)的公式为：
[0034][0035]
其中，x和y分别为基因x和基因y表达值的列向量，p(x)和p(y)是x和y的边缘概率分布函数。
[0036]
进一步的，步骤3在不同网格中利用互信息i(x；y)获得基因间最大互信息系数mic(x；y)的公式为：
[0037][0038]
其中min(x，y)代表取x与y的最小值，其中x为网格x方向被分段的数量，y为网格y方向被分段的数量。
[0039]
有益效果：本方法针对大规模基因调控网络中关键调控基因的稀疏性，首先构建
关键调控基因识别算法降低候选基因维度，排除冗余基因表达信息，有效提高了计算效率；之后再利用关键调控基因重要性融合算法和基因调控关系约简算法来准确推断基因间的调控关系，可以提高推断调控关系的准确度。
附图说明
[0040]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0041]
图1为本发明推断方法流程图；
[0042]
图2为使用本发明的基因间最大互信息系数仿真效果图。
具体实施方式
[0043]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0044]
本实施例提供了一种面向大规模基因调控网络的智能推断方法，如图1，其特征在于，包括：
[0045]
步骤1，在基因表达数据点分布的两维空间中，使用m*n的网格划分数据空间，使落在第(x，y)格子中的数据点的频率作为联合概率分布函数p(x，y)的估计值o(x，y)；
[0046]
具体的，获取联合概率分布函数p(x，y)的估计值o(x，y)的公式为：
[0047][0048]
其中，g为第(x，y)个网格中的数据点数，g为总的数据点数。
[0049]
步骤2，利用联合概率分布函数p(x，y)的估计值o(x，y)确定基因x和基因y表达值之间的互信息i(x；y)；
[0050]
具体的，确定基因x和基因y表达值之间的互信息i(x；y)的公式为：
[0051][0052]
其中，x和y分别为基因x和基因y表达值的列向量，p(x)和p(y)是x和y的边缘概率分布函数。
[0053]
因为m*n的网格划分数据点的方式不止一种，所以取使互信息最大的网格划分。然后使用归一化因子，将互信息的值映射到(0，1)区间。最后，找到使归一化互信息最大的网格分辨率，作为最大互信息系数的度量值。其中网格的分辨率限制为m*n＜b，b＝data_size0.6，其中b代表网格分辨率最大值。
[0054]
步骤3，在不同网格中利用互信息i(x；y)计算最大互信息系数mic(x；y)，筛选出调控基因集x中关于靶基因y的关键调控基因集r；
[0055]
具体的，步骤3可分为：
[0056]
步骤3.1：给定基因真实表达数据，其中基因真实表达数据共有s个样本和g个基因，依次选取一个基因作为靶基因j，剩下的基因作为候选调控基因；
[0057]
步骤3.2：获取给定靶基因j与候选调控基因之间的最大信息系数mic(x；y)，选取阈值剔除冗余调控基因，得到给定靶基因的关键调控基因列表r
j
；
[0058]
步骤3.3：重复步骤3.1和步骤3.2，得到所有靶基因的关键调控基因集r，其中r＝{r1，...，r
j
...，r
g
}，从而实现对大规模基因表达数据降维的作用。
[0059]
其中，计算最大互信息系数mic(x；y)的公式为：
[0060][0061]
其中min(x，y)代表取x与y的最小值，其中x为网格x方向被分段的数量，y为网格y方向被分段的数量。
[0062]
面向大规模基因表达数据时，若不能在推断基因调控关系之前识别出关键的调控基因，意味着数据集中包含的所有基因均为潜在的调控基因，这必然会导致推断算法的计算效率和准确性不高。但利用本方法步骤1
‑
步骤3的关键调控基因识别算法就能够有效提高基因调控网络推断算法的计算效率和准确性。
[0063]
利用以上关键调控基因识别算法可以确定靶基因j的关键调控基因集合r，在此基础上，采用关键调控基因重要性融合算法进而确定集合r中的调控基因对基因j的调控强度。关键调控基因重要性融合算法是一种基于信息增益算法来计算候选调控基因对靶基因存在调控关系的置信度，定义为调控作用的强度，即调控强度。
[0064]
关键调控基因重要性融合算法，是推断基因调控关系的核心算法。关键调控基因重要性融合算法采用三种高效的梯度提升树模型来计算调控基因对基因j的调控强度，具体见步骤4。
[0065]
步骤4，利用关键调控基因重要性融合算法计算关键调控基因集r中的调控基因对靶基因的重要性，获得重要性矩阵i；
[0066]
其中，步骤4具体为：
[0067]
步骤4.1，整合基因真实表达数据中的时间序列数据和平稳状态数据及关键调控基因集r，获得用于建模的基因表达数据集；
[0068]
步骤4.2，将基因表达数据集分解为p个子集，每一个子集均有一个靶基因；
[0069]
步骤4.3，在每个子集中，分别运用xgboost、lightgbm和catboost算法建立树模型学习非线性函数，计算得到关键调控基因对靶基因的重要性值；
[0070]
其中，步骤4.3中在每个子集中，分别运用xgboost、lightgbm和catboost算法建立树模型学习非线性函数，计算得到关键调控基因对靶基因的重要性值具体为：
[0071]
步骤4.3.1，利用xgboost算法获得第一重要性值f
xgb
；
[0072]
步骤4.3.2，利用lightgbm算法获得第二重要性值f
lg
；
[0073]
步骤4.3.3，利用catboost算法获得第三重要性值f
cb
；
[0074]
步骤4.3.4，由第一重要性列表f
xgb
、第二重要性列表f
lg
和第三重要性列表f
cb
构成关键调控基因对靶基因的重要性值。
[0075]
这三种算法都可以使用多核cpu进行加速学习，从而方便地实现关键调控基因重要性融合算法的多线程并行计算。
[0076]
步骤4.4，将关键调控基因对靶基因的重要性值排序，得到一个有序的关键调控基因的重要性列表f1。
[0077]
步骤4.4中将关键调控基因对靶基因的重要性值累加后排序，得到一个有序的关键调控基因的第一重要性列表f1具体公式为：
[0078]
f1＝α1·
f
xgb
α2·
f
lg
α3·
f
cb
[0079]
其中，α1、α2、α3为权重系数。
[0080]
在构建回归树过程中，计算出调控基因对靶基因的重要性值，然后将三种重要性乘以相应的权重系数α并累加排序，得到调控基因对靶基因的重要性列表f。权重系数α可根据经验尝试更改，默认值可设为1。
[0081]
将重要性值作为调控基因对靶基因的调控强度，最终得到所有基因间的调控关系。在构建回归树过程中，需使用现有的贪婪算法进行增益计算和树的分裂。
[0082]
步骤4.5，重复4.3和4.4，依据基因编号将p个重要性列表合成为重要性矩阵i；
[0083]
其中，重要性矩阵i＝{i1，...，i
j
，...，i
p
}。在重要性矩阵中，元素的位置坐标反映了基因调控的方向，元素的数值大小反映了基因的调控强度。
[0084]
步骤5，在重要性矩阵i中保留控制强度超过设定阈值的调控关系，控制强度超过设定阈值的调控关系构成基因调控网络。
[0085]
步骤5具体为：首先从重要性矩阵i中识别出所有控制强度超过设定阈值的调控关系。阈值计算规则为淘汰调控强度较小的调控关系，保留调控强度较大的调控关系。定义i
i，j
为调控基因i对靶基因j的重要性值，t为阈值。如果i
i，j
≤t，则令i
i，j
＝0，将所有i
i，j
为0的调控关系剔除，即阈值规则会淘汰重要性较小的直接作用的调控关系，从而实现基因约简。若i
i，j
≠0，意味着基因i对基因j存在调控作用，且调控的强度为i
i，j
。如此便得到了基因调控网络。
[0086]
本发明计算两个基因间的最大互信息系数时，需要选取最优网格得到最大的互信息系数，这一过程是完全是算法自动完成的。
[0087]
图2为本发明仿真图，图中底部分别为x轴和y轴，数值都是从0到100，代表100个基因，整张图是关于y＝x对称的。z轴表示相关系数大小，数值越大代表越相关，可以设置阈值删除靶基因的一些低相关候选调控基因，得到关键调控基因。
[0088]
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。