Docking算法需要每个原子带有电荷并且需要标记原子的属性。这些信息通常未包含在PDB文件中。我们需要在对蛋白和小分子的PDB文件预处理,生成PDBQT文件同时包含以上信息和PDB文件中的原子坐标信息。进一步地对于“柔性配体docking”,我们还需要定义配体的柔性部分和刚性部分。所有这些都可以通过软件AutoDock Tools (adt)来完成。
=====准备受体蛋白=======
PDB文件(1hsg.pdb)中包含了蛋白、配体和水分子(从显示可以看出配体在受体的中心区域);首先提取出蛋白的坐标,即以关键字ATOM和TER开头的行 ,存储到文件1hsg_prot.pdb。
注:在windows下,我们可以手动选择,或者利用Excel的筛选功能。在Linux下,使用命令egrep "^(ATOM|TER)" 1hsg.pdb >1hsg_prot.pdb
2. 启动AutoDockTools
3. 依次点选File-ReadMolecule-1hsg_prot.pdb加载蛋白分子。
按住左键拖动旋转分子结构;点击中键滚动缩放;按住右键移动晶体位置。
4. 更改展示方式:依次点选Color-By Atom Type-All Geometries-OK。
5. 加氢:晶体结构中通常缺少氢原子的坐标 (因为氢原子电子少,且质子核对电子吸引能力弱,因此很难定位,具体见http://www.uh.edu/~chembi/ChemSocRev_Jones_critical.pdf)。但是在docking过程中,氢原子尤其是极性氢原子对计算静电作用是必须的。因此我们需要给蛋白加上氢原子,依次点选Edit-Hydrogen-Add-Polar only-OK。这时氢原子会以白色短线形式出现。
6. 存储对蛋白的每个原子所做的修改和原子类型判断:依次点选Grid-Macromolecule-Choose-1HSG_protein-Slect Molecule。ADT会弹出一个信息框包含程序所做的处理,比如合并非极性氢原子,计算原子局部电荷和判断原子类型,并提示保存Save-1hsg_prot.pdbqt。打开文件,查看最后两列,分别为每个原子的电量和类型 。如下图所示分别是原始的pdb文件和新生成的pdbqt文件。
7. 在受体蛋白定义配体结合的3D搜索空间: 如果我们事先不知道结合位点,理论上可以定义一个长方体盒子包含整个蛋白或者随便一个特定区域 。
依次点选Grid-Grid box将会在蛋白上画出一个长方体,并且有一个弹出框。在弹出框中,拖拽刻度线查看长方体的变化,完成设置。
在这个例子中,我们知道结合位点,就选取以其为中心的一个小空间。设置Spacing (angstrom)为1埃 (这实际是一个换算系数, 相当于步长; 默认为0.375,是C-C单键长度的1/4,最大为1。spacing值与(各个维度上的点的数目+1)的乘机就是长方体Grid box的大小)。在我们调整的过程中,可以看到随着这个数值的变大,立方体也被放大了。另外我们设置x,y,z center为16,25,4,number of points in (x,y,z)-dimension为30,30,30(最大为126,必须为偶数,AutoDock会自动再每一维再加一个点)。记下我们设置的这些点,下面会用到。
在刻度转盘处点击右键会弹出一个窗口,输入数字回车即可设置GRID的中心坐标和大小。较大的number of points in (xyz)-dimension和较小的Spacing会增加搜索的精度,同时需要花费更多的计算时间。
8. 设置受体的柔性残基:在ADT中依次点选Flexible Residues-Input-Choose Macromolecule-1hsg_prot; select-select from string-Residue: ARG8-Add-Dismiss, 8号ARG氨基酸残基就被选中了。
再依次点选Flexible Residues-Choose Torsions in Currently Selected Residues将选择的残基标记为柔性残基并设置可扭转的数量。在分子显示窗口中分别点击两个残基上CA和CB原子之间的建,使之变为非扭转的(紫色显示),这样两个残基中的32个键中有6个是可扭转的。这里设置配体的柔性残基或者使CA-CB的键为刚性都是可选操作。
Flexible Residues-Output-Save Flexible PDBQT保存柔性残基文件。Flexible Residues-Output-Save Rigid PDBQT保存柔性残基文件。
9. 关掉grid和删除protein:Grid Options-File-Close w/out saving; Edit-Delete-Delete Molecule-1hsg_prot-Continue。
=====准备配体=====
蛋白结构类似,配体的结构也缺少氢原子,我们需要添加氢原子并且定义哪些键是可以旋转的以用于柔性docking。从PDB结构中提取配体的原子位置。indinavir的配体残基名字为MK1,以HETATM开头的行表示非核心多聚体的成分 (heteroatoms)。
Linux系统下,运行grep "^HETATM.*MK1" 1hsg.pdb >indinavir.pdb
Windows系统下,直接拷贝到文件indinavir.pdb
3. 将结构读入ADT;依次点选File-ReadMolecule-indinavir.pdb;Color-By Atom Type-All Geometreies-OK。
4. 晶体结构中通常缺少氢原子 (因为氢原子电子少,且质子核对电子吸引能力弱,因此很难定位)。但是在docking过程中,氢原子,尤其是极性氢原子对计算静电作用是必须的。因此我们需要给配体加上氢原子,Edit-Hydrogen-Add-Polar only-OK(之所以选择Polar only是因为vina的官方视频里面是这么选择的)。这时氢原子会以白色短线形式出现。
5. 在ADT中定义此化合物为配体,以便ADT为其计算局部电荷(partialcharges)和设置可旋转配体键。依次点选Ligand-input-Choose-indinavir-Select Molecule for AutoDock4。这时会有一个弹出框显示ADT所做的操作,包括合并非极性氢(只在添加了的情况下)、计算电荷电量和设置旋转键。然后点选Ligand-Output-Save as PDBQT存储结果。
6. 查看ADT检测出的旋转键,依次点选Ligand-TorsionTree-Choose Torsions,可以看到Number of rotatable bonds=14/32。
====准备docking的配置文件====
docking配置文件包含了输入的受体(蛋白)、配体(化合物)和搜索参数的信息,为一个文本文件,名字任意,可以为conf.txt,内容如下
receptor = 1hsg_prot.pdbqt
ligand = indinavir.pdbqt
num_modes = 50
out = dockingResult.pdbqt
log = docking.log
center_x = 16
center_y = 25
center_z = 4
size_x = 30
size_y = 30
size_z = 30
seed = 2009
receptor和ligand为输入文件的名字,与conf.txt在同一目录下; out为输出文件的名字;log为输出日志文件的名字。centerhe和size定义搜索空间的位置和大小。num_modes设置最多显示的结合模型,鉴于只输出符合能量值要求的结果,最后输出的结合模型数量可能少于这一数值。seed设置随机数生成的种子,可以为任意整数。如果想重现之前的分析结果就需要使用相同的seed。
====Docking 小分子化合物indinavir到HIV-1蛋白酶======
使用AutoDockVina执行docking预测。
在windows命令行提示符或linux终端下运行命令
vina--config conf.txt,大概需要几分钟时间。
2. 输出结果包含两个文件,构象文件dockingResult.pdbqt和日志文件docking.log。
dockingResult.pdbqt: 包含所有docking的模式,通常第一个为结合最好的构象,但如果前几个能量值相差不大时也有例外。
docking.log: 日志文件,包含结合能量值(第一列,越低越稳定,默认由低到高排序,所以第一个为最好的构象)、每个构象与第一个构象的距离、每个构象与第一个构象的差别。
====用PyMOL可视化docking结果====
打开PyMOL,依次点选File-Open文件类型选择All Files-选取结果dockingResult.pdbqt文件、原始蛋白和配体的pdb文件、原教程的pdbqt文件。
这个是原始的蛋白受体和配体信息:
受体蛋白信息:
配体信息:
预测的配体信息:
本文使用 文章同步助手 同步