CN114783545B - 基于gpu加速的分子对接方法和装置 - Google Patents
基于gpu加速的分子对接方法和装置 Download PDFInfo
- Publication number
- CN114783545B CN114783545B CN202210449861.2A CN202210449861A CN114783545B CN 114783545 B CN114783545 B CN 114783545B CN 202210449861 A CN202210449861 A CN 202210449861A CN 114783545 B CN114783545 B CN 114783545B
- Authority
- CN
- China
- Prior art keywords
- molecular
- conformations
- conformation
- random
- gpu
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000003032 molecular docking Methods 0.000 title claims abstract description 59
- 238000000034 method Methods 0.000 title claims abstract description 26
- 230000001133 acceleration Effects 0.000 title claims abstract description 22
- 239000003446 ligand Substances 0.000 claims abstract description 48
- 230000003993 interaction Effects 0.000 claims abstract description 18
- 238000010845 search algorithm Methods 0.000 claims abstract description 17
- 210000001503 joint Anatomy 0.000 claims abstract description 10
- 239000003814 drug Substances 0.000 claims abstract description 8
- 229940079593 drug Drugs 0.000 claims abstract description 8
- 230000006870 function Effects 0.000 claims description 64
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000003041 virtual screening Methods 0.000 claims description 9
- 238000007637 random forest analysis Methods 0.000 claims description 7
- 238000004891 communication Methods 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 4
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 claims description 4
- 239000010931 gold Substances 0.000 claims description 4
- 229910052737 gold Inorganic materials 0.000 claims description 4
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 abstract description 2
- 238000012827 research and development Methods 0.000 abstract description 2
- 125000004429 atom Chemical group 0.000 description 33
- 230000006872 improvement Effects 0.000 description 10
- 150000001875 compounds Chemical class 0.000 description 8
- 238000012216 screening Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 125000004433 nitrogen atom Chemical group N* 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 150000003384 small molecules Chemical class 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000009510 drug design Methods 0.000 description 1
- 238000009509 drug development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000003092 force field based scoring function Methods 0.000 description 1
- 238000000338 in vitro Methods 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002547 new drug Substances 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/50—Molecular design, e.g. of drugs
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F9/00—Arrangements for program control, e.g. control units
- G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
- G06F9/46—Multiprogramming arrangements
- G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
- G06F9/5005—Allocation of resources, e.g. of the central processing unit [CPU] to service a request
- G06F9/5027—Allocation of resources, e.g. of the central processing unit [CPU] to service a request the resource being a machine, e.g. CPUs, Servers, Terminals
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/70—Machine learning, data mining or chemometrics
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/90—Programming languages; Computing architectures; Database systems; Data warehousing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D10/00—Energy efficient computing, e.g. low power processors, power management or thermal management
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Computing Systems (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Software Systems (AREA)
- Chemical & Material Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Crystallography & Structural Chemistry (AREA)
- Physics & Mathematics (AREA)
- Databases & Information Systems (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Artificial Intelligence (AREA)
- General Engineering & Computer Science (AREA)
- Medicinal Chemistry (AREA)
- Pharmacology & Pharmacy (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于GPU加速的分子对接方法及装置，包括：CPU端进行数据准备以及设置OpenCL环境，其中数据准备包括对接盒子及对接盒子的参数准备，随后将数据传输到GPU端；GPU端计算对接盒子中格点处原子和受体的原子的相互作用力并生成网格缓存表；CPU端生成配体的多个随机分子构象并传到GPU端，GPU端利用蒙特卡罗的迭代局部搜索算法优化随机分子构象并打分，得到构象并传输至回CPU端；CPU端将构象进行优化和排序；依次对配体库中所有配体进行上述优化步骤；使用打分函数ΔvinaRF20对CPU端中的分子构象进行打分并排序。本发明通过ΔvinaRF20对构象重新打分并排序，提高了输出构象的准确性；提出了一种基于OpenCL的AutoDock Vina异构并行加速结构，降低了对接时间和药物研发成本。
Description
技术领域
本发明涉及一种基于GPU加速的分子对接方法和装置，属于硬件加速药物设计领域。
背景技术
新药研发通常需要针对特定的靶点测试一个规模庞大的配体库，可以通过分子对接软件预测蛋白质受体与小分子配体之间的结合模式及其结合亲和力，对结合后的化合物的结合亲和力进行排序，减少进行后续体内和体外测试的配体数量。
AutoDock Vina是分子对接领域最为常用的开源对接软件，在CASF-2016的评分函数对比评测中，AutoDock Vina是所有分子对接工具中对接能力最高的软件。AutoDockVina使用了基于梯度的优化器BFGS和基于蒙特卡罗的迭代局部搜索算法，有效地提升了搜索效率和化合物结合模式预测的准确性；AutoDock Vina使用CPU的多线程性，有效提升了AutoDock Vina的运行速度，减少分子对接所花费的时间。分子对接最关键的部分为打分函数，打分函数用于确定小分子配体在蛋白质受体上的结合位点和结合模式，并根据受体和配体的复杂结构估计它们之间的结合亲和力，筛选虚拟小分子配体库，以进一步确定苗头化合物。AutoDock Vina中的打分函数是基于力场的打分函数，主要计算了范德华力以及库仑力。大量研究表明，AutoDock Vina中的打分函数在对接和筛选能力测试中表现相对较好，但在打分能力测试中表现较差。
ΔvinaRF20是基于随机森林的打分函数。ΔvinaRF20扩大了训练集规模，训练集中增加弱结合亲和力的复合物结构和CSAR decoy set中的复合物结构，提高了模型的泛化性能；采用随机森林对AutoDock Vina的打分函数进行参数化修正，将AutoDock Vina中优秀的对接能力与随机森林的优势相结合，提升了预测分子对接软件生成构象的结合亲和力的准确性。在CASF-2013、CASF-2007、CASF-2016的所有包含打分、排序、对接和筛选能力的基准测试中，ΔvinaRF20表现得最为突出。CASF-2016评估了多个主流打分函数，ΔvinaRF20的打分能力、排序能力排名第一，对接能力和筛选能力测试排名第二，对接成功率高达90％。
在虚拟筛选时，化合物的规模非常重要。随着虚拟筛选规模的增加，所得构象的平均得分也更好，从而提高了筛选出具有更高结合亲和力的分子的几率，这会导致更高的真实命中率和实验结合亲和力。
由于类药小分子化合物的整个化学空间可望达到1060，研究人员对庞大数量的化合物数据库进行筛选时需要很多计算资源和时间，研究人员预计使用单个中央处理单元筛选100亿种化合物需要3000多年。AutoDock Vina用于真实情景下的药物虚拟筛选时，对接速度较慢、化合物结合模式预测的准确性有待提高。
有鉴于此，需要对现有的分子对接方法和装置进行改进，以解决上述问题。
发明内容
本发明的目的在于提供一种基于GPU加速的分子对接方法及装置，以解决现有分子对接方法及装置的打分能力差、对接速度慢等问题。
为实现上述目的，本发明提供了一种基于GPU加速的分子对接方法，用于将受体和配体库中的配体进行分子对接，包括如下步骤：
S1、CPU端进行数据准备以及设置OpenCL环境，其中，数据准备包括对接盒子及对接盒子的参数准备，并将准备完成的数据传输到GPU端；
S2、所述对接盒子中包括若干个格点，GPU端的内核函数一用于计算所述格点和所述受体中的原子之间的相互作用力，并生成网格缓存表；
S3、CPU端生成配体的多个随机分子构象并传到GPU端，GPU端的内核函数二利用蒙特卡罗的迭代局部搜索算法优化随机分子构象并打分，得到多个构象，并将多个构象传输至CPU端；
S4、CPU端将S3输出的多个构象进行优化和排序，得到分子构象；
S5、重复S3和S4，直至所述配体库中所有配体均完成S3和S4；
S6、使用ΔvinaRF20对CPU端中的所有分子构象进行打分并排序，保留排名靠前的若干个分子构象，实现药物的虚拟筛选。
作为本发明的进一步改进，S1中数据准备包括：
S11、读取文件，包括：
读取受体文件，所述受体文件的格式为.pdbqt；
读取对接盒子参数文件，所述对接盒子参数文件包括：
所述受体文件的路径；
所述配体库的路径，所述配体库包括若干个所述配体；
所述对接盒子的大小，分别为size_x、size_y、size_z；
所述对接盒子的中心位点的三维坐标,分别为center_x、center_y、center_z；
线程数的大小为thread；
S12、所述OpenCL环境的设置，包括：
查询平台，使用所述OpenCL中的clGetPaltformIDs()函数获取可用的OpenCL环境的数量，并初始化一个可用的OpenCL环境；
查询设备，使用所述OpenCL中的clGetDeviceIDs()函数获取平台中设备的数量，并初始化一个可用的设备；
创建上下文，所述上下文中包括内存、程序以及内核对象，使用clCreateContext()函数对初始化后的设备创建上下文，用于设备与程序以及设备与设备之间的通信；
创建命令队列，使用cl_queue_properties()的属性值创建命令队列，所述命令队列能够操作所述上下文中的多个对象；
创建程序对象，使用clCreateProgramWithSource()函数对OpenCL环境中的源代码创建程序对象，使用clBuildProgram()函数为指定的设备创建程序对象；
创建内核对象，使用clCreateKernel()函数创建两个内核对象；
设置内核参数，使用clSetKernelArg()函数对所述内核对象设置内核参数，便于执行所述内核对象；
执行内核，使用clEnqueueNDRangeKernel()函数调用所述命令队列，使得所述内核对象在设备进行排队；
S13、生成分子内能量表，用于计算所述分子构象的分子内能量；
S14、根据随机数种子生成随机数查找表。
作为本发明的进一步改进，S2具体包括：
将所述对接盒子划分成多个格子，所述格点位于所述格子的顶角处，划分后所述格子的大小为l*m*n，得到所述格点的量化坐标，所述内核函数一中包括与若干个所述格点相对应的若干个线程；
将所述量化坐标转换为受体中的真实坐标；
以所述格点为中心建立正方体，所述正方体内包含所述受体的原子，假设所述格点处存在原子，计算所述格点处原子与受体原子的相互作用力，将所述相互作用力进行存储，得到所述网格缓存表，并将所述网格缓存表存放在全局内存中，供所述GPU线程共享。
作为本发明的进一步改进，所述正方体的大小为所述格点与所述配体原子之间的距离小于/>时，计算所述格点处原子与所述受体原子的相互作用力，所述格点处原子的类型共有17种。
作为本发明的进一步改进，S3具体包括：
S31、随机分子构象的表示方式如下：
随机分子构象在搜索空间中的三维坐标，Position＝(a,b,c)；
分子构象的四元数，orientation＝(R0,R1,R2,R3)；
分子构象中可旋转键的扭转角度，torsion＝(t0,t1,t2)；
S32、每个随机分子构象对应一个线程，每个线程中均进行蒙特卡罗的迭代局部搜索算法，线程数为8000；
S33、将内核函数二中每个线程得到的分子构象均传到CPU端。
作为本发明的进一步改进，S32包括：
S321、对随机分子构象进行扰动；
S322、对随机分子构象进行BFGS，包括：
a、将海森矩阵初始化为单位矩阵；
b、对随机分子构象进行三维空间的转换，得到随机分子构象中每个原子的空间坐标；
c、使用打分函数计算随机分子构象的能量，随机分子构象的能量包含分子间能量与分子内能量，通过所述网格缓存表计算分子间能量，通过分子内能量表计算分子内能量；
d、计算构象中每个原子的空间坐标的导数，将其转换成随机分子构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数；
e、使用Armijo-Goldstein准则进行黄金搜索得到新的分子构象，初始化步长alpha＝1，计算新的分子构象能量并计算分子构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数，使用更新前后的分子构象的能量判断步长是否需要更新，若需要更新则alpha＝alpha/2，继续进行下一轮黄金搜索，否则结束搜索；
f、判断更新后分子构象的导数是否满足‖导数‖≥1e-5,若满足条件则更新海森矩阵；
g、判断e、f循环次数是否达到10次，若满足则输出为BFGS最终得到的分子构象；
S323、使用接受准则判断是否接收所述S322得到的分子构象；
S324、判断S322和S323的循环次数是否达到设定值，循环次数达到设定值后输出本线程的分子构象。
作为本发明的进一步改进，S4包括：
S41、使用C++中内置函数sort对所有分子构象进行排序，其中分子构象的能量越低，排名越靠前；
S42、选取排名靠前的若干个分子构象，使用BFGS对分子构象进行优化，使用打分函数计算优化后的分子构象的分数，根据分子构象在搜索空间中的位置,计算两个分子构象之间的均方根误差,若均方根误差的数值小于等于1，则仅保留两个分子构象中能量较低的构象。
作为本发明的进一步改进，S6包括：
S61、使用ΔvinaRF20对CPU端得到的所有分子构象重新打分，所述Δvina RF20的函数为：
pKd(ΔvinaRF)＝pKd(Vina)+ΔpKd(RF)；
其中，pKd(Vina)＝-0.73349E(Vina)；E(Vina)为AutoDock Vina打分函数所得的分数，ΔpKd(RF)为随机森林模型得到的修正项；
S62、按照ΔvinaRF20的打分结果对输出分子构象重新排序，并保留排名靠前的若干个分子构象。
为实现上述目的，本发明提供了一种基于GPU加速的分子对接装置，以执行前述的基于GPU加速的分子对接方法，包括：
数据准备模块，用于配置OpenCL环境、对接盒子的参数文件以及生成多个随机分子构象；
网格缓存生成模块，用于生成网格缓存表；
执行模块，用于优化所述随机分子构象；
精细搜索模块，对执行模块生成的分子构象进行优化；
打分模块，使用ΔvinaRF20对精细搜索模块输出的分子构象重新打分并排序。
作为本发明的进一步改进，所述执行模块每次接受一个配体的多个随机分子构象；
执行模块使用数据准备模块生成的随机分子构象进行蒙特卡罗的迭代局部搜索算法，其中需要计算构象的总能量，总能量包括分子间能量和分子内能量，分子间能量通过对网格缓存表进行三线性插值得出，分子内能量使用数据准备模块的分子内能量表计算得出；
BFGS单元，用于优化分子构象；
接受准则单元，用于判断是否接受BFGS单元优化后的分子构象。
本发明的有益效果是：本发明通过使用基于随机森林的打分函数ΔvinaRF20对分子构象进行重新打分并排序，使用OpenCL实现AutoDock Vina的GPU和CPU的异构并行结构，将运行时间较长的蒙特卡罗的迭代局部搜索算法以及后续复杂的受体与配体间的分子间能量计算时所需的网格缓存表交给GPU计算，不同配体的分子间能量计算均使用相同共享的网格查找表，剔除了CPU端的冗余计算，大大降低了时间和药物研发的成本，提高了输出分子构象的准确性，使得AutoDock Vina更适合用于真实情景下的药物虚拟筛选。
附图说明
图1是本发明的基于GPU加速的分子对接方法的步骤图。
图2是图1中基于GPU加速的分子对接方法的流程图。
图3是图1中S1的流程图。
图4是图1中S2的流程图。
图5是本发明中对接盒子的示意图。
图6是图1中S3的流程图。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明进行详细描述。
本发明提供了一种基于GPU加速的分子对接方法及装置，并且基于OpenCL实现AutoDock Vina的GPU和CPU的异构并行结构，减少了CPU和GPU之间的数据传输次数，降低了通信成本；通过引入ΔvinaRF20对AutoDock Vina输出的分子构象进行打分并重新排序，提高了输出的分子构象的准确性，使得AutoDock Vina更适合用于真实情景下的药物虚拟筛选。
需要说明的是，以下实施例中所用的GPU型号为NVIDIA GeForce GTX 1060 6GB，CPU型号为Intel(R)Core(TM)i7-8700 CPU@3.21GHz，当然，在其他实施例中，GPU的型号和CPU的型号可以根据实际情况进行选择，此处不作限制。
请参阅图1至图6所示，本发明的基于GPU加速的分子对接方法包括以下步骤：
S1：CPU端进行数据准备以及设置OpenCL环境，其中，数据准备包括对接盒子及对接盒子的参数准备，并将准备完成的数据传输到GPU端。
请参阅图3所示，S1具体包括：
S11：读取文件，包括：读取受体文件，受体文件的格式为.pdbqt。
具体的受体文件如表1和表2所示，由于受体的类型较多，故本实施例以2bm2为例，但不应以此为限：
表1
表2
读取对接盒子参数文件，对接盒子参数文件包括：
受体文件的路径；
配体库的路径，配体库包括所有配体；
对接盒子的大小，分别为size_x、size_y、size_z；
对接盒子的中心位点的三维坐标,分别为center_x、center_y、center_z。
对接盒子的参数文件如表3所示：
表3
其中，receptor：2bm2的受体文件所在路径；
ligand_directory：2bm2的所有配体文件所在路径；
output_directory：2bm2的输出文件所在路径；
center_x、center_y、center_z：对接盒子的中心位点的三维坐标；
size_x、size_y、size_z：对接盒子的大小；
thread：线程数的大小，默认值为8000。
本实施例中，所用线程数为8000，搜索步数为6，受体文件为2bm2_protein.pdbqt，配体文件夹为test，配体文件夹中包含30个配体文件，当然，在其他实施例中，具体的参数也可以根据实际情况进行设置。
S12：OpenCL环境设置，包括；
查询平台，使用OpenCL中clGetPaltformIDs()函数获取可用的OpenCL环境的数量，根据所得环境数量继续调用一次该函数初始化一个可用环境；
查询设备，使用OpenCL中clGetDeviceIDs()函数获取平台中的设备数量，根据所得设备数量继续调用一次该函数初始化一个可用设备；
创建上下文，上下文中包含内存、程序、内核对象等，使用clCreateContext()函数对初始化后的设备创建上下文，用于设备与程序以及设备与设备之间的通信；
创建命令队列，使用cl_queue_properties()的属性值创建命令队列，命令队列能够操作上下文中的多个对象；
创建程序对象，使用clCreateProgramWithSource()函数对OpenCL环境中的源代码创建程序对象，使用clBuildProgram()函数为指定的设备创建程序对象；
创建内核对象，使用clCreateKernel()函数创建两个内核对象；
设置内核参数，使用clSetKernelArg()函数对内核对象设置内核参数，便于执行内核对象；
执行内核，使用clEnqueueNDRangeKernel()函数调用命令队列，使得内核对象在设备上进行排队。
S13：生成分子内能量表，用于计算分子构象的分子内能量。
S14：根据随机数种子生成随机数查找表。
S15：将S11至S14的数据传至GPU端。
S2：对接盒子中包括若干个格点，GPU端的内核函数一用于计算格点和受体中的原子之间的相互作用力，并生成网格缓存表，并将网格缓存表存放至全局内存。具体的，假设格点处存在原子，则原子的类型共有17种。
请参阅图4和图5所示，S2具体包括：
S21：将对接盒子1划分成多个格子2，格点3位于格子2的顶角处，划分后格子2的大小为l*m*n，优选的，划分后格子2的大小为以得到格点3的量化坐标，内核函数一中包括与若干个格点3相对应的若干个线程，内核函数一中的线程数由格点3的数量决定，即线程与格点3的数量相对应，假设格点3处存在原子，每个线程计算其对应格点3处的原子作用力。
每个格点3在一个线程中进行以下操作：
S22：将格点处的量化坐标转换成其在受体4中的真实坐标，部分格点的量化坐标和真实坐标如表4所示：
表4
S23：以格点3为中心建立正方体5，优选的，正方体5的大小为该正方体5内包含受体4的原子，计算在该格点3处的原子与受体4的原子之间的相互作用力，将相互作用力进行存储，得到网格缓存表，并将网格缓存表存放在全局内存中，供所述GPU线程共享。
具体的，在正方体5内的受体4的原子中，如果该原子离格点3的距离小于则计算该格点3处与该原子之间的相互作用力，由于格点处的原子的类型为17种，格点处的原子可以得到17个相互作用力的数值；
部分格点3处得到17个相互作用力的数值如表5所示：
表5
其中，1-17是相互作用力的序号，(0.0.0)、(20.20.20)、(20.22.24)、(25.25.25)、(30.31.32)为部分格点的量化坐标，其他数值为相互作用力数值。
S24：所有线程完成以上操作后，按照原子类型，可以得到17张网格缓存表。将网格缓存表存放在全局内存中，供所有的GPU线程共享。
S3：CPU端生成配体的多个随机分子构象并传到GPU端，GPU端的内核函数二利用蒙特卡罗的迭代局部搜索算法优化随机分子构象并打分，得到多个构象，并将多个构象传输至CPU端。
具体的，CPU端将一个配体的多个随机分子构象传到GPU端的内核函数二中，GPU端的内核函数二使用CPU端使用蒙特卡罗的迭代局部搜索算法优化随机分子构象并打分，得到分数排名靠前的9个构象。
构象所得分数与构象的总能量有关，构象的总能量包含分子间能量与分子内能量，分子间能量由内核函数一所得的网格缓存表计算得出，分子内能量由分子内能量表计算得出，GPU端的内核函数二将构象传回CPU端。
构象表示公式为：
C＝{x,y,z,a,b,c,d,ψ1,ψ2,…,ψNrot}；
其中，C：构象，x、y、z：构象在搜索空间中的位置，a、b、c、d：构象的四元数，ψ1、ψ2、…、ψNrot：可旋转键的扭转角度，Nrot：配体可旋转键的个数。
配体文件夹中的一个配体文件如表6和表7所示：
表6
表7
如图6所示，GPU的内核函数二中每个线程生成最佳构象的流程图，S3具体包括：
S31：随机分子构象的表示方式如下：
随机分子构象在搜索空间中的三维坐标，Position＝(a,b,c)；
构象的四元数，orientation＝(R0,R1,R2,R3)；
构象中可旋转键的扭转角度，torsion＝(t0,t1,t2)；
S32、每个随机分子构象对应一个线程，每个线程中均进行蒙特卡罗的迭代局部搜索算法，本实施例中的线程数为8000。每个线程的搜索步数由启发式公式得出，公式为：max(1,floor(0.24*Natom+0.29*Nrot-3.41))，即取1与floor(0.24*Natom+0.29*Nrot-3.41)中的最大值作为搜索步数，其中Natom为配体的原子数目，Nrot为配体中可旋转键的数目，floor(0.24*Natom+0.29*Nrot-3.41)表示对0.24*Natom+0.29*Nrot-3.41得到的数值向下取整。
每个线程中均进行蒙特卡罗的迭代局部搜索算法，具体步骤为：
S321：对随机分子构象进行扰动，随机改变以下三组中的一组即可。
①对构象在搜索空间中的三维坐标进行基于均匀分布随机数产生的扰动；
②对构象的四元数进行基于正态分布随机数产生的扰动；
③对构象中可旋转键的扭转角度中某一具体维度值替换为一个符合均匀分布的随机数。
S322：对构象进行BFGS，具体包含以下内容：
a、初始化海森矩阵：将海森矩阵初始化为单位矩阵。
b、对随机分子构象进行三维空间的转换：即根据构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度以得到构象中每个原子的空间坐标。
c、使用打分函数计算构象的能量，构象的能量包含分子间能量与分子内能量，公式为e＝einter+eintra。遍历构象中的配体原子，使用三线性插值法查找全局内存中的网格缓存表近似计算配体与受体的分子间能量，即为einter；遍历各个有相互作用力的配体原子，使用分子内能量表精确计算配体的分子内能量，即为eintra。
d、计算构象中每个原子的空间坐标的导数，将其转换成构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数。
e、初始化步长alpha＝1，并使用Armijo-Goldstein准则进行黄金搜索算法得到新的构象，计算新的构象能量并计算构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数，使用更新前后的构象的能量判断步长是否需要更新，若需要更新则alpha＝alpha/2，继续进行下一轮黄金搜索算法，否则返回alpha值。
f、判断更新后构象的导数是否满足‖导数‖≥1e-5,若满足条件则更新海森矩阵。
g、使用Armijo-Goldstein准则判断e、f循环次数是否达到10次，若满足则输出为BFGS最终得到的构象。
S323：使用接受准则判断是否接受更新后的构象，具体包含以下内容：
第一次运行S322生成的构象一定被接受，如果此时S322生成的新的构象能量小于上一次生成的构象能量，此新构象一定被接受；如果此时S322生成的新构象能量大于上一次生成的构象能量，以一定概率接受该构象。
接受准则公式如下：
接受准则：
其中，P：接受优化后构象的概率，e0：扰动前构象的自由能，eopt：构象优化后的自由能。
S324：判断S322、S323的循环次数是否达到设定值，设定值为搜索步数。
S33、将内核函数二中每个线程得到的构象均传到CPU端。
S4：CPU端将S3中的多个构象进行优化和排序，得到分子构象。
S41、使用C++中内置函数sort对所有分子构象进行排序，其中分子构象的能量越低，排名越靠前；
S42、选取排名靠前的若干个分子构象，使用BFGS对分子构象进行优化，使用打分函数计算优化后的分子构象的分数，根据分子构象在搜索空间中的位置,计算两个分子构象之间的均方根误差,若均方根误差的数值小于等于1，则仅保留两个分子构象中能量较低的分子构象。
配体文件得到的分子构象数据如表8和表9所示。
表8
/>
表9
其中，REMARK VINA RESULT:-5.8 0.000 0.000中-5.8为该分子构象的分数。CPU端得到的分子构象数据存放在test_out文件夹中。
S5、将下一个配体的多个随机分子构象传到GPU端重复S3和S4，直至所述配体库中所有配体均经过S3和S4。
S6、使用ΔvinaRF20对CPU端中的所有分子构象进行打分并进行排序，保留排名靠前的若干个分子构象，实现药物的虚拟筛选。优选的，保留排名靠前的9个分子构象。
S61、使用ΔvinaRF20对CPU端得到的所有分子构象重新打分，所述Δvina RF20的函数为：
pKd(ΔvinaRF)＝pKd(Vina)+ΔpKd(RF)；
其中，pKd(Vina)＝-0.73349E(Vina)；E(Vina)为AutoDock Vina打分函数所得的分数，ΔpKd(RF)为随机森林模型得到的修正项。
S62、按照ΔvinaRF20的打分结果对输出的分子构象重新排序，并保留排名靠前的若干个分子构象。
本实施例运行后得到的9个构象的名称和分数如表10所示：
表10
2bm2_176的数据如表11、表12和表13所示：
表11
/>
表12
/>
表13
本申请还公开了一种新的基于GPU加速的分子对接装置，包括：
数据准备模块，用于配置OpenCL环境、对接盒子参数文件以及生成多个随机分子构象等；
网格缓存生成模块，用于生成网格缓存表；
执行模块，用于优化随机分子构象；
执行模块每次接受一个配体的多个随机分子构象；
执行模块使用数据准备模块生成的随机分子构象进行蒙特卡罗的迭代局部搜索算法，其中需要计算构象的总能量，总能量包括分子间能量和分子内能量，其中分子间能量通过对网格缓存表进行三线性插值得出，分子内能量使用数据准备模块的分子内能量表得出。
执行模块还包括：
BFGS单元，用于优化分子构象；
接受准则单元，用于判断是否接受BFGS单元得出的构象；
精细搜索模块，对执行模块生成的分子构象进行优化；
打分模块，对精细搜索模块输出的构象使用ΔvinaRF20重新打分排序。
综上所述，本发明提供了一种使用OpenCL实现AutoDock Vina的GPU和CPU的异构并行结构，使用两个GPU内核并行计算，将蒙特卡罗的迭代局部搜索算法以及后续复杂的受体与配体间的分子间能量计算时所需的网格缓存表交给GPU计算，缩短运行时间，剔除了CPU端的冗余计算，降低了通信成本；使用ΔvinaRF20对AutoDock Vina输出构象进行打分并重新排序，保留排名在前的多个构象，使得AutoDock Vina更适合用于真实情景下的药物虚拟筛选，即使用特定受体对大型配体数据库进行虚拟筛选并提高了输出分子构象的准确性。
以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。
Claims (9)
1.一种基于GPU加速的分子对接方法，用于将受体和配体库中的配体进行分子对接，其特征在于，包括如下步骤：
S1、CPU端进行数据准备以及设置OpenCL环境，其中，数据准备包括对接盒子及对接盒子的参数准备，并将准备完成的数据传输到GPU端；
S2、所述对接盒子中包括若干个格点，GPU端的内核函数一用于计算所述格点和所述受体中的原子之间的相互作用力，并生成网格缓存表；
S3、CPU端生成配体的多个随机分子构象并传到GPU端，GPU端的内核函数二利用蒙特卡罗的迭代局部搜索算法优化随机分子构象并打分，得到多个构象，并将多个构象传输至CPU端；
S4、CPU端将S3输出的多个构象进行优化和排序，得到分子构象；
S5、重复S3和S4，直至所述配体库中所有配体均完成S3和S4；
S6、使用ΔvinaRF20对CPU端中的所有分子构象进行打分并排序，保留排名靠前的若干个分子构象，实现药物的虚拟筛选；
其中，S1中数据准备包括：
S11、读取文件，包括：
读取受体文件，所述受体文件的格式为.pdbqt；
读取对接盒子参数文件，所述对接盒子参数文件包括：
所述受体文件的路径；
所述配体库的路径，所述配体库包括若干个所述配体；
所述对接盒子的大小，分别为size_x、size_y、size_z；
所述对接盒子的中心位点的三维坐标,分别为center_x、center_y、center_z；
线程数的大小为thread；
S12、所述OpenCL环境的设置，包括：
查询平台，使用所述OpenCL中的clGetPaltformIDs()函数获取可用的OpenCL环境的数量，并初始化一个可用的OpenCL环境；
查询设备，使用所述OpenCL中的clGetDeviceIDs()函数获取平台中设备的数量，并初始化一个可用的设备；
创建上下文，所述上下文中包括内存、程序以及内核对象，使用clCreateContext()函数对初始化后的设备创建上下文，用于设备与程序以及设备与设备之间的通信；
创建命令队列，使用cl_queue_properties()的属性值创建命令队列，所述命令队列能够操作所述上下文中的多个对象；
创建程序对象，使用clCreateProgramWithSource()函数对OpenCL环境中的源代码创建程序对象，使用clBuildProgram()函数为指定的设备创建程序对象；
创建内核对象，使用clCreateKernel()函数创建两个内核对象；
设置内核参数，使用clSetKernelArg()函数对所述内核对象设置内核参数，便于执行所述内核对象；
执行内核，使用clEnqueueNDRangeKernel()函数调用所述命令队列，使得所述内核对象在设备进行排队；
S13、生成分子内能量表，用于计算所述分子构象的分子内能量；
S14、根据随机数种子生成随机数查找表。
2.根据权利要求1所述的基于GPU加速的分子对接方法，其特征在于，S2具体包括：
将所述对接盒子划分成多个格子，所述格点位于所述格子的顶角处，划分后所述格子的大小为l*m*n，得到所述格点的量化坐标，所述内核函数一中包括与若干个所述格点相对应的若干个线程；
将所述量化坐标转换为受体中的真实坐标；
以所述格点为中心建立正方体，所述正方体内包含所述受体的原子，假设所述格点处存在原子，计算所述格点处原子与受体原子的相互作用力，将所述相互作用力进行存储，得到所述网格缓存表，并将所述网格缓存表存放在全局内存中，供所述GPU线程共享。
3.根据权利要求2所述的基于GPU加速的分子对接方法，其特征在于：所述正方体的大小为所述格点与所述配体原子之间的距离小于/>时，计算所述格点处原子与所述受体原子的相互作用力，所述格点处原子的类型共有17种。
4.根据权利要求1所述的基于GPU加速的分子对接方法，其特征在于，S3具体包括：
S31、随机分子构象的表示方式如下：
随机分子构象在搜索空间中的三维坐标，Position＝(a,b,c)；
分子构象的四元数，orientation＝(R0,R1,R2,R3)；
分子构象中可旋转键的扭转角度，torsion＝(t0,t1,t2)；
S32、每个随机分子构象对应一个线程，每个线程中均进行蒙特卡罗的迭代局部搜索算法，线程数为8000；
S33、将内核函数二中每个线程得到的分子构象均传到CPU端。
5.根据权利要求4所述的基于GPU加速的分子对接方法，其特征在于，S32包括：
S321、对随机分子构象进行扰动；
S322、对随机分子构象进行BFGS，包括：
a、将海森矩阵初始化为单位矩阵；
b、对随机分子构象进行三维空间的转换，得到随机分子构象中每个原子的空间坐标；
c、使用打分函数计算随机分子构象的能量，随机分子构象的能量包含分子间能量与分子内能量，通过所述网格缓存表计算分子间能量，通过分子内能量表计算分子内能量；
d、计算构象中每个原子的空间坐标的导数，将其转换成随机分子构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数；
e、使用Armijo-Goldstein准则进行黄金搜索得到新的分子构象，初始化步长alpha＝1，计算新的分子构象能量并计算分子构象在搜索空间中的三维坐标、四元数、所有可旋转键的扭转角度的导数，使用更新前后的分子构象的能量判断步长是否需要更新，若需要更新则alpha＝alpha/2，继续进行下一轮黄金搜索，否则结束搜索；
f、判断更新后分子构象的导数是否满足‖导数‖≥1e-5,若满足条件则更新海森矩阵；
g、判断e、f循环次数是否达到10次，若满足则输出为BFGS最终得到的分子构象；
S323、使用接受准则判断是否接收所述S322得到的分子构象；
S324、判断S322和S323的循环次数是否达到设定值，循环次数达到设定值后输出本线程的分子构象。
6.根据权利要求1所述的基于GPU加速的分子对接方法，其特征在于，S4包括：
S41、使用C++中内置函数sort对所有分子构象进行排序，其中分子构象的能量越低，排名越靠前；
S42、选取排名靠前的若干个分子构象，使用BFGS对分子构象进行优化，使用打分函数计算优化后的分子构象的分数，根据分子构象在搜索空间中的位置,计算两个分子构象之间的均方根误差,若均方根误差的数值小于等于1，则仅保留两个分子构象中能量较低的构象。
7.根据权利要求1所述的基于GPU加速的分子对接方法，其特征在于，S6包括：
S61、使用ΔvinaRF20对CPU端得到的所有分子构象重新打分，所述ΔvinaRF20的函数为：
pKd(ΔvinaRF)＝pKd(Vina)+ΔpKd(RF)；
其中，pKd(Vina)＝-0.73349E(Vina)；E(Vina)为AutoDockVina打分函数所得的分数，ΔpKd(RF)为随机森林模型得到的修正项；
S62、按照ΔvinaRF20的打分结果对输出分子构象重新排序，并保留排名靠前的若干个分子构象。
8.一种基于GPU加速的分子对接装置，以执行如权利要求1-7中任一项所述的基于GPU加速的分子对接方法，其特征在于，包括：
数据准备模块，用于配置OpenCL环境、对接盒子的参数文件以及生成多个随机分子构象；
网格缓存生成模块，用于生成网格缓存表；
执行模块，用于优化所述随机分子构象；
精细搜索模块，对执行模块生成的分子构象进行优化；
打分模块，使用ΔvinaRF20对精细搜索模块输出的分子构象重新打分并排序。
9.根据权利要求8所述的基于GPU加速的分子对接装置，其特征在于：所述执行模块每次接受一个配体的多个随机分子构象；
执行模块使用数据准备模块生成的随机分子构象进行蒙特卡罗的迭代局部搜索算法，其中需要计算构象的总能量，总能量包括分子间能量和分子内能量，分子间能量通过对网格缓存表进行三线性插值得出，分子内能量使用数据准备模块的分子内能量表计算得出；
BFGS单元，用于优化分子构象；
接受准则单元，用于判断是否接受BFGS单元优化后的分子构象。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210449861.2A CN114783545B (zh) | 2022-04-26 | 2022-04-26 | 基于gpu加速的分子对接方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210449861.2A CN114783545B (zh) | 2022-04-26 | 2022-04-26 | 基于gpu加速的分子对接方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114783545A CN114783545A (zh) | 2022-07-22 |
| CN114783545B true CN114783545B (zh) | 2024-03-15 |
Family
ID=82433107
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210449861.2A Active CN114783545B (zh) | 2022-04-26 | 2022-04-26 | 基于gpu加速的分子对接方法和装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114783545B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116646029B (zh) * | 2023-05-23 | 2024-02-13 | 正则量子(北京)技术有限公司 | 基于量子近似优化算法分子对接方法、装置、介质及设备 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107992718A (zh) * | 2017-11-28 | 2018-05-04 | 江苏理工学院 | 分子对接方法及系统 |
| CN114373509A (zh) * | 2021-12-21 | 2022-04-19 | 南京邮电大学 | 一种基于GPU加速AutoDock Vina的方法 |
| CN114860441A (zh) * | 2022-05-05 | 2022-08-05 | 南京邮电大学 | 一种面向药物虚拟筛选的多线程协作的gpu加速方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130141443A1 (en) * | 2011-12-01 | 2013-06-06 | Michael L. Schmit | Software libraries for heterogeneous parallel processing platforms |
-
2022
- 2022-04-26 CN CN202210449861.2A patent/CN114783545B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107992718A (zh) * | 2017-11-28 | 2018-05-04 | 江苏理工学院 | 分子对接方法及系统 |
| CN114373509A (zh) * | 2021-12-21 | 2022-04-19 | 南京邮电大学 | 一种基于GPU加速AutoDock Vina的方法 |
| CN114860441A (zh) * | 2022-05-05 | 2022-08-05 | 南京邮电大学 | 一种面向药物虚拟筛选的多线程协作的gpu加速方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114783545A (zh) | 2022-07-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhou et al. | Accelerating graph analytics on CPU-FPGA heterogeneous platform | |
| AU2019284011B2 (en) | Data processing method and related products | |
| Anderson et al. | General purpose molecular dynamics simulations fully implemented on graphics processing units | |
| CN103902702B (zh) | 一种数据存储系统和存储方法 | |
| Zhong et al. | Medusa: Simplified graph processing on GPUs | |
| CN104685499B (zh) | 执行过滤和投影操作的方法、系统和计算机可读介质 | |
| Bian et al. | Wide table layout optimization based on column ordering and duplication | |
| CN114783545B (zh) | 基于gpu加速的分子对接方法和装置 | |
| Baker et al. | An architecture for efficient hardware data mining using reconfigurable computing systems | |
| JP2014146319A (ja) | インメモリデータベースにおける高効率のゲノムリードアラインメント | |
| Chen et al. | Flow-guided file layout for out-of-core pathline computation | |
| CN110459263A (zh) | 一种基于bfgs算法的计算机药物筛选方法 | |
| Browne et al. | Forest packing: Fast parallel, decision forests | |
| JP2023007366A (ja) | 分子構造取得方法、装置、電子デバイス及び記憶媒体 | |
| CN102663207A (zh) | 一种利用gpu加速介观体系物理问题求解的方法 | |
| Liu et al. | Parallel reconstruction of neighbor-joining trees for large multiple sequence alignments using CUDA | |
| Gayatri et al. | Rapid exploration of optimization strategies on advanced architectures using testsnap and lammps | |
| Wu et al. | Multipredicate join algorithms for accelerating relational graph processing on GPUs | |
| Duan et al. | Hph: Hybrid parallelism on heterogeneous clusters for accelerating large-scale dnns training | |
| WO2022006771A1 (zh) | 分子力场多目标拟合算法库系统及工作流程方法 | |
| CN114373509A (zh) | 一种基于GPU加速AutoDock Vina的方法 | |
| Lupescu et al. | Design of hashtable for heterogeneous architectures | |
| Merelo-Guervós et al. | Scaling in distributed evolutionary algorithms with persistent population | |
| Wang et al. | Hardware-accelerated hypergraph processing with chain-driven scheduling | |
| CN105573834B (zh) | 一种基于异构平台的高维词汇树构建方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |