基于可变半径高斯球的分子体积计算方法以及用于药物虚拟筛选的分子三维相似度打分方法与流程

1.本发明涉及计算机辅助药物研发技术领域，具体涉及一种基于可变半径高斯球的分子体积计算方法以及用于药物虚拟筛选的分子三维相似度打分方法。


背景技术：

2.药物研发具有投入大、风险高、周期长的特点，通常一个药物研发周期在10年以上，研发投入在数亿美金，并且呈现逐年上升趋势。药物筛选是药物发现的关键环节，而高通量药物虚拟筛选能大大降低筛选时间和成本，对于加速药物研发具有重要意义。
3.在药物虚拟筛选中，目前常用到基于分子三维相似度打分方法进行分子的筛选排序。这类打分方法包含分子形状的相似性和药效团的相似性(通常是选用数据库分子与模版分子或药效团模型的相似度打分)，并以此通过一个简单的加权函数构成一个相似度的综合打分，这个打分函数的有效性决定了筛选的有效性和计算速度。然而，以上采用的打分方法存在的显著问题是准确性较差，使得药物虚拟筛选存在较高的假阳性或假阴性；另一方面，由于分子体积与分子形状相关，传统的药物分子三维体积的计算方法是通过一组相互叠合的高斯球模型来计算的，例如公开号为cn102799779a、专利名称为“一种药物分子形状的比较方法”的发明专利，其技术方案中计算分子三维体积的方法是：根据分子中各原子的类型得到相应的范德华半径，将分子三维结构信息转换为一组代表该分子中各原子的高斯球，每个高斯球的半径的取值等于相应原子的范德华半径，且每个高斯球的位置与相应原子的坐标相同；在此基础上计算分子三维体积，然后对两分子形状的相似度进行比较。然而，以上计算方法存在的缺陷是，分子中表示原子体积的高斯球半径是取相应原子的范德华半径，这样做的结果并没有消除由于高斯球相互重叠带来的误差，这也使得药物虚拟筛选的精度受到影响。


技术实现要素：

