AutoDock Vina与含Zn蛋白对接流程经验分享

cpplh

最近研究涉及到小分子与含Zn蛋白的对接，但是常用的能处理金属蛋白对接的软件如薛定谔，MOE等通常需要购买。幸运的是，寻找后发现免费的AutoDock vina提供了相应流程（https://autodock-vina.readthedocs.io/en/latest/docking_zinc.html）。由于Zn²⁺配位具有明显的方向性和几何约束，官方教程推荐使用 AutoDock4Zn，也就是先用 AutoGrid4 生成 AD4Zn affinity maps，再用 Vina 的 --scoring ad4 进行对接。官方文档中说明了 AutoDock4Zn 是针对 zinc-coordinating ligands 扩展的力场，包含能量和几何两方面的 Zn 配位描述，并基于 292 个含锌晶体复合物优化。简要流程如下：

1. 受体准备→配体准备→Zn pseudo atom添加→GPF 生成→AutoGrid4 生成 maps→Vina 读取 maps 进行对接

复制代码

本文将从复现官方对接流程以及使用中碰到的一些问题这两方面进行分享，使用系统为VMware虚拟机-Ubuntu20.04.4
该流程分享基于官方流程改编而来，如有错误，还请批评指教
一、官方对接流程
1. Miniconda的安装

下载安装Miniconda

1. wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
   
2. bash Miniconda3-latest-Linux-x86_64.sh
   

复制代码

结束后，关闭然后重新打开终端，输入conda --version看有没有弹出conda版本，弹出类似conda XX.X.X就说明安装成功了。
2. 创建环境与安装vina
根据官网介绍，Python bindings 和 Vina 可执行程序是两套安装流程，安装 Python bindings 不一定包含 vina 可执行文件，反之亦然。为了降低出错概率，本文采用 conda-forge 统一安装。并且建议不要把所有东西装进 base 环境。Vina、Meeko、RDKit、Prody 等依赖较多，单独建环境可以减少版本冲突。安装命令如下

1. conda create -n vina_clean -c conda-forge python=3.10 numpy scipy rdkit vina meeko gemmi autogrid prody molscrub -y
   

复制代码

第一次输完conda create -n vina_clean -c conda-forge python=3.10 numpy scipy rdkit vina meeko gemmi autogrid prody molscrub -y后可能会弹出condaTosNonInteractiveError: Terms of Service have not been accepted for the following channels. Please accept or remove them before proceeding:，这是因为命令里有 -y，属于非交互安装，conda 不会自动弹出让你输入 yes，所以它要求你先手动接受它提示的两个条款。
按照提示输入如下命令即可，执行完后再重新输入上面的命令重新安装。

1. conda tos accept --override-channels --channel https://repo.anaconda.com/pkgs/main
   
2. conda tos accept --override-channels --channel https://repo.anaconda.com/pkgs/r
   

复制代码

其中-n vina_clean 代表创建一个叫vina_clean 的虚拟环境，可以根据自己需要命名。RDKit 是一个开源化学信息学工具包，可以用于读取 SDF/MOL2 等小分子结构，判断键级、价态、芳香性和形式电荷，检查是否存在 implicit H等。Meeko负责把小分子和受体结构转换成 Vina 可读的 PDBQT 文件。Gemmi 是结构生物学和晶体学相关的文件处理库，主要用于读取和处理 PDB、mmCIF、CIF、MTZ 等结构和晶体学文件。 autogrid用于生成受体周围的 affinity maps，随后vina通过读取这些.map 文件开始对接。Prody 是蛋白质结构、动力学和序列分析相关的 Python 包，Meeko通过其读取分析pdb文件。molscrub 提供 scrub.py 工具，主要用于 docking 前的小分子前处理，可以用于常规小分子加氢并基于UFF力场做一些简单优化。

3. 安装 ADFRsuite
直接去https://ccsb.scripps.edu/adfr/downloads/下载第一个linux版本安装包即可。随后按照官方下载界面右侧的installation mechanism输入命令即可安装成功，建议在输入安装命令前先输入如下命令

1. conda deactivate

复制代码

输入conda deactivate是因为建议不要把 ADFRsuite 加入conda 环境，这样可以避免和先前安装在conda环境里的 autogrid4、python 混淆。后续调用 ADFR 工具时，直接写完整路径即可，如：~/ADFRsuite-1.0/bin/pythonsh以及~/ADFRsuite-1.0/bin/autogrid4。其中pythonsh是ADFRsuite 自带的 Python 解释器，用来运行 ADFRsuite里的老脚本，autogrid4同之前的autogrid一样，用于生成 maps。

