预测小分子CYP代谢酶的代谢位点的方法及系统

摘要

一种预测小分子CYP代谢酶的代谢位点的方法及系统，包括：收集与代谢位点相关的分子结构数据，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集；将二维分子结构转化为三维分子结构；将分子三维结构中的原子类型的矢量化；将分子三维结构中距离的矢量化；将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率；在不同CYP代谢酶亚型上根据为代谢位点的概率预测代谢位点；上述方法及系统，对分子与不同亚型的CYP代谢发生化学反应的代谢位点预测，该预测结果可以帮助药物化学家设计或者优化分子结构。

说明书

技术领域

本发明涉及计算机技术领域，特别涉及一种预测小分子CYP代谢酶的代谢位点的方法及系统。

背景技术

潜在药物的ADMET(吸收，分布，代谢，排泄，毒性)是成药性属性预测重要的组成部分。肝脏是药物的主要清除器官，肝脏清除分成肝脏代谢和胆汁排泄两种方式。肝脏富含药物Ⅰ相代谢和Ⅱ相代谢所需的各种酶，其中以P450酶最为重要。P450酶是由多种类型的P450酶所组成的一个大家族，根据氨基酸的排序的雷同性，P450酶可以分为不同几个大类，每个大类又可以细分成几个小类。在人体中重要的P450酶有CYP1A2、CYP2A6、CYP2B6、CYP2C8、CYP2C9、CYP2C19、CYP2D6、CYP2E1、CYP3A4和CYP3A5。P450酶存在有明显的种属差异，药物在动物和人体内的代谢途径和代谢产物可能是不同的。多态性(polymorphisms)是P450酶的一个重要特征，是导致药物反应的个体差异的一个重要原因。所谓的多态性，是指同一种属的不同个体间某一P450酶的量存在较大的差异。量高的个体代谢速度就快，称为快代谢型(extensive metabolizer)；量低的个体代谢速度就慢，称为慢代谢型(poormetabolizer)。人体内许多P450酶表现出多态性，其中以CYP2D6和CYP2C19的多态性最为典型。另外，P450酶具有可诱导和可抑制性。也就是说，P450酶的量和活性会受到药物(或其他外源物)的影响，可能会影响药物本身的代谢，并可能会引起代谢性药物相互作用。在药物设计的过程中，潜在药物可能因为代谢速率过快，导致没办法在靶向器官达到有效药物作用浓度，而没办法发挥药效。潜在药物也有可能因为代谢速度过慢，而导致药物累积在代谢器官-肝脏，从而导致药物在肝脏毒性。所以在药物设计的早期阶段，根据药物在动物体内的代谢表现，优化分子结构以达到有效浓度而不致于积累肝脏毒性。

发明内容

基于此，有必要提供一种提高预测效果的预测小分子CYP代谢酶的代谢位点的方法。

同时，提供一种提高预测效果的预测小分子CYP代谢酶的代谢位点的系统。

一种预测小分子CYP代谢酶的代谢位点的方法，包括：

整理数据：收集与代谢位点相关的分子结构数据，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集；

结构转化：将二维分子结构转化为三维分子结构；

原子类型的矢量化：将分子三维结构中的原子类型矢量化，将分子中所有原子的矢量依据原子排序矢量化形成矩阵M；

距离的矢量化：将分子三维结构中的距离矢量化，计算分子中两两原子距离，形成矩阵D；

构建模型：在图卷积神经网络模型中，输入原子类型的矢量化与距离的矢量化，将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率；

预测代谢位点：在不同CYP代谢酶亚型上预测代谢位点原子，在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名靠前的原子预测为代谢位点。

在优选的实施例中，还包括：添加标签：根据预测结果对原子标记0,1标签值，0代表不是代谢位点，1代表是代谢位点。

在优选的实施例中，所述整理数据步骤中，基于不同的代谢酶亚型将收集来的数据整理成为9种不同亚型的数据集：CYP 1A2、CYP 2A6、CYP 2B6、CYP 2C19、CYP 2C8、CYP2C9、CYP 2D6、CYP 2E1、CYP 3A4，对分子结构进行去盐、去溶剂处理，整理成为标准的模型输入格式，若分子设置有标签则在分子的标签值上，将原子转化为0,1标签值，0代表不是代谢位点，1代表是代谢位点；将分子结构数据与标签值以设定格式文件进行存储，该设定格式文件信息包括：a分子中每一个原子的类型与坐标，b原子与原子连接键的类型，c标签值，分子中原子是否为潜在代谢位点。

在优选的实施例中，所述原子类型的矢量化步骤中，将原子的质子数作为原子类型判断，将分子中的原子转化为矢量表示，采用0,1组成的行向量表示，设定原子的矢量长度为n,若原子的质子数为x则该原子对应的矢量的第x为1，其余n-1位为0，将分子中所有原子依照原子排序矢量化成为矩阵M，矩阵M为m*n矩阵，若分子中的原子数小于m则使用全0矢量进行补足；所述矩阵D中的每一个元素$D_{i,j}$表示为第i个原子与第j个原子的距离，该矩阵D中所有对角线元素全为0。

在优选的实施例中，所述原子的矢量的长度为78，若原子的质子数为m则该原子对应的矢量的第m为1，其余77位为0，矩阵M为100*78，包括100个行向量与78个列向量，若分子中的原子数小于100，则使用全0矢量进行补足；矩阵D为100*100矩阵，包括100个行向量和100个列向量。

在优选的实施例中，所述构建模型中，图卷积神经网络模型的超参数包括：图卷积半径r为3-5；图卷积神经网络模型训练过程中的超参数包括：mini_batch大小为32-128，学习率为0.0003-0.001；学习率为0.0003-0.001，卷积过程中的激发函数采用的sigmoid函数，损失函数采用对于原子是否是代谢位点的二分类交叉熵；所述更新顶点信息包括：中心原子的信息、以中心原子为中心设定半径范围内临近的化学键信息以及临近的原子信息；

经过设定层图卷积之后的中心原子的信息如下：

第一次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为1的所有临近的化学键的信息，c.半径为1的所有临近的原子信息；

第二次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为2的所有临近的化学键的信息，c.半径为2的所有临近的原子信息；

第三次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为3的所有临近的化学键的信息，c.半径为3的所有临近的原子信息；……以此类推。

在优选的实施例中，所述图卷积神经网络模型为空间域图卷积神经网络模型，空间域图卷积分解包括：消息传递过程，状态更新过程，所述消息传递过程包括：以节点为中心，将节点周围的原子信息与化学键的信息汇聚在一起；所述状态更新过程包括：依据中心节点的信息和基于消息传递过程汇聚的信息，综合起来更新中心节点的信息，所述空间域图神经网络模型输入作为节点信息矩阵的原子类型的矢量化矩阵M、距离的矢量化矩阵D，M矩阵的第j个原子的矢量化信息用c

所述消息传递过程中的函数形式：

其中，W、b为图卷积神经网络中需要训练的参数，v

所述状态更新过程中的函数形式，第j个节点被更新的公式

该公式表示了中心节点的状态更新：以中心节点本身的信息(c

针对于计算某个原子是不是代谢位点的概率，采用一层全连接神经网络计算节点原子是代谢位点的概率：

Probability＝sigmoid(c

在优选的实施例中，所述预测代谢位点步骤中在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名前三的原子预测为代谢位点。

一种预测小分子CYP代谢酶的代谢位点的系统，包括：

整理数据模块：收集与代谢位点相关的分子结构数据，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集；

结构转化模块：将二维分子结构转化为三维分子结构；

原子类型的矢量化模块：将分子三维结构中的原子类型的矢量化，将分子中所有原子的矢量依据原子排序矢量化形成矩阵M；

距离的矢量化模块：将分子三维结构中距离的矢量化，计算分子中两两原子距离，形成矩阵D；

构建模型模块：在图卷积神经网络模型中，输入原子类型的矢量化与距离的矢量化，将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率；

预测代谢位点模块：在不同CYP代谢酶亚型上预测代谢位点原子，在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名靠前的原子预测为代谢位点。

在优选的实施例中，所述整理数据模块中，对分子结构进行去盐、去溶剂处理，整理成为标准的模型输入格式，若分子设置有标签则在分子的标签值上，将原子转化为0,1标签值，0代表不是代谢位点，1代表是代谢位点；将分子结构数据与标签值以设定格式文件进行存储，该设定格式文件信息包括：a分子中每一个原子的类型与坐标，b原子与原子连接键的类型，c标签值，分子中原子是否为潜在代谢位点；

所述原子类型的矢量化模块中，将原子的质子数作为原子类型判断，将分子中的原子转化为矢量表示，采用0,1组成的行向量表示，设定原子的矢量长度为n,若原子的质子数为x则该原子对应的矢量的第x为1，其余n-1位为0，将分子中所有原子依照原子排序矢量化成为矩阵M，矩阵M为m*n矩阵，若分子中的原子数小于m则使用全0矢量进行补足；

所述距离的矢量化模块中，所述矩阵D中的每一个元素$D_{i,j}$表示为第i个原子与第j个原子的距离，该矩阵D中所有对角线元素全为0；

所述构建模型模块中，图卷积神经网络模型的超参数包括：图卷积半径r为3-5；图卷积神经网络模型训练过程中的超参数包括：mini_batch大小为32-128，学习率为0.0003-0.001；卷积过程中的激发函数采用的sigmoid函数，损失函数采用对于原子是否是代谢位点的二分类交叉熵；所述更新顶点信息包括：中心原子的信息、以中心原子为中心设定半径范围内临近的化学键信息以及临近的原子信息；

经过设定层图卷积之后的中心原子的信息如下：

第一次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为1的所有临近的化学键的信息，c.半径为1的所有临近的原子信息；

第二次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为2的所有临近的化学键的信息，c.半径为2的所有临近的原子信息；

第三次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为3的所有临近的化学键的信息，c.半径为3的所有临近的原子信息；……以此类推。

上述预测小分子CYP代谢酶的代谢位点的方法及系统，重点关注在药物分子在体内代谢过程中，分子与不同亚型的CYP代谢发生化学反应的代谢位点预测，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集，在不同CYP代谢酶亚型上预测代谢位点原子提高预测效果，通过在图卷积神经网络模型中，输入原子类型的矢量化与距离的矢量化，将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率，进一步提高预测效果，该预测结果可以帮助药物化学家设计或者优化分子结构。

附图说明

图1为本发明一实施例的预测小分子CYP代谢酶的代谢位点的方法的部分流程图；

图2为本发明一实施例的分子在CYP代谢酶亚型1A2上的预测结果示意；

图3为本发明一实施例的分子在CYP代谢酶亚型2A6上的预测结果示意；

图4为本发明一实施例的分子在CYP代谢酶亚型2B6上的预测结果示意；

图5为本发明一实施例的分子在CYP代谢酶亚型2C19上的预测结果示意；

图6为本发明一实施例的分子在CYP代谢酶亚型2C8上的预测结果示意；

图7为本发明一实施例的分子在CYP代谢酶亚型2C9上的预测结果示意；

图8为本发明一实施例的分子在CYP代谢酶亚型2D6上的预测结果示意；

图9为本发明一实施例的分子在CYP代谢酶亚型2E1上的预测结果示意；

图10为本发明一实施例的分子在CYP代谢酶亚型3A4上的预测结果示意。

具体实施方式

如图1所示，本发明一实施例的预测小分子CYP代谢酶的代谢位点的方法，包括：

步骤S101，整理数据：收集与代谢位点相关的分子结构数据，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集；

步骤S103，结构转化：将二维分子结构转化为三维分子结构；

步骤S105，原子类型的矢量化：将分子三维结构中的原子类型矢量化，将分子中所有原子的矢量依据原子排序矢量化形成矩阵M；

步骤S107，距离的矢量化：将分子三维结构中的距离矢量化，计算分子中两两原子距离，形成矩阵D；

步骤S109，构建模型：在图卷积神经网络模型中，输入原子类型的矢量化与距离的矢量化，将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率；

步骤S111，预测代谢位点：在不同CYP代谢酶亚型上预测代谢位点原子，在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名靠前的原子预测为代谢位点。

进一步，本实施例的预测小分子CYP代谢酶的代谢位点的方法，还包括：添加标签：根据预测结果对原子标记0,1标签值，0代表不是代谢位点，1代表是代谢位点。

本实施例中，可以从drugbank,ChemBL等公开数据集中收集与代谢位点相关的分子结构数据。

进一步，本实施例的整理数据步骤中，基于不同的代谢酶亚型将收集来的数据整理成为9种不同亚型的数据集：CYP 1A2、CYP 2A6、CYP 2B6、CYP 2C19、CYP 2C8、CYP 2C9、CYP 2D6、CYP 2E1、CYP 3A4。

由于从公开数据集中收集的数据，数据格式不一致，所以需要对于分子结构进行去盐，去溶剂的等处理，整理成为一个标准的模型输入格式。

若分子设置有标签则在分子的标签值上，将原子转化为0,1标签值，0代表不是代谢位点，1代表是代谢位点；将分子结构数据与标签值以设定格式文件进行存储，该设定格式文件信息包括：a分子中每一个原子的类型与坐标，b原子与原子连接键的类型，c标签值，分子中原子是否为潜在代谢位点。

优选的，本实施例中将一个样本中的数据与标签值以一个sdf文件的形式存储。Sdf包含的信息包括：a.分子中每一个原子的类型与坐标；b.原子与原子连接键的类型；c.标签值，分子中那个原子是潜在的代谢位点。最后以不同代谢酶亚型为一个文件夹包含这个sdf文件。模型的输入为a与b，模型的输出标签值为c。

步骤S103，结构转化：将二维分子结构转化为三维分子结构，可以使用公开的化学信息学软件工具包rdkit，采用rdkit中的mmff94二维分子结构转化为三维分子结构。

进一步，优选的，二维分子结构转化为三维分子结构的工具包可以采用的是rdkit，1.先定义好rdkit中的力场，采用这样的函数prop＝AllChem.MMFFGetMoleculeProperties(m,mmffVariant＝"MMFF94s"))定义“MMFF94s”力场，2.使用“MMFF94s”力场来随机初始化分子的三维坐标，rdkit.chem.AllChem.MMFFGetMoleculeForceField(m,prop,confId＝id)，3.因为随机的三维坐标可能因为键长或者不合理结构，需要优化一下初始化的分子三维结构，使用“MMFF94s”力场优化分子结构的三维构象。采用rdkit中的函数ff.Minimize()来优化三维构象。

进一步，本实施例的原子类型的矢量化步骤中，将原子的质子数作为原子类型判断，将分子中的原子转化为矢量表示，采用0,1组成的行向量表示，设定原子的矢量长度为n,若原子的质子数为x则该原子对应的矢量的第x为1，其余n-1位为0，将分子中所有原子依照原子排序矢量化成为矩阵M，矩阵M为m*n矩阵，若分子中的原子数小于m则使用全0矢量进行补足。

进一步，本实施例的原子的矢量的长度为78，若原子的质子数为m则该原子对应的矢量的第m位为1，其余77位为0，矩阵M为100*78，包括100个行向量与78个列向量，若分子中的原子数小于100，则使用全0矢量进行补足。

具体的，对于样本分子中每一个原子的矢量化，使用原子的质子数作为重要的原子类型的判断，比如氧(O)为8，碳(C)为6等等，统计所有分子中所有的原子类型，因为收集的分子都是有机分子，没有大于78质子数的原子。然后所有不同类型的原子依据化学周期表排列：[H,He,Li,Be,B,C,N,O…]，一共排列长度为78。对于一个分子也有一个依据分子结构对于分子内原子的排序，如：C,C,O,C,C,N。将分子中的第一个C原子转化为一个one-hot矢量，如[0,0,0,0,0,1,0,0…]。该矢量长度为78，该矢量只有第6位为1，其他位置都为0。依据此规则将分子中的第二个C原子也整理成[0,0,0,0,0,1,0,0…]，同理将分子中第二个O原子也整理成[0,0,0,0,0,0,0,1…]的矢量，以此类推。最后将一个分子中所有原子类型的矢量依据原子排序矢量化成为一个矩阵M，矩阵M的大小为100*78。矩阵M第一个维度是分子中原子的个数，第二个维度是每一个原子矢量化的长度。如果有些分子的原子数不足100，将采用使用全0的矢量补足矩阵M。

步骤S107，距离的矢量化步骤中，矩阵D中的每一个元素$D_{i,j}$表示为第i个原子与第j个原子的距离，该矩阵D中所有对角线元素全为0。

进一步，优选的，矩阵D为100*100矩阵，包括100个行向量和100个列向量。

进一步，本实施例的构建模型中，图卷积神经网络模型的超参数包括：图卷积半径r为3-5、mini_batch大小为32-128，学习率为0.0003-0.001，卷积过程中的激发函数采用的sigmoid函数，损失函数采用对于原子是否是代谢位点的二分类交叉熵。

进一步，本实施例中，在图卷积神经网络模型的输入为原子类型的矢量化与距离的矢量化。模型的输出为每一个原子是代谢位点的概率值。模型将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量化作为图卷积神经网络中边信息矢量化矩阵，依据图卷积神经网络中规则，来更新顶点信息。优选的，本实施例的图神经网络模型超参数包括：图卷积半径r为5。卷积过程中的激发函数采用的sigmoid函数。损失函数采用对于原子是否是代谢位点的二分类交叉熵。本实施例的图神经网络模型训练过程中的超参数包括：mini_batch大小为32，学习率为0.0003。

进一步，本实施例的更新顶点信息包括：中心原子的信息、以中心原子为中心设定半径范围内临近的化学键信息以及临近的原子信息；

经过设定层图卷积之后的中心原子的信息如下：

第一次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为1的所有临近的化学键的信息，c.半径为1的所有临近的原子信息；

第二次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为2的所有临近的化学键的信息，c.半径为2的所有临近的原子信息；

第三次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为3的所有临近的化学键的信息，c.半径为3的所有临近的原子信息；……以此类推。

进一步，本实施例的图卷积神经网络模型为空间域图卷积神经网络模型，空间域图卷积分解包括：消息传递过程，状态更新过程，所述消息传递过程包括：以节点为中心，将节点周围的原子信息与化学键的信息汇聚在一起；所述状态更新过程包括：依据中心节点的信息和基于消息传递过程汇聚的信息，综合起来更新中心节点的信息，所述空间域图神经网络模型输入作为节点信息矩阵的原子类型的矢量化矩阵M、距离的矢量化矩阵D，M矩阵的第j个原子的矢量化信息用c

消息传递过程中的函数形式：

其中，W、b为图卷积神经网络中需要训练的参数，即W

状态更新过程中的函数形式，第j个节点被更新的公式

该公式表示了中心节点的状态更新：以中心节点本身的信息(c

针对于计算某个原子是不是代谢位点的概率，采用一层全连接神经网络计算节点原子是代谢位点的概率：

Probability＝sigmoid(c

进一步，预测代谢位点步骤中在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名前三的原子预测为代谢位点。

参见图2-图10，在9个不同CYP代谢酶亚型上预测最有可能的代谢位点原子。

图2为分子中在CYP代谢酶亚型1A2上的预测结果。图2的左上角是sdf文件中原子排序的分子结构信息，图2的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图2右侧是具体每一个原子预测为代谢位点概率的具体数值。

图3为分子中在CYP代谢酶亚型2A6上的预测结果。图3的左上角是sdf文件中原子排序的分子结构信息，图3的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图3右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图4为分子中在CYP代谢酶亚型2B6上的预测结果。图4的左上角是sdf文件中原子排序的分子结构信息，图4的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图4右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图5为分子中在CYP代谢酶亚型2C19上的预测结果。图5的左上角是sdf文件中原子排序的分子结构信息，图5的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图5右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图6为分子中在CYP代谢酶亚型2C8上的预测结果。图6的左上角是sdf文件中原子排序的分子结构信息，图6的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图6右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图7为分子中在CYP代谢酶亚型2C9上的预测结果。图7的左上角是sdf文件中原子排序的分子结构信息，图7的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图7右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图8为分子中在CYP代谢酶亚型2D6上的预测结果。图8的左上角是sdf文件中原子排序的分子结构信息，图8的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图8右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图9为分子中在CYP代谢酶亚型2E1上的预测结果。图9的左上角是sdf文件中原子排序的分子结构信息，图9的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图9右侧是具体每一个原子预测为代谢位点的概率的具体数值。

图10为分子中在CYP代谢酶亚型3A4上的预测结果。图10的左上角是sdf文件中原子排序的分子结构信息，图10的左下角是预测的结果中，标注为阴影圆点处的原子最有可能的前3代谢位点。图10右侧是具体每一个原子预测为代谢位点的概率的具体数值。

本发明一实施例的预测小分子CYP代谢酶的代谢位点的系统，包括：

整理数据模块：收集与代谢位点相关的分子结构数据，基于不同代谢酶亚型将收集的数据整理成为不同亚型的数据集；

结构转化模块：将二维分子结构转化为三维分子结构；

原子类型的矢量化模块：将分子三维结构中的原子类型矢量化，将分子中所有原子的矢量依据原子排序矢量化形成矩阵M；

距离的矢量化模块：将分子三维结构中的距离矢量化，计算分子中两两原子距离，形成矩阵D；

构建模型模块：在图卷积神经网络模型中，输入原子类型的矢量化与距离的矢量化，将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量作为图卷积神经网络中边信息矢量化矩阵，根据图卷积神经网络模型更新顶点信息，计算原子为代谢位点的概率；

预测代谢位点模块：在不同CYP代谢酶亚型上预测代谢位点原子，在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名靠前的原子预测为代谢位点。

进一步，本实施例的预测小分子CYP代谢酶的代谢位点的系统，还包括：添加标签模块：根据预测结果对原子标记0,1标签值，0代表不是代谢位点，1代表是代谢位点。

本实施例中，可以从drugbank,ChemBL等公开数据集中收集与代谢位点相关的分子结构数据。

进一步，本实施例的整理数据步骤中，基于不同的代谢酶亚型将收集来的数据整理成为9种不同亚型的数据集：CYP 1A2、CYP 2A6、CYP 2B6、CYP 2C19、CYP 2C8、CYP 2C9、CYP 2D6、CYP 2E1、CYP 3A4。

由于从公开数据集中收集的数据，数据格式不一致，所以需要对于分子结构进行去盐，去溶剂的等处理，整理成为一个标准的模型输入格式。

若分子设置有标签则在分子的标签值上，将原子转化为0,1标签值，0代表不是代谢位点，1代表是代谢位点；将分子结构数据与标签值以设定格式文件进行存储，该设定格式文件信息包括：a分子中每一个原子的类型与坐标，b原子与原子连接键的类型，c标签值，分子中原子是否为潜在代谢位点。

优选的，本实施例中将一个样本中的数据与标签值以一个sdf文件的形式存储。Sdf包含的信息包括：a.分子中每一个原子的类型与坐标；b.原子与原子连接键的类型；c.标签值，分子中那个原子是潜在的代谢位点。最后以不同代谢酶亚型为一个文件夹包含这个sdf文件。模型的输入为a与b，模型的输出标签值为c。

本实施例的结构转化模块：将二维分子结构转化为三维分子结构，可以使用公开的化学信息学软件工具包rdkit，采用rdkit中的mmff94二维分子结构转化为三维分子结构。

进一步，优选的，二维分子结构转化为三维分子结构的工具包可以采用的是rdkit，1.先定义好rdkit中的力场，采用这样的函数prop＝AllChem.MMFFGetMoleculeProperties(m,mmffVariant＝"MMFF94s"))定义“MMFF94s”力场，2.使用“MMFF94s”力场来随机初始化分子的三维坐标，rdkit.chem.AllChem.MMFFGetMoleculeForceField(m,prop,confId＝id)，3.因为随机的三维坐标可能因为键长或者不合理结构，需要优化一下初始化的分子三维结构，使用“MMFF94s”力场优化分子结构的三维构象。采用rdkit中的函数ff.Minimize()来优化三维构象。

进一步，本实施例的原子类型的矢量化模块中，将原子的质子数作为原子类型判断，将分子中的原子转化为矢量表示，采用0,1组成的行向量表示，设定原子的矢量长度为n,若原子的质子数为x则该原子对应的矢量的第x为1，其余n-1位为0，将分子中所有原子依照原子排序矢量化成为矩阵M，矩阵M为m*n矩阵，若分子中的原子数小于m则使用全0矢量进行补足。

进一步，本实施例的原子的矢量的长度为78，若原子的质子数为m则该原子对应的矢量的第m位为1，其余77位为0，矩阵M为100*78，包括100个行向量与78个列向量，若分子中的原子数小于100，则使用全0矢量进行补足。

具体的，对于样本分子中每一个原子的矢量化，使用原子的质子数作为重要的原子类型的判断，比如氧(O)为8，碳(C)为6等等，统计所有分子中所有的原子类型，因为收集的分子都是有机分子，没有大于78质子数的原子。然后所有不同类型的原子依据化学周期表排列：[H,He,Li,Be,B,C,N,O…]，一共排列长度为78。对于一个分子也有一个依据分子结构对于分子内原子的排序，如：C,C,O,C,C,N。将分子中的第一个C原子转化为一个one-hot矢量，如[0,0,0,0,0,1,0,0…]。该矢量长度为78，该矢量只有第6位为1，其他位置都为0。依据此规则将分子中的第二个C原子也整理成[0,0,0,0,0,1,0,0…]，同理将分子中第二个O原子也整理成[0,0,0,0,0,0,0,1…]的矢量，以此类推。最后将一个分子中所有原子类型的矢量依据原子排序矢量化成为一个矩阵M，矩阵M的大小为100*78。矩阵M第一个维度是分子中原子的个数，第二个维度是每一个原子矢量化的长度。如果有些分子的原子数不足100，将采用使用全0的矢量补足矩阵M。

距离的矢量化模块中，矩阵D中的每一个元素$D_{i,j}$表示为第i个原子与第j个原子的距离，该矩阵D中所有对角线元素全为0。

进一步，优选的，矩阵D为100*100矩阵，包括100个行向量和100个列向量。

进一步，本实施例的构建模型中，图卷积神经网络模型的超参数包括：图卷积半径r为3-5、mini_batch大小为32-128，学习率为0.0003-0.001，卷积过程中的激发函数采用的sigmoid函数，损失函数采用对于原子是否是代谢位点的二分类交叉熵。

进一步，本实施例中，在图卷积神经网络模型的输入为原子类型的矢量化与距离的矢量化。模型的输出为每一个原子是代谢位点的概率值。模型将原子类型的矢量作为图卷积神经网络中顶点信息矢量化矩阵，将距离的矢量化作为图卷积神经网络中边信息矢量化矩阵，依据图卷积神经网络中规则，来更新顶点信息。优选的，本实施例的图神经网络模型超参数包括：图卷积半径r为5。卷积过程中的激发函数采用的sigmoid函数。损失函数采用对于原子是否是代谢位点的二分类交叉熵。本实施例的图神经网络模型训练过程中的超参数包括：mini_batch大小为32，学习率为0.0003。

进一步，本实施例的更新顶点信息包括：中心原子的信息、以中心原子为中心设定半径范围内临近的化学键信息以及临近的原子信息；

经过设定层图卷积之后的中心原子的信息如下：

第一次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为1的所有临近的化学键的信息，c.半径为1的所有临近的原子信息；

第二次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为2的所有临近的化学键的信息，c.半径为2的所有临近的原子信息；

第三次卷积过后，更新的顶点信息包括：a.中心原子的信息，b.半径为3的所有临近的化学键的信息，c.半径为3的所有临近的原子信息；……以此类推。

进一步，本实施例的图卷积神经网络模型为空间域图卷积神经网络模型。空间域图卷积分解包括：消息传递过程，状态更新过程，所述消息传递过程包括：以节点为中心，将节点周围的原子信息与化学键的信息汇聚在一起；所述状态更新过程包括：依据中心节点的信息和基于消息传递过程汇聚的信息，综合起来更新中心节点的信息，所述空间域图神经网络模型输入作为节点信息矩阵的原子类型的矢量化矩阵M、距离的矢量化矩阵D，M矩阵的第j个原子的矢量化信息用c

消息传递过程中的函数形式：

其中，W、b为图卷积神经网络中需要训练的参数，即W

状态更新过程中的函数形式，第j个节点被更新的公式

该公式表示了中心节点的状态更新：以中心节点本身的信息(c

针对于计算某个原子是不是代谢位点的概率，采用一层全连接神经网络计算节点原子是代谢位点的概率：

Probability＝sigmoid(c

进一步，预测代谢位点步骤中在不同CYP代谢酶亚型上根据为代谢位点的概率取出排名前三的原子预测为代谢位点。

参见图2-图10，在9个不同CYP代谢酶亚型上预测最有可能的代谢位点原子。

本发明重点关注在药物分子在体内代谢过程中，分子与不同亚型的CYP代谢发生化学反应的代谢位点预测。该预测结果可以帮助药物化学家设计或者优化分子结构。

以上述依据本申请的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项申请技术思想的范围内，进行多样的变更以及修改。本项申请的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。