摘要:本文演示了如何用Pymol与Openbabel从蛋白复合物结构出发准备AutoDock Vina计算所需的蛋白结构文件、配体文件、对接参数文件,并进行VINA分子对接计算。本文获得的这些文件,可进一步用于Cloudam VirtualFlow的高通量虚拟使用。

前言

抗白血病药物格列卫是ABL激酶抑制剂,1IEP晶体结构为格列卫(STI, GLEEVEC)与ABL激酶的复合物结构。在本练习中,演示了用Pymol从1IEP复合物结构出发准备起始的受体结构文件,并用在云E算力平台上预装的OpenBabel与分子对接软件AutoDock VINA完成:1)蛋白结构准备;2)配体结构准备;3)分子对接参数设置;4)分子对接;5)结果文件的格式转化;6)可视化分析。本文准备的文件可进一步用于Virtualflow进行高通量虚拟筛选。

本教程的目的

本教程的主要目的是初步了解:

  1. 如何用PYMOL准备蛋白结构
  2. 如何用OpenBabel转化PDB格式为PDBQT
  3. 如何准备Vina的对接参数文件
  4. 如何提交一个简答的VINA对接计算

用到的软件

  1. PYMOL
  2. 用于蛋白的准备,需要预先安装于本地机器上。

    Pymol的获取:https://github.com/schrodinger/pymol-open-source

  3. OpenBabel
  4. 用于格式转化,源代码与使用手册:http://openbabel.org

    OpenBabel已经预装于Cloudam云E算力平台,用module命令加载:

    1
    2
    3
    
    module add Anaconda3/2020.02
    source activate
    conda activate rdkit
  5. Vina家族软件
  6. 分子对接计算引擎VINA已经预装于Cloudam云E算力平台的virtualflow,可以复制到本地使用:

    1
    2
    
    export VFVS_ROOT=/public/software/.local/easybuild/software/virtualflow
    cp -r  $VFVS_ROOT/vfvs/tools .

    更简单的办法是通过PATH直接使用:

    1
    2
    
    export VFVS_ROOT=/public/software/.local/easybuild/software/virtualflow
    export PATH=$VFVS_ROOT/vfvs/tools/bin:$PATH

    Vina家族对接软件可执行文件都在VFTools/bin目录里。更多的AutoDock软件包介绍见:https://ccsb.scripps.edu/autodock以及Stefano等人[1]的使用教程。

本教程需要了解的预备技能

  • 对Linux的shell命令有初步了解
  • 对分子可视化软件Pymol有初步了解

操作步骤

Cloudam教程 | AutoDock Vina分子对接软件的使用-墨灵格的博客

Figure 1. 用Pymol与OpenBabel进行结构准备的操作流程

1. 下载蛋白

在Pymol命令行键入:

1
fetch 1iep

2. 将蛋白A链单独保存

由于1IEP有两个链,我们仅需要其中一个(Chain A),在Pymol命令行键入:

1
create 1iepA, Chain A

3. 删除一个对象

为了操作简单,删除掉一个对象:

1
delete 1iep

4. 删除水分子

1
remove resn HOH

5. 加极性氢

1
h_add elem O or elem N

6. 将复合物的配体复制为一份单独的对象

1
create ligand,RESN STI

7. 将复合物的蛋白复制为一份单独的对象

1
create protein, 1iepA and not RESN STI

8. 保存蛋白为本地PDB文件

1
save d:/work/docking/protein.pdb protein

9. 保存配体为本地PDB文件

1
save d:/work/docking/ligand.pdb ligand

10.将配体(STI)周围4.2A残基定义为结合位点并创建对象

1
create pocket, (byres (RESN STI round 4.2))

该操作的目的是可以对结合位点进行更好的渲染,不是必须的。

11. 将受体蛋白PDB格式转为PDBQT格式

本步骤将将第8步保存的受体文件protein.pdb转为PDBQT文件,用于后续分子对接计算用。

1
obabel -ipdb protein.pdb -opdbqt -xr -O protein.pdbqt

12. 将配体分子格式转化为PDBQT格式

以ZINC下载的mol2文件为例(在实际研究中,可能是其它软件比如chem3d保存的一个sdf或mol2文件),用openbabel转为配体用的PDBQT文件,用于后续分子对接计算用。