4.本发明的目的之一在于针对现有技术中的不足，而提供一种筛选精度高的基于可变半径高斯球的分子体积计算方法。
5.本发明的目的之二在于针对现有技术中的不足，而提供一种用于药物虚拟筛选的分子三维相似度打分方法，该方法在药物筛选准确率上有显著的改进，同时保持着筛选高通量的计算速度。
6.本发明的目的通过以下技术方案实现：
7.提供一种基于可变半径高斯球的分子体积计算方法，该方法包括以下步骤：
8.读入分子的三维结构信息，所述三维结构信息包括分子中每个原子的类型及其坐标数值；
9.根据分子中各原子的类型得到原子的范德华半径，将三维结构信息转换为一组代表分子中各原子的高斯球，每个高斯球的半径不等于原子的范德华半径，而是小于其范德
华半径，高斯球半径的计算公式如下：
[0010][0011]
其中：对于原子i，其范德华半径为ri，体积为vi，该原子与其它原子总的重叠体积为oi；α为可调参数，取值在0～1.0之间，通过虚拟筛选的测试数据调整到最佳分类效果；通过以上公式计算得到调整后的高斯球半径r
i
，高斯球的位置坐标与原子的坐标相同；
[0012]
计算分子中高斯球组的叠合体积，其中第ij个高斯球组包括第i个原子对应的高斯球和第j个原子对应的高斯球，第ij个高斯球组的叠合体积为v
ij
；计算分子的体积即为其自身的叠合体积n是分子中原子的总数。
[0013]
本发明还提供一种用于药物虚拟筛选的分子三维相似度打分方法，包括以下步骤：
[0014]
步骤一、基于可变半径高斯球的分子体积计算
[0015]
该步骤是采用权利要求1的计算方法；
[0016]
步骤二、基于可变半径高斯球，计算用于相似性比较的两分子的各类特征参数
[0017]
分别读取用于相似性比较的两分子的拓扑结构和三维结构信息，获取各类特征参数，所述特征参数包括：两分子的原子个数差异(f1)；两分子的可旋转化学键个数(f2)；两分子的体积差异(f3)；两分子的形状相似度(f4)；两分子氢键受体的相似度(f5)；两分子氢键供体的相似度(f6)；两分子芳香环的相似度(f7)；两分子的疏水中心相似度(f8)；两分子的正电基团相似度(f9)；以及两分子的负电基团相似度(f10)；其中：
[0018]
f1的计算方式是通过读入两分子各自的拓扑结构信息，然后取两分子原子总数差值的绝对值；
[0019]
f2的计算方式是在f1计算方式的基础上，判断每个化学键是否为可旋转键，得到两分子各自的可旋转键总数，然后取两分子可旋转键总数差值的绝对值；
[0020]
f3的计算方式是在f1计算方式的基础上，根据两分子中各原子的类型得到原子的范德华半径，采用步骤一的方法，计算两分子各自自身的叠合体积为n是分子中原子的总数；然后取两分子自身叠合体积的差值的绝对值；
[0021]
f4的计算方式是在f3计算方式的基础上，计算两分子在多种叠合情况下的分子间叠合体积其中v
ij
为第一分子中的第i个原子与第二个分子中的第j个原子的叠合体积，n是第一个分子中原子的总数，m是第二个分子中原子的总数，选择其中的最大值作为最大分子间体积；计算两分子的形状相似度其中v
a
是第一个分子的自身叠合体积，v
b
是第二个分子的自身叠合体积；
[0022]
f5的计算方式是在f1计算方式的基础上，找出两分子中氢键受体位置；计算两分子中各自氢键受体的叠合体积其中f
ij
为第i个氢键受体与第j个氢键受体之间
的叠合体积；计算两分子在多种叠合情况下的分子间氢键受体的叠合体积其中f
ij
为第一分子中的第i个氢键受体与第二个分子中的第j个氢键受体的叠合体积，n是第一个分子中氢键受体的总数，m是第二个分子中氢键受体的总数，选择其中的最大值作为最大分子间氢键受体的叠合体积；计算两分子的氢键受体相似度其中p
a
是第一个分子中氢键受体的自身叠合体积，p
b
是第二个分子中氢键受体的自身叠合体积；
[0023]
f6的计算方式与f5的计算方式相同，只需要将两分子中的氢键受体替换为氢键供体；
[0024]
f7的计算方式与f5的计算方式相同，只需要将氢键受体替换为芳香环；
[0025]
f8的计算方式与f5的计算方式相同，只需要将氢键受体替换为疏水中心；
[0026]
f9的计算方式与f5的计算方式相同，只需要将氢键受体替换为正电基团；
[0027]
f10的计算方式与f5的计算方式相同，只需要将氢键受体替换为负电基团；
[0028]
步骤三、训练深度学习模型
[0029]
采用dud
‑
e数据集，数据集中有102个生物靶点信息，每个靶点都有对应的活性分子集与decoy分子集，对每个靶点的数据做如下处理：
[0030]
选取每个靶点的活性分子集中的晶体结构分子作为模板分子，分别与该靶点的活性分子集的其他分子，以及decoy分子集的所有分子，均按照步骤二中的计算方式，计算模板分子与其他分子中两两分子的f1～f10特征参数，使每个靶点都计算得到一套特征参数数据；
[0031]
使用深度学习方法进行建模，将上述计算得到的每个靶点的特征参数数据作为输入数据，而分子的活性与否作为二分类的目标函数，模型优化的方向是使得所有靶点在分子活性预测的误差最小化，从而使auc值的平均值最大；完成训练后，即得到最终的深度学习模型；
[0032]
步骤四、深度学习模型的外部验证
[0033]
采用muv数据集验证深度学习模型的泛化能力，选取muv数据集中的10个生物靶点信息，每个靶点都有对应的活性分子集与decoy分子集；选取每个靶点的活性分子集中的晶体结构分子作为模板分子，分别与该靶点的活性分子集的其他分子，以及decoy分子集的所有分子，均按照步骤二中的计算方式，计算模板分子与其他分子中两两分子的f1～f10特征参数；将这些特征参数输入已经训练完成的深度学习模型，计算得到每个靶点虚拟筛选的auc值。
[0034]
上述技术方案中，步骤二中，所述三维结构信息包括分子中原子总数、化学键总数、每个原子的类型及其坐标数值。
[0035]
上述技术方案中，步骤三中，采用类5倍交叉验证的方式，对每个靶点计算虚拟筛选的auc值。
[0036]
本发明的有益效果：
[0037]
本发明的一种用于药物虚拟筛选的分子三维相似度打分方法，包括步骤一、基于可变半径高斯球的分子体积计算方法；步骤二、基于可变半径高斯球，计算用于相似性比较的两分子的各类特征参数，该特征参数包括两分子的原子个数差异、可旋转化学键个数、体
积差异、形状相似度、氢键受体的相似度、氢键供体的相似度、芳香环的相似度、疏水中心相似度、正电基团以及负电基团相似度；步骤三、训练深度学习模型；以及步骤四、深度学习模型的外部验证，采用muv数据集验证该模型的泛化能力。与现有技术相比，本发明具有以下优点：
[0038]
(1)在分子体积的计算方法中，本发明改进了高斯球半径的取值方法，即高斯球半径不等于相应原子的范德华半径，而是小于其范德华半径，且改变的比例与该原子所在的化学环境相关，计算得到调整后的高斯球半径以消除由于高斯球相互重叠带来的误差；以上对高斯球半径取值方法做出的调整是根据该原子所对应高斯球与其它原子所对应高斯球的重叠程度来确定的，重叠度愈高，相应半径的减少则越大，而半径减少的程度是重叠度的一个单调递减函数。通过药物虚拟筛选测试结果显示，本发明有效提升了药物虚拟筛选的精度；
[0039]
(2)本发明选取了多维度相似性作为特征参数(包括两分子的原子个数差异、可旋转化学键个数、体积差异、形状相似度、氢键受体的相似度、氢键供体的相似度、芳香环的相似度、疏水中心相似度、正电基团以及负电基团相似度)，经深度神经网络的变换，给出最终的综合打分即auc值。auc值评价指标是评价药物虚拟筛选方法准确性的常用标准，auc值取值范围在0～1之间，数值越接近1，表明筛选方法越准确。经实验证明，本发明相比现有技术在药物筛选准确率上有显著的改进，同时保持着筛选高通量的计算速度。因此，采用本发明的方法，能够尽可能多地将具有潜在生物活性的分子从含有大量分子的数据库中筛选出来，而筛选方法越准确，就越容易找到具有潜在活性的分子，这样就能尽可能地不漏掉任何潜在活性分子，从而解决了假阳性或假阴性的问题。
具体实施方式
[0040]
下面结合具体实施例对本发明作进一步的详细说明，但并不构成对本发明的任和限制。
[0041]
实施例1：基于可变半径高斯球的分子体积计算方法
[0042]
本实施例的基于可变半径高斯球的分子体积计算方法，包括以下步骤：
[0043]
读入分子的三维结构信息，该三维结构信息包括分子中每个原子的类型及其坐标数值。
[0044]
根据分子中各原子的类型得到原子的范德华半径，将三维结构信息转换为一组代表分子中各原子的高斯球，即每个原子用一个高斯球表示，每个高斯球的半径不等于原子的范德华半径，而是小于其范德华半径，调整后的高斯球半径的计算公式如下：
[0045][0046]
其中：对于原子i，其范德华半径为ri，体积为vi，该原子与其它原子总的重叠体积为oi；α为可调参数，具体取值在0～1.0之间，通过虚拟筛选的测试数据调整到最佳分类效果；通过以上公式计算得到调整后的高斯球半径r
i
，高斯球的位置坐标与原子的坐标相同；
[0047]
计算分子中高斯球组的叠合体积，其中第ij个高斯球组包括第i个原子对应的高斯球和第j个原子对应的高斯球，第ij个高斯球组的叠合体积为v
ij
；计算分子的体积即为其自身的叠合体积n是分子中原子的总数。
[0048]
本发明改进了高斯球半径的取值方法，即高斯球半径不等于相应原子的范德华半径，而是其范德华半径的变化，以上对高斯球半径取值方法做出的调整是根据该原子所对应高斯球与其它原子所对应高斯球的重叠程度来确定的，重叠度愈高，相应半径的减少则越大，而半径减少的程度是重叠度的一个单调递减函数。因此，该计算方法能够消除由于高斯球相互重叠带来的误差，提升药物虚拟筛选的精度。
[0049]
实施例2：一种用于药物虚拟筛选的分子三维相似度打分方法
[0050]
本实施例以dud
‑
e数据集中的靶点ada17为例，dud
‑
e数据集中共有102个生物靶点信息，每个靶点都有对应的活性分子集与decoy分子集。其中，ada17靶点的数据集中含有1,341个活性分子集和35,900个decoy分子集。分别选取晶体结构中的活性分子为模板分子(以下称为“分子a”)和活性分子集中的第一个分子(以下称为“分子b”)。
[0051]
一种用于药物虚拟筛选的分子三维相似度打分方法的具体实施方式，包括以下步骤：
[0052]
步骤一、基于可变半径高斯球的分子体积计算方法：
[0053]
该步骤与实施例1的完全相同。
[0054]
步骤二、基于可变半径高斯球，计算用于相似性比较的两分子的各类特征参数：
[0055]
读取分子a和分子b的拓扑结构和三维结构信息，该三维结构信息包括分子中原子总数、化学键总数、每个原子的类型及其坐标数值，然后按以下步骤，分别计算得到用于相似性比较的特征参数：
[0056]
步骤1、取分子a与b的原子总数差值的绝对值，获得第一个特征参数f1。
[0057]
步骤2、判断分子a与b中每个化学键是否为可旋转键，分别得到分子a与b的可旋转键总数，取两者总数差值的绝对值，获得第二个特征参数f2。
[0058]
步骤3、根据分子a与b中各原子的类型得到原子的范德华半径，每个原子用一个高斯球表示，采用步骤一的方法计算调整后的高斯球半径，高斯球的位置坐标与原子的坐标相同；然后计算分子a中任意两两原子对应的一组高斯球(以下称为高斯球组)的叠合体积，其中第ij个高斯球组包括分子a中第i个原子对应的高斯球和第j个原子对应的高斯球，第ij个高斯球组的叠合体积为v
ij
；计算分子a自身的叠合体积为n是分子a中原子的总数；使用同样的方法计算分子b自身的叠合体积为m是分子b中原子的总数；然后取分子a和b自身叠合体积的差值的绝对值，得到f3。
[0059]
步骤4、计算分子a与b在多种叠合情况下的分子间叠合体积其中v
ij
为分子a中的第i个原子与分子b中第j个原子的叠合体积，n是分子a中原子的总数，m是分子b中原子的总数，选择其中的最大值作为最大分子间体积；计算两分子的形状相似度
其中v
a
是分子a的自身叠合体积，v
b
是分子b的自身叠合体积(即步骤3中计算得到)。
[0060]
步骤5、找出分子a与b中所有氢键受体的位置；计算分子a中氢键受体的叠合体积其中f
ij
为分子a中第i个氢键受体与第j个氢键受体之间的叠合体积，n为分子a中氢键受体的总数；使用同样的方式计算分子b中氢键受体的叠合体积其中f
ij
为分子b中第i个氢键受体与第j个氢键受体之间的叠合体积，m为分子b中氢键受体的总数；计算分子a与b在多种叠合情况下的分子间氢键受体的叠合体积其中f
ij
为分子a中的第i个氢键受体与分子b中的第j个氢键受体的叠合体积，n是分子a中氢键受体的总数，m是分子b中氢键受体的总数，选择其中的最大值作为分子a与b之间氢键受体的叠合体积；计算分子a与b的氢键受体相似度其中p
a
是分子a中氢键受体的自身叠合体积，p
b
是分子b中氢键受体的自身叠合体积。
[0061]
步骤6、与步骤5的计算方式一致，只需要将氢键供体替换为氢键受体，即可得到分子a与b的氢键受体相似度f6。
[0062]
步骤7、与步骤5的计算方式一致，只需要将氢键供体替换为芳香环，即可得到分子a与b的芳香环相似度f7。
[0063]
步骤8、与步骤5的计算方式一致，只需要将氢键供体替换为疏水中心，即可得到分子a与b的疏水中心相似度f8。
[0064]
步骤9、与步骤5的计算方式一致，只需要将氢键供体替换为正电基团，即可得到分子a与b的正电基团相似度f9。
[0065]
步骤10、与步骤5的计算方式一致，只需要将氢键供体替换为负电基团，即可得到分子a与b的负电基团相似度f10。
[0066]
由此得到分子a与活性分子集中第一个分子即分子b的特征参数f1～f10，共10个。
[0067]
步骤三、训练深度学习模型：
[0068]
同样以dud
‑
e数据集中靶点ada17为例，采用与步骤二相同的计算方式，计算分子a与活性分子集中第二个分子(以下称为“分子c”)的特征参数，也得到10个相应的特征参数。
[0069]
以此类推，分别计算分子a与活性分子集中第三个、第四个、第五个.......第n个分子(n为自然数)，直到活性分子集中的所有活性分子都已计算，得到特征参数的数量为10*1341＝13,410个特征参数。
[0070]
接着，采用与步骤一相同的计算方式，计算分子a与decoy分子集中每个分子的特征参数，得到10*35900＝359,000个特征参数。至此，靶标ada17的所有特征参数计算完毕。
[0071]
然后，采用以上同样的计算方式，计算dud
‑
e数据集中其他101个靶点的特征参数集。至此，dud
‑
e数据集中的102个靶点，每个靶点都有一套特征参数数据。
[0072]
最后，将以上102套特征参数数据作为深度学习模型的输入特征数据，而分子的活
性与否作为二分类的目标函数，模型优化的方向是使得所有靶点在分子活性预测的误差最小化，从而使auc值的平均值最大。训练过程中，采用类5倍交叉验证的方式，对每个靶点计算虚拟筛选的auc值(见表1)，并取所有102个靶点auc值的平均值aucaver，模型优化的方向是使得aucaver的值最大。完成训练后，即可得到最终的深度学习模型。
[0073]
表1.dud
‑
e数据集采用传统方法和本实施例的方法计算得到的auc值
[0074]
[0075]
[0076][0077]
步骤四、深度学习模型的外部验证：
[0078]
采用muv数据集验证深度学习模型的泛化能力。以muv数据集中的靶点466为例，466靶点数据集中含有31个活性分子集和15000个decoy分子集。选取晶体结构中的活性分子为模板分子(以下称为分子a)。采用步骤二中的计算方式，得到分子a分别与活性分子集和decoy分子集中的每个分子的特征参数，总共得到10*(31 15000)＝150310个特征参数的数据集合。将这些特征参数输入已经训练完成的深度学习模型中，即可计算得到该靶点虚拟筛选的auc值(见表2)。
[0079]
表2.muv数据集采用传统方法和本实施例的方法计算得到的auc值
[0080][0081]
表1和表2中均以传统打分方法计算得到的auc值作为对比例，该传统打分方法指的是背景技术所提到的方法即选用数据库分子与模版分子或药效团模型的相似度打分，并以此通过一个简单的加权函数构成一个相似度的综合打分方法，其中药物分子三维体积的计算方法中，每个高斯球半径的取值等于相应原子的范德华半径，在此基础上计算得到分子三维体积，然后对两分子形状的相似度进行比较。
[0082]
药物虚拟筛选的auc值评价指标是评价筛选方法准确性的常用标准，auc值取值范围在0
‑
1之间，数值越接近1，表明筛选方法越准确。
[0083]
如表1和表2所示的实验结果证明，本发明的打分方法相比现有技术在药物筛选准确率上有显著的改进，同时保持着筛选高通量的计算速度。这样能够尽可能多地将具有潜在生物活性的分子从含有大量分子的数据库中筛选出来，而筛选方法越准确，就越容易找到具有潜在活性的分子，而不漏掉任何潜在活性分子，从而解决了假阳性或假阴性的问题。
[0084]
以上所举实施例为本发明的较佳实施方式，仅用来方便说明本发明，并非对本发明作任何形式上的限制，任何所属技术领域中具有通常知识者，若在不脱离本发明所提技术特征的范围内，利用本发明所揭示技术内容所作出局部更动或修饰的等效实施例，并且未脱离本发明的技术特征内容，均仍属于本发明技术特征的范围内。