4. 官方案例下载
官方示例文件可以从https://github.com/ccsb-scripps/AutoDock-Vina点击右上角Code-Download ZIP下载后用unzip 命令解压。其中示例文件位于 AutoDock-Vina-develop/example/docking_with_zinc_metalloproteins/data ，Zn 专用脚本 zinc_pseudo.py ，prepare_gpf4zn.py 和AD4Zn.dat位于 AutoDock-Vina-develop /example/autodock_scripts。分别用于添加TZ pseudo atom以及生成Zn专用的GPF文件，而AD4Zn.dat包含了 Zn/TZ 相关参数。可以创建一个文件夹用于运行该示例，本文创建的文件夹路径为~/docking/zinc_docking_test。随后将上述文件全部拷入其中。

5. 准备受体
需要先提前处理好蛋白文件（如加H，去水，去配体等等，其中去水去配体可以使用pymol，加H可以看第二章），然后再执行下面命令：

1. conda activate vina_clean
   
2. cd ~/docking/zinc_docking_test
   
3. mk_prepare_receptor.py -i proteinH.pdb -o protein -p

复制代码

这一步可能会出现关于 FPP 的 warning，例如 FPP not in residue_templates。这通常不是报错，而是因为受体中含有辅因子 FPP，Meeko 会尝试根据结构信息为其建立模板。只要最终生成 protein.pdbqt，该步骤即可继续。
然后添加TZ原子以方便给配体定位

1. ~/ADFRsuite-1.0/bin/pythonsh zinc_pseudo.py -r protein.pdbqt -o protein_tz.pdbqt

复制代码

会提示Wrote 1 TZ atoms on protein_tz.pdbqt.说明该例子只添加了一个 Zn pseudo atom（TZ原子），这是因为Zn的其他配位点已被蛋白质中的残基所占据。

6. 准备配体

1. scrub.py 1s63_ligand.sdf -o 1s63_ligandH.sdf
   
2. mk_prepare_ligand.py -i 1s63_ligandH.sdf -o 1s63_ligand.pdbqt
   
3. 
   

复制代码

scrub.py用于给配体加氢并进行基于UFF力场的能量极小化。mk_prepare_ligand.py用于把配体转成 Vina 可读取的 PDBQT 格式

（别的配体加H方法见第二章）

7. 生成Zn专用GPF文件
GPF文件为grid parameter file，即格点参数文件，可以用来定义对接盒子以及参数

1. ~/ADFRsuite-1.0/bin/pythonsh prepare_gpf4zn.py -l 1s63_ligand.pdbqt -r protein_tz.pdbqt -o protein_tz.gpf -p npts=40,30,50 -p gridcenter=18,134,-1 -p parameter_file=AD4Zn.dat
   
2. 
   
3. -l：ligand，输入配体 PDBQT。
   
4. -r：receptor，输入带 TZ 的受体 PDBQT。
   
5. -o：输出的 GPF 文件命名。
   
6. -p npts=40,30,50：对接盒子在 x、y、z 方向上的格点数量。
   
7. -p gridcenter=18,134,-1：对接盒子中心坐标。
   
8. -p parameter_file=AD4Zn.dat：指定 Zn 专用参数文件。
   

复制代码

-p 用于设置 box center、size 以及 Zn 专用参数文件；AD4Zn.dat 包含 TZ atom type 定义，AutoGrid 需要在当前工作目录找到它，或者在 GPF 中就给出其路径。注意npts 单位不是 Å ，而是格点数。GPF 里默认格点间距是 0.375 Å，因此实际对接盒子大小为：

x 方向：40 × 0.375 = 15.0 Å
y 方向：30 × 0.375 = 11.25 Å
z 方向：50 × 0.375 = 18.75 Å
盒子的位置和尺寸可以通过AutoDockTools的grid box功能确定，b站上有非常多的教学视频可以看


8. 用 AutoGrid4 生成 affinity maps

1. ~/ADFRsuite-1.0/bin/autogrid4 -p protein_tz.gpf -l protein_tz.glg
   
2. 
   

复制代码

正常的话会生成很多 map 文件，例如：protein_tz.A.map，protein_tz.C.map，protein_tz.NA.map，protein_tz.OA.map，protein_tz.e.map，protein_tz.d.map，protein_tz.maps.fld这些 maps 相当于预先计算好的受体三维能量场。后续 Vina 不再直接根据受体坐标计算普通的 vina score，而是读取这些 AutoDock4/AutoDock4Zn maps。

9. 对接

1. vina --ligand 1s63_ligand.pdbqt --maps protein_tz --scoring ad4 --exhaustiveness 32 --out 1s63_ligand_ad4_out.pdbqt
   

复制代码

--ligand：输入配体 PDBQT。
--maps protein_tz：读取以 protein_tz 为前缀的 AutoDock4 maps。
--scoring ad4：使用 AutoDock4 scoring，而不是默认 vina scoring。
--exhaustiveness 32：搜索穷举度，越大搜索越充分，但耗时也越长，默认为8，官方教程上写的是32。
--out：输出 docking 结果命名。
当使用 AutoDock4 force field 时，只需提供 affinity maps 和 ligand，并用 --scoring ad4 指定使用 AutoDock4 force field即可。
还有一些有用的参数设置如--num_modes 20 --energy_range 5 --cpu 16。其中--num_modes 20代表最多输出 20 个对接构象 ，--energy_range代表只输出与最低能量构象能量差不超过5 kcal/mol的结构，其余不输出，--cpu 16表示计算的时候cpu用16核

10. 结果
官方教程说最佳 pose 的预测结合自由能应约为 -13 kcal/mol，我计算出来的如下，可以看到几乎一致

1. mode   affinity (kcal/mol)
   
2. 1       -12.94
   
3. 2       -12.37
   
4. 3       -12.25
   
5. ...
   

复制代码

二、替代方法与常见报错
1. 原子类型不支持
Autodock支持原子类型有限，仅支持C, H, O, N, S, P ,卤素以及部分金属（Zn/Mg/Ca/Mn/Fe）。
某些蛋白可能包含meeko不支持的原子，如含有K离子时可能会遇到如下报错：

1. Element K doesn't have an implemented covalent radius

复制代码

此时若不支持的元素对对接位点影响较小，可以先删除再做对接。

2. 存在多构象残基
蛋白中的一些链可能柔性较强，导致pdb文件中部分残基占有度不为1，存在多种构象，就可能遇到如下报错：