1
obabel -imol2 ZINC000019632618.mol2 –opdbqt -O ZINC000019632618.pdbqt
Cloudam教程 | AutoDock Vina分子对接软件的使用-墨灵格的博客

Figure 2. 对接参数文件的准备

注意OpenBabel对PDBQT的分支或可旋转键的定义可能会发生错误,以式1化合物为例来说明。

Cloudam教程 | AutoDock Vina分子对接软件的使用-墨灵格的博客

式1化合物,具有三个柔性键

lig_90.sdf是它的一个初始3D结构,用OpenBabel 3进行格式转化:

1
obabel -isdf lig_90.sdf –opdbqt -O lig_90_ob.pdbqt

会发现,Openbabel生成的pdbqt文件(lig_90_ob.pdbqt)仅包含两个可旋转键:C1_C12与C12_N18,具体见PDBQT文件的第4行与第5行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
REMARK  Name = lig_90
REMARK  2 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    1  A    between atoms: C_1  and  C_12
REMARK    2  A    between atoms: C_12  and  N_18
REMARK                            x       y       z     vdW  Elec       q    Type
REMARK                         _______ _______ _______ _____ _____    ______ ____
ROOT
ATOM      1  C   UNL     1      -6.646  -0.766   0.406  0.00  0.00    -0.027 A 
ATOM      2  C   UNL     1      -7.226   0.465   0.208  0.00  0.00    +0.009 A 
ATOM      3  C   UNL     1      -8.581   0.631   0.388  0.00  0.00    +0.050 A 
ATOM      4  C   UNL     1      -9.389  -0.424   0.758  0.00  0.00    +0.162 A 
ATOM      5  O   UNL     1     -10.699  -0.156   0.932  0.00  0.00    -0.488 OA
ATOM      6  C   UNL     1     -11.436  -0.874   1.946  0.00  0.00    +0.232 C 
ATOM      7  C   UNL     1     -11.914  -2.258   1.457  0.00  0.00    +0.082 C 
ATOM      8  C   UNL     1     -10.816  -2.990   0.656  0.00  0.00    +0.232 C 
ATOM      9  O   UNL     1      -9.502  -2.805   1.229  0.00  0.00    -0.488 OA
ATOM     10  C   UNL     1      -8.810  -1.692   0.920  0.00  0.00    +0.162 A 
ATOM     11  C   UNL     1      -7.447  -1.830   0.758  0.00  0.00    +0.054 A 
ENDROOT
BRANCH   1  12
ATOM     12  C   UNL     1      -5.154  -1.022   0.240  0.00  0.00    +0.139 C 
BRANCH  12  13
ATOM     13  N   UNL     1      -4.292   0.168   0.332  0.00  0.00    -0.309 N 
ATOM     14  H   UNL     1      -4.590   0.893   0.726  0.00  0.00    +0.149 HD
ATOM     15  C   UNL     1      -2.926   0.047   0.165  0.00  0.00    +0.215 A 
ATOM     16  N   UNL     1      -2.222  -0.392   1.230  0.00  0.00    -0.222 NA
ATOM     17  C   UNL     1      -0.899  -0.464   1.024  0.00  0.00    +0.131 A 
ATOM     18  C   UNL     1      -0.325  -0.124  -0.172  0.00  0.00    +0.054 A 
ATOM     19  C   UNL     1      -1.136   0.298  -1.191  0.00  0.00    +0.131 A 
ATOM     20 BR   UNL     1       1.534  -0.234  -0.411  0.00  0.00    -0.047 Br
ATOM     21  N   UNL     1      -2.467   0.400  -1.053  0.00  0.00    -0.222 NA
ENDBRANCH  12  13
ENDBRANCH   1  12
TORSDOF 2

胺基嘧啶片段没对应的分支(BRANCH)没有出现在文件里,键N18_C17也没有出现在可旋转键列表里,这与ADT的prepare_ligand4.py准备的结果不同,可能会导致在分子对接过程中N18_C17键是刚性的而没有机会获得正确的结合模式。我建议尽量用ADT的prepare_ligand4.py来准备结构,或者用其它更准确的软件。

此外,电荷也是错误的:没有将删掉的氢上电荷合并到键连的重原子上,这将导致在采用力场打分错误(不影响VINA的docking,因为vina不使用用户给的电荷)。更多的PDBQT格式说明,请参见:What is the format of a PDBQT file?

