一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法

1.本发明属于生物信息学领域，具体涉及一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法。

背景技术：

2.在药物从设计到上市的整个研究过程中，毒性通常是导致失败的主要原因。如今，化合物的毒性评价已经与其药效学研究具有同等重要的地位。然而，传统的药物毒性测试方法通常需要使用动物模型。这些动物试验不仅耗费大量时间，而且成本很高，还涉及到动物福利等问题。近年来，虽然可以通过体外高通量筛选等方法以较快的速度评价化合物的毒性，但是基于实验的评价方法依然代价大、周期长。而一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，利用结构信息对化合物的毒性进行预测，就可以保证在进行体外和体内试验之前对化合物进行大规模毒性筛选，在药物开发早期就筛除有风险的化合物，将能更好地解决候选药物具有毒性的问题。

技术实现要素：

3.本发明的目的是提供一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，能够有效地提高对化合物致突变性、心脏毒性等毒性的预测正确率。
4.为实现上述目的，本发明采用如下的技术方案：一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，包括以下步骤：
5.步骤1：对数据进行预处理，生成分子图；
6.步骤2：提取分子的异构边特征；
7.步骤3：将异构边特征嵌入分子图，并给予权重；
8.步骤4：构建图注意力网络并使用新的分子图建模，使用训练好的模型进行预测。
9.所述步骤1具体为：
10.步骤1.1：对用于毒性预测模型训练集、测试集、验证集数据中化合物的smiles字符串进行预处理，包括去盐，去重，去除混合物、聚合物，去除含有重金属的化合物，去除不能生成分子图的化合物；
11.步骤1.2：根据smiles字符串初始化化合物的分子图。
12.所述步骤2具体为：
13.步骤2.1：以分子中的某一个化学键为起点，出分子中所有的三元组，三元组包括邻居原子、键、目标原子；
14.步骤2.2：根据三元组中的元素组成确定分子异构边，即如果存在两个或多个三元组中，邻居原子和目标原子相同，键不同，则这样的键构成分子中的异构边。
15.所述步骤3具体为：
16.步骤3.1：生成异构边特征矩阵，将异构边特征嵌入分子图；
17.步骤3.2：对目标原子和邻居原子的状态向量使用leaky_relu激活函数进行非线性激活得到异构边特征向量和原子对特征向量，然后利用softmax函数给予异构边特征向量一定权重，计算方法如下：
18.h
vu
＝leaky_relu(w
·
[gv,gu])
ꢀꢀꢀ
(15)
[0019]evu
＝leaky_relu(w
·
[hv,hu])
ꢀꢀꢀ
(16)
[0020][0021]
其中，gv,hv是目标原子的状态向量，gu,hu是邻居原子的状态向量；w是可训练的权值矩阵；h
vu
是激活函数激活之后的异构边特征向量，e
vu
是经过leaky_relu激活函数激活之后的普通边特征向量；n(v)是目标原子的所有邻居原子的集合；a
vu
是经过softmax函数加权之后得到的邻居原子u对目标原子v的权值。
[0022]
所述步骤4具体为：
[0023]
步骤4.1：将邻居原子u对目标原子v的权值a
vu
，使用elu激活函数进行非线性激活(断句完成)，得到目标原子v的上下文向量cv，计算方法如下：
[0024][0025]
其中w是可训练的权值矩阵；hu是邻居原子的状态向量；n(v)是目标原子的所有邻居原子的集合。
[0026]
步骤4.2：邻居原子u利用图注意力机制不断地为目标原子v聚合信息，计算方法如下：
[0027][0028]
其中m
k-1
是迭代k-1次的消息传递函数；是邻居原子u迭代k-1次的状态向量；是目标原子v迭代k-1次的状态向量；n(v)是目标原子的所有邻居原子的集合；是迭代k-1次的目标原子的聚合上下文特征向量。
[0029]
步骤4.3：更新目标原子的状态向量，获取目标原子v的前一个状态向量以及聚合上下文向量作为输入，将前一个状态更新为当前状态计算方法如下：
[0030][0031]
其中gru
k-1
是迭代k-1次的更新函数；是迭代k-1次的目标原子的聚合上下文特征向量；是目标原子v迭代k-1次的状态向量；是目标原子v迭代k次的状态向量。
[0032]
步骤4.4：将迭代k次的目标原子v的状态向量作为输入，置于网络的输出层，预测分子的毒性，计算方法如下：
[0033][0034]
其中w是根据原子特征向量的维度，设置的降维矩阵；是目标原子v迭代k次的状态向量；output是输出层；molfeature是最终预测的结果。
[0035]
所述毒性预测模型如下：
[0036]hvu
＝leaky_relu(w
·
[gv,gu])
ꢀꢀꢀ
(22)
[0037]evu
＝leaky_relu(w
·
[hv,hu])
ꢀꢀꢀ
(23)
[0038][0039][0040][0041][0042][0043]
其中，计算异构边特征向量为公式(8)，计算原子对特征向量为公式(9)，计算邻居原子对目标原子的权值为公式(10)，计算目标原子的上下文特征向量为公式(11)，计算目标原子的聚合上下文特征向量为公式(12)，计算目标原子的最新状态向量为公式(13)，根据目标原子的最新状态向量对化合物进行毒性预测为公式(14)。
[0044]
本发明与现有技术相比具有以下有益效果：
[0045]
1、本发明的建模流程简洁明了，无需计算冗余的分子指纹数据，而是使用结构更加清晰的分子图建模，对于相同的化合物，预测正确率更高。
[0046]
2、本发明使用图注意力神经网络对分子图建模，相比于传统的机器学习算法建模，模型预测的准确率显著提高。
[0047]
2、本发明在建模过程中考虑到了分子中的异构边特征，并且将异构边特征融合进分子图中，提高了异构边在整个分子图中的权重，更好地提高了对化合物致突变性、心脏毒性等毒性的预测正确率。
[0048]
3、本发明在建模过程中使用的训练集、测试集、验证集数据来源广泛，模型的泛化性较好。
附图说明
[0049]
图1建模流程图；
[0050]
图2化合物致突变性模型在外部验证集上的roc图；
[0051]
图3化合物herg心脏毒性模型在外部验证集1上的rmse图；
[0052]
图4化合物herg心脏毒性模型在外部验证集2上的rmse图；
[0053]
图5化合物herg心脏毒性模型在外部验证集3上的rmse图；
[0054]
图6是本发明的化合物herg心脏毒性模型在外部验证集4上的rmse图。
具体实施方式
[0055]
下面结合附图及实例对本发明做进一步说明。
[0056]
实施例1：
[0057]
下面结合附图及实例对本发明做进一步说明。
[0058]
如图1所示，本发明的具体实施例流程为：
[0059]
步骤1：对于心脏毒性分子，首先从chembl数据库中获得11512个心脏毒性分子化合物作为训练集，从pubchem bioassay数据集、herg central数据库、open source malaria数据库和文献中选取共1161个心脏毒性分子化合物作为外部验证集。使用rdkit化学工具包去盐，将化合物上传到ochem网站去重，去除不能生成分子图结构的化合物，最后得到训练数据9200个，外部验证数据1146个。对于致突变性分子，从ames致突变性基准数据集中获得致突变性分子化合物作为训练集，从ccris数据库、国家毒理学计划ntp数据库和桑尼塔伤寒沙门菌高等研究所isssty数据库中获得致突变性分子化合物作为外部验证集，过滤了混合物、聚合物、含有重金属的化合物以及碳原子少于3个的化合物，去除不能生成分子图结构的化合物，最后得到训练数据6506个，外部验证数据1471个。
[0060]
步骤2：对于每一个smiles字符串，生成相应的分子图，使用rdkit化学工具包出分子中所有的键以及键两端的原子，构成邻居原子，键，目标原子三元组，根据所有三元组确定分子中的异构边，即邻居原子和目标原子相同，键不同的三元组，其中不同的键构成分子中的异构边。
[0061]
步骤3：将分子键初始化为全0矩阵，异构边初始化为1，将异构边特征嵌入分子图，对目标原子和邻居原子的状态向量使用leaky_relu激活函数进行非线性激活得到异构边特征向量和原子对特征向量，然后利用softmax函数给予异构边特征向量一定权重，计算方法如下：
[0062]hvu
＝leaky_relu(w
·
[gv,gu])
ꢀꢀꢀ
(1)
[0063]evu
＝leaky_relu(w
·
[hv,hu])
ꢀꢀꢀ
(2)
[0064][0065]
其中gv,hv和gu,hu分别是目标原子和邻居原子的状态向量；w是可训练的权值矩阵；h
vu
和e
vu
分别是经过leaky_relu激活函数激活之后的异构边特征向量和原子对特征向量；n(v)是目标原子的所有邻居原子的集合；a
vu
是经过softmax函数加权之后得到的邻居原子u对目标原子v的重要性，即为权值。
[0066]
步骤4：将邻居原子u对目标原子v的重要性，即权值a
vu
，使用elu激活函数进行非线性激活，得到目标原子v的上下文向量cv，计算方法如下：
[0067][0068]
其中w是可训练的权值矩阵；hu是邻居原子的状态向量；n(v)是目标原子的所有邻居原子的集合。
[0069]
邻居原子u利用图注意力机制不断地为目标原子v聚合信息，计算方法如下：
[0070][0071]
其中m
k-1
是迭代k-1次的消息传递函数；是邻居原子u迭代k-1次的状态向
量；是目标原子v迭代k-1次的状态向量；n(v)是目标原子的所有邻居原子的集合；是迭代k-1次的目标原子的聚合上下文特征向量。
[0072]
更新目标原子的状态向量，获取目标原子v的前一个状态向量以及聚合上下文向量作为输入，将前一个状态更新为当前状态计算方法如下：
[0073][0074]
其中gru
k-1
是迭代k-1次的更新函数；是迭代k-1次的目标原子的聚合上下文特征向量；是目标原子v迭代k-1次的状态向量；是目标原子v迭代k次的状态向量。
[0075]
将迭代k次的目标原子v的状态向量作为输入，置于网络的输出层，预测分子的毒性，计算方法如下：
[0076][0077]
其中w是根据原子特征向量的维度，设置的降维矩阵；是目标原子v迭代k次的状态向量；output是输出层；molfeature是最终预测的结果。
[0078]
进一步的，所述毒性预测模型具体为：
[0079]hvu
＝leaky_relu(w
·
[gv,gu])
ꢀꢀꢀ
(8)
[0080]evu
＝leaky_relu(w
·
[hv,hu])
ꢀꢀꢀ
(9)
[0081][0082][0083][0084][0085][0086]
其中，计算异构边特征向量为公式8，计算原子对特征向量为公式9，计算邻居原子对目标原子的权值为公式10，计算目标原子的上下文特征向量为公式11，计算目标原子的聚合上下文特征向量为公式12，计算目标原子的最新状态向量为公式13，根据目标原子的最新状态向量对化合物进行毒性预测为公式14。
[0087]
实例1：
[0088]
给定一个化合物，取自ccris数据库，cas号：820-75-7，要预测其致突变性，首先根据其smiles字符串生成分子图，使用rdkit化学工具包出化合物中所有的异构边，然后初始化异构边特征并嵌入分子图，根据其分子图数据，使用本发明的致突变性模型进行预测。根据当前状态判断结果为1，说明此化合物有致突变性。预测结果与实验结果相同。
[0089]
实例2：
[0090]
给定一个化合物，取自chembl数据库，id号：chembl100415，要预测其herg心脏毒
性，首先根据其smiles字符串生成分子图，使用rdkit化学工具包出化合物中所有的异构边，然后初始化异构边特征并嵌入分子图，根据其分子图数据，使用本发明的herg心脏毒性模型进行预测。根据当前状态判断pic50值为6.89，而其真实pic50值为6.72。预测结果与实验结果相近。
[0091]
本发明的有效性验证：对于化合物致突变性的预测，如表1所示，本发明在外部验证集上的auc为0.842。本发明与其他对于化合物致突变性的预测方法进行了比较，其中在mutagenpred-gcnns方法中，外部验证集的最好auc为0.838。在专利号为cn112466399a的专利中，外部验证集上的auc为0.797。本发明在化合物致突变性的预测上有所提高。对于化合物herg引起的心脏毒性的预测，如表2所示，本发明在外部验证集1到4上的rmse分别为0.948、0.908、0.810、1.060，验证集平均值为0.9315，本发明与其他对于化合物herg引起的心脏毒性的预测方法进行了比较，其中在tssf-herg方法中外部验证集1到4上的rmse分别为0.876、0.661、1.117、1.117，验证集平均值为0.94275，本发明虽然在验证集1和2上的性能有所降低，但是在验证集3和4上的性能有所提高。本发明在化合物herg心脏毒性的预测上有所提高。
[0092]
表1在致突变性数据集上与其他模型的性能比较
[0093]
化合物致突变性的预测方法外部验证集auc本发明0.842mutagenpred-gcnns0.838专利号：cn112466399a0.797
[0094]
表2在心脏毒性数据集上与其他模型的性能比较
[0095]
化合物herg心脏毒性的预测方法验证集1验证集2验证集3验证集4验证集平均值本发明0.9480.9080.8101.0600.9315tssf-herg0.8760.6611.1171.1170.94275
。

技术特征：

1.一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，包括以下步骤：步骤1：对数据进行预处理，生成分子图；步骤2：提取分子的异构边特征；步骤3：将异构边特征嵌入分子图，并给予权重；步骤4：构建图注意力网络并使用新的分子图建模，使用训练好的模型进行预测。2.根据权利要求1所述的一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，所述步骤1具体为：步骤1.1：对用于毒性预测模型训练集、测试集、验证集数据中化合物的smiles字符串进行预处理，包括去盐，去重，去除混合物、聚合物，去除含有重金属的化合物，去除不能生成分子图的化合物；步骤1.2：根据smiles字符串初始化化合物的分子图。3.根据权利要求1所述的一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，所述步骤2具体为：步骤2.1：以分子中的某一个化学键为起点，出分子中所有的三元组，三元组包括邻居原子、键、目标原子；步骤2.2：根据三元组中的元素组成确定分子异构边，即如果存在两个或多个三元组中，邻居原子和目标原子相同，键不同，则这样的键构成分子中的异构边。4.根据权利要求1所述的一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，所述步骤3具体为：步骤3.1：生成异构边特征矩阵，将异构边特征嵌入分子图；步骤3.2：对目标原子和邻居原子的状态向量使用leaky_relu激活函数进行非线性激活得到异构边特征向量和原子对特征向量，然后利用softmax函数给予异构边特征向量一定权重，计算方法如下：h
vu
＝leaky_relu(w
·
[g
v
,g
u
])
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)e
vu
＝leaky_relu(w
·
[h
v
,h
u
])
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)其中，g
v
,h
v
是目标原子的状态向量，g
u
,h
u
是邻居原子的状态向量；w是可训练的权值矩阵；h
vu
是激活函数激活之后的异构边特征向量，e
vu
是经过leaky_relu激活函数激活之后的普通边特征向量；n(v)是目标原子的所有邻居原子的集合；a
vu
是经过softmax函数加权之后得到的邻居原子u对目标原子v的权值。5.根据权利要求1所述的一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，所述步骤4具体为：步骤4.1：将邻居原子u对目标原子v的权值a
vu
，使用elu激活函数进行非线性激活(断句完成)，得到目标原子v的上下文向量c
v
，计算方法如下：
其中w是可训练的权值矩阵；h
u
是邻居原子的状态向量；n(v)是目标原子的所有邻居原子的集合。步骤4.2：邻居原子u利用图注意力机制不断地为目标原子v聚合信息，计算方法如下：其中m
k-1
是迭代k-1次的消息传递函数；是邻居原子u迭代k-1次的状态向量；是目标原子v迭代k-1次的状态向量；n(v)是目标原子的所有邻居原子的集合；是迭代k-1次的目标原子的聚合上下文特征向量。步骤4.3：更新目标原子的状态向量，获取目标原子v的前一个状态向量以及聚合上下文向量作为输入，将前一个状态更新为当前状态计算方法如下：其中gru
k-1
是迭代k-1次的更新函数；是迭代k-1次的目标原子的聚合上下文特征向量；是目标原子v迭代k-1次的状态向量；是目标原子v迭代k次的状态向量。步骤4.4：将迭代k次的目标原子v的状态向量作为输入，置于网络的输出层，预测分子的毒性，计算方法如下：其中w是根据原子特征向量的维度，设置的降维矩阵；是目标原子v迭代k次的状态向量；output是输出层；molfeature是最终预测的结果。6.根据权利要求2所述的一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，所述毒性预测模型如下：h
vu
＝leaky_relu(w
·
[g
v
,g
u
])
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)e
vu
＝leaky_relu(w
·
[h
v
,h
u
])
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)(9)(9)(9)(9)其中，计算异构边特征向量为公式(8)，计算原子对特征向量为公式(9)，计算邻居原子对目标原子的权值为公式(10)，计算目标原子的上下文特征向量为公式(11)，计算目标原子的聚合上下文特征向量为公式(12)，计算目标原子的最新状态向量为公式(13)，根据目标原子的最新状态向量对化合物进行毒性预测为公式(14)。

技术总结

一种通过嵌入分子异构边特征的图注意力网络预测化合物毒性的方法，包括以下步骤：步骤1、对数据进行预处理，生成分子图；步骤2、提取分子的异构边特征；步骤3、将异构边特征嵌入分子图，并给予一定的权重；步骤4、构建图注意力网络并使用新的分子图建模，将训练好的模型用于预测。本发明考虑到了分子中的异构边特征，利用图注意力网络对分子图建模，将异构边特征融合进分子图中，更好地提高了对化合物致突变性、心脏毒性等毒性的预测正确率。心脏毒性等毒性的预测正确率。