1. Residues with alternate location: ['A:709’]

复制代码

此时可以打开pdb文件看一下对于709号残基哪种构象占有度更高，如果占有度相同，可以比较下哪种更合理，或者干脆两种构象都做个对接看看哪个更合理，如果对口袋影响很小，那随便用哪个都无所谓。此时可以用如下命令解决该报错：

1. mk_prepare_receptor.py -i 蛋白.pdb -o 蛋白 -p --default_altloc A           #代表所有有多种构象的残基均用A构象
   
2. mk_prepare_receptor.py -i protein.pdb -o protein -p --wanted_altloc A:709=A               #代表只有残基709用A构象

复制代码

3. 蛋白质加H
对于含金属蛋白的加H和常规蛋白加H略有不同，这是因为蛋白中的一些残基会和金属离子配位，此时不能再往这些配位原子上加H。这里我推荐可以使用ChimeraX（https://www.cgl.ucsf.edu/chimerax/）来加H，导入蛋白质后，可以用如下命令加H：

1. addh hbond true metalDist 3.95
   

复制代码

其中hbond true表示加H的时候考虑氢键，可以判断 H 应该加在哪个位置、朝哪个方向，但会增加些许耗时；metalDist 3.95表示当具有强电负性的重原子 X（如氧原子或氮原子） 与金属离子 M 之间的距离小于3.95 Å 时将禁止为重原子 X 添加氢原子。详细解释可以看https://www.cgl.ucsf.edu/chimerax/docs/user/commands/addh.html。
4. 配体处理
官方是采用scrub.py脚本加H，但是这种方法通常只能处理一些常规有机小分子，最后还是采用的UFF力场优化结构，对于虚拟筛选很方便，但是对于需要对接的分子数量不多的情况下我认为没必要使用这个脚本。仅对于我个人而言，做完分子对接后，通常需要用动力学验证，因而就需要计算小分子的RESP电荷，那么前提自然就需要优化小分子结构，此时小分子的键长包括H的位置都是DFT方法计算出来的，比自动加H更加可靠。其次，很多分子如果是靠比如羟基跟金属离子配位的话，配位完之后氧上面的H通常会解离，此时拿结构为氧负的分子做分子对接应该更可靠更接近实际，但是scrub.py会给所有原子加H。对于这种情况我推荐使用Avogadro处理（https://avogadro.cc/docs/tools/draw-tool.html），对于数据库中下载的分子，将配体导入Avogadro后，首先检查分子的形式键级，如果不对，可以按ctrl+2或者右上角的黄色铅笔按钮进入分子编辑模式然后点击不对的那根键，每点一下键级就会加一，到三键后再点就重新变为单键，调整完形式键级后，点击构建-氢原子-添加氢原子来加氢，加完后再点击黄色铅笔然后看左边的“绘制“窗口，取消勾选调整氢原子数目，再在可视化界面里右键点击需要去掉的H原子。最后点击分析-属性-原子属性，选择去掉了H原子的杂原子（原子属性中选中的那一行原子在可视化窗口会被一个透明的篮色球体包住），更改其形式电荷，如氧原子去掉H后形式电荷需改为-1。对于DFT计算后的结构，只需导入进来后更改形式电荷即可。要注意的是对于某些结构比如吡啶环跟卤代烷取代后生成的吡啶阳离子结构，导入程序后往往没法正确判断形式键级与形式电荷，一定要检查然后修改。所有东西都改完后点击文件-导出-分子-选择SDF格式，因为其能记录形式键级与形式电荷，并且生成小分子pdbqt的mk_prepare_ligand.py脚本仅支持SDF，mol2等极少数格式。之所以形式键级与形式电荷如此重要，是因为mk_prepare_ligand.py处理小分子结构时会根据这两者判断分子中该有多少H原子，并且要求每个H原子都必须有明确的坐标，假设没改形式电荷，比如O去掉H后仍然是0，那么脚本就会认为这个O只有一个键，然后形式电荷为0，那么他就应该还有一个被隐藏起来的H没有给出具体坐标，即认为小分子的H没有加完整，就会有如下报错：

1. RDKit molecule has implicit Hs. Need explicit Hs

复制代码

即说明结构中有隐藏起来的H，此时极大概率是由于形式电荷或者形式键级没设置对导致的，这一点也可以通过SDF末尾检查，如M  CHG  1   5  -1就代表体系里有1个带电原子，为5号原子，带1个负电。

5.  可旋转键的检查
上述提到的mk_prepare_ligand.py工作流程可以如下大致表示

1. mk_prepare_ligand.py
   
2. 输入 SDF / MOL2
   
3.         ↓
   
4. RDKit 读取分子
   
5.         ↓
   
6. 检查价态、键级、形式电荷、氢原子、坐标
   
7.         ↓
   
8. Meeko 创建 MoleculeSetup
   
9.         ↓
   
10. 分配 AutoDock atom types
    
11.         ↓
    
12. 分配原子电荷（Gasteiger电荷）
    
13.         ↓
    
14. 识别可旋转键
    
15.         ↓
    
16. 建立 torsion tree
    
17.         ↓
    
18. 将非极性 H 合并（AutoDock采用联合原子模型，会将非极性H并入所连重原子上，极性H保留）
    
19.         ↓
    
20. 输出 PDBQT

复制代码

可以看到当中有识别可旋转键并建立 torsion tree，用过AutoDockTools的应该知道在将小分子选为配体后可以在Ligand-torsion tree-choose torsions中查看可以旋转的键，但有时候程序判断的会不准，对于mk_prepare_ligand.py也是如此。因而有必要生成完之后做一个检查，打开pdbqt文件，可以看到类似于BRANCH 1 2的信息，这就表示1号原子与2号原子之间的键可以旋转，在文件末尾可以看到类似于TORSDOF 10的信息，代表总共有10个可旋转键。也可以用如下命令快速查看：

1. grep -E "^(BRANCH|TORSDOF)" 配体.pdbqt

复制代码

需要注意的是，流程中有一步会将非极性 H 合并，因而不能拿原来的文件的原子序号与pdbqt中的原子序号对应，一个很简单的办法就是将pdbqt文件导入Avogadro，然后左侧勾选标签，在标签设置中将原子标签那一栏换成原子序号来与pdbqt文件对照检查可旋转键。如果需要固定某些键，可以使用如下命令：

1. mk_prepare_ligand.py -i 配体.sdf -o 配体_fixed.pdbqt --rigidify_bonds_smarts "[CX3](=O)[NX3]" --rigidify_bonds_indices 1 3

复制代码

--rigidify_bonds_smarts "[CX3](=O)[NX3]" --rigidify_bonds_indices 1 3 表示固定“[CX3](=O)[NX3]”中的1，3号原子，如固定酰胺键中的C-N键。
命令具体细节见https://meeko.readthedocs.io/en/develop/lig_cli_options.html，smarts写法可以看https://www.daylight.com/dayhtml/doc/theory/theory.smarts.html和https://rdkit.readthedocs.io/en/latest/RDKit_Book.html。我还搜索出来一个SMARTSeditor 图形工具，据说可以交互式创建和编辑 SMARTS，但还没有尝试，网址为https://www.zbh.uni-hamburg.de/en/forschung/amd/software/smartseditor.html。