进一步分析发现所有包含类似胍基片段的N-C键都被OpenBabel识别为不可旋转键,比如式2所示的常见2-氨基嘧啶片段。

Cloudam教程 | AutoDock Vina分子对接软件的使用-墨灵格的博客

式2片段结构,具有三个柔性键

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
(base) [gkxiao@master demo]$ cat S2.smi
C1(NC2=NC(C3=CN=CC=C3)=CC=N2)=CC=CC=C1 S2
(base) [gkxiao@master demo]$ obabel -ismi S2.smi -opdbqt --gen3d --best
REMARK  Name = S2
REMARK  2 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    1  A    between atoms: C_1  and  N_2
REMARK    2  A    between atoms: C_5  and  C_6
REMARK                            x       y       z     vdW  Elec       q    Type
REMARK                         _______ _______ _______ _____ _____    ______ ____
ROOT
ATOM      1  N   UNL     1       2.812  -0.545   0.128  0.00  0.00    +0.000 N
ATOM      2  C   UNL     1       3.296  -1.851   0.224  0.00  0.00    +0.000 A
ATOM      3  N   UNL     1       3.942  -2.233   1.334  0.00  0.00    +0.000 NA
ATOM      4  C   UNL     1       4.393  -3.501   1.427  0.00  0.00    +0.000 A
ATOM      5  C   UNL     1       4.180  -4.407   0.410  0.00  0.00    +0.000 A
ATOM      6  C   UNL     1       3.504  -3.932  -0.691  0.00  0.00    +0.000 A
ATOM      7  N   UNL     1       3.060  -2.669  -0.809  0.00  0.00    +0.000 NA
ATOM      8  H   UNL     1       2.532  -0.383  -0.830  0.00  0.00    +0.000 HD
ENDROOT
BRANCH   1   9
ATOM      9  C   UNL     1       3.354   0.587   0.763  0.00  0.00    +0.000 A
ATOM     10  C   UNL     1       4.730   0.740   0.977  0.00  0.00    +0.000 A
ATOM     11  C   UNL     1       5.237   1.885   1.596  0.00  0.00    +0.000 A
ATOM     12  C   UNL     1       4.371   2.898   1.998  0.00  0.00    +0.000 A
ATOM     13  C   UNL     1       3.003   2.773   1.765  0.00  0.00    +0.000 A
ATOM     14  C   UNL     1       2.501   1.632   1.133  0.00  0.00    +0.000 A
ENDBRANCH   1   9
BRANCH   4  15
ATOM     15  C   UNL     1       5.110  -3.815   2.675  0.00  0.00    +0.000 A
ATOM     16  C   UNL     1       6.100  -4.795   2.757  0.00  0.00    +0.000 A
ATOM     17  N   UNL     1       6.791  -5.119   3.879  0.00  0.00    +0.000 NA
ATOM     18  C   UNL     1       6.496  -4.421   4.994  0.00  0.00    +0.000 A
ATOM     19  C   UNL     1       5.546  -3.417   5.038  0.00  0.00    +0.000 A
ATOM     20  C   UNL     1       4.853  -3.116   3.867  0.00  0.00    +0.000 A
ENDBRANCH   4  15
TORSDOF 2

可以看到,键C1-N2与C5-C6是可旋转的,而键N2-C3被视为不可旋转:这要求初始的N2-C3一定是对的才有可能在后续的docking中得到正确的pose。OpenBabel的这个行为很危险,尤其是激酶抑制剂分子对接虚拟筛选时,比如PDB 1IEP的redocking计算(测试文件:1iep.zip),使用OpenBabel准备的配体文件(sti_ob.pdbqt)进行的Vina docking不能获得正确的pose,而使用adt准备的配体文件(sti_adt.pdbqt)却可以重现结合模式。

接下来开始准备vina的对接参数文件。

13. 通过共晶配体定义结合位点的中心

在第10步,我们将共晶配体提取出来,并保存为PDB格式。可以利用该配体来定义结合位点的中心坐标:

1
obabel –ipdb ligand.pdb –obox

其中输出如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
HEADER    CORNERS OF BOX
REMARK    CENTER (X Y Z)          15.190     53.903     16.917
REMARK    DIMENSIONS (X Y Z)      10.664     18.739     15.526
ATOM      1  DUA BOX     1       9.858  44.533   9.154
ATOM      2  DUA BOX     1      20.522  44.533   9.154
ATOM      3  DUA BOX     1      20.522  44.533  24.680
ATOM      4  DUA BOX     1       9.858  44.533  24.680
ATOM      5  DUA BOX     1       9.858  63.272   9.154
ATOM      6  DUA BOX     1      20.522  63.272   9.154
ATOM      7  DUA BOX     1      20.522  63.272  24.680
ATOM      8  DUA BOX     1       9.858  63.272  24.680
CONECT    1    2    4    5
CONECT    2    1    3    6
CONECT    3    2    4    7
CONECT    4    1    3    8
CONECT    5    1    6    8
CONECT    6    2    5    7
CONECT    7    3    6    8
CONECT    8    4    5    7

我们关注第2行,给出了配体的质心,我们利用该质心用于定义Vina的x,y,z中心坐标。

14. Vina的参数文件

Vina计算所需的参数可以通过键入vina命令来获得。参数可以从命令行输入,也可以从参数文件读入。假设我们的参数文件为config.txt,部分参数设置如下:

1
2
3
4
5
6
7
8
9
receptor = protein.pdbqt   #第11步准备的蛋白文件
center_x = 15.190          #第13步获得
center_y = 53.903          #第13步获得
center_z = 16.917          #第13步获得
size_x = 30.0
size_y = 30.0
size_z = 30.0
exhaustiveness = 8
cpu = 1

第7行exhaustiveness的值越大,AutoDock Vina搜索的越充分,计算所需的时间与之成线性增加;如果第8行不指定CPU的数量,则使用全部的cpu核心进行计算,所用CPU数越多,完成计算所需的挂钟时间越短。

15. 用VINA进行分子对接

现在我们用vina将第12步准备的配体对接到蛋白结合位点里,用纯命令行方式输入参数:

1
2
3
4
5
6
7
8
9
10
11
12
vina --receptor protein.pdbqt \
         --ligand ZINC000019632618.pdbqt \
         --center_x 15.1 \
         --center_y 53.9 \
         --center_z 16.9 \
         --size_x 30 \
         --size_y 30 \
         --size_z 30 \
         --out vina.pdbqt \
         --log vina.log \
         --exhaustiveness 20 \
         --num_modes 20

如果使用第14步准备好的参数文件,则命令可得到简化:

1
2
3
4
vina --config config.txt \
     --ligand ZINC000019632618.pdbqt \
     --out vina.pdbqt \
     --log vina.log

其中--out为输出文件(pose预测结果);其中--log为日志文件。预测的pose可以直接重命名为pdb格式,用视图软件分析,比如pymol或Flare Viewer。其中Flare Viewer比起Pymol更加简单易用。

Figure 3. Flare分析STI与BCR-ABL的相互作用

避免误解的说明

本文采用PyMol进行初始的蛋白结构准备,仅是出于完成AutoDock Vina计算任务最基本的要求,并没有全面考虑蛋白结构准备的各个方面,比如:原子或残基缺失的修复、质子化状态与互变异构体的处理、侧链构象的选择、水分子的保留与否、金属离子的处理等等。正式的研究中应该充分考虑这些问题。PyMol并非蛋白结构准备的必须,也并非优选工具。有很多优秀的软件可以更可靠的进行结构准备。

OpenBabel用于将蛋白的PDB格式文件与配体的mol2格式文件转化为PDBQT格式文件非常方便,但是官方标配是AutoDockTools。并不了解OpenBabel与AutoDockTools两者之间的差异。读者需自行了解,但是OpenBabel无疑是非常方便的。

接下来可以做什么

  • Virutalflow进行超高通量虚拟筛选
  • 以预测的pose为起点用FEP计结合亲和力
  • 以预测的pose为起点进一步用分子动力学模拟考察结合的稳定性
  • 以预测的pose为起点用基于结构的方法生成药效团模型

文献

  • Forli, S.; Huey, R.; Pique, M. E.; Sanner, M. F.; Goodsell, D. S.; Olson, A. J. Computational Protein–Ligand Docking and Virtual Drug Screening with the AutoDock Suite. Nat. Protoc. 2016, 11 (5), 905–919. https://doi.org/10.1038/nprot.2016.051.