AutoDock使用札记(1)分子对接
从最初学习AutoDock到现在已经有4年的时间。据报道，AutoDock是引用最多的分子对接软件，今年又有了最新的版本。使用多了，就想总结一下。用了一周的时间，终于总结完了。觉得和朋友分享一下就更好了。不完善之处还希望大家能给出意见！AutoDock在2.0以前使用的是模拟退火算法（Simulated Annealing Algorithm），在其后的版本中改为了拉马克遗传算法(Lamarckian Genetic Algorithm,LGA)。Autodock和Dock一样，都是基于格点（grid）的计算，首先用围绕受体活性位点的氨基酸残基形成box，然后用不同类型的原子作为探针（probe）进行扫描，计算格点能量，然后对配体在box范围内进行构象搜索（conformational search），最后根据配体的不同构象(conformation),方向（orientation）和位置（position）进行评分（scoring），排序（ranking）。
AutoDock虽然提供了诸如ADT(AutoDock Tools)，BDT的图形界面工具，但是我们仍然认为其是一个字符界面运行的分子对接软。由于运算速度不快，所以在字符界面下对于运算的管理更为方便，可以将长时间的运算放在后台运行。对于输入文件的准备问题，Autodock有其自带的一些程序，比如addsol等，同时ADT(AutoDock Tools)也提供了一些python脚本（路径：MGLTools/MGLToolsPckgs/AutoDockTools/Utilities24），可以代替这些程序来进行输入文件的处理。
本手册按照以下顺序编写：
Autodock分子对接流程
AutoDock结果分析
AutoDock用于虚拟筛选
AutoDock Tools简介
限于篇幅，本贴先介绍分子对接及结果分析，虚拟筛选等另外开贴介绍。
※※※※※※ 手册以外的话 ※※※※※※
1、autodock只能对单分子文件对接(一个文件只含有一个分子），不能对多分子文件筛选，所以，对于含有多个分子的文件（比如multi_mol2以及sdf），需要将之分解为多个单分子文件。有以下方法实现：
i)使用babel：
$babel -imol2 specs_p0.0.mol2 -o specs.mol2 -split
将生成specs0001.mol2,specs0002.mol2 ...
ii)使用splitmol（DOCK软件包中程序，只能对mol2和pdb文件进行分解）：
$mkdir specs_p0.0
$cd specs_p0.0
$splitmol -mol2 ../specs_p0.0.mol2
将生成split00001,split00002,...
2、利用上面命令分解得到的单分子文件，外面的名称和分子内部的名称不符，可以通过编写教本批量重命名：
适用于ZINC数据库中的mol2数据集：
for i in split*
do
mv $i $(grep ZINC $i).mol2
done
提取包含ZINC所在行的字段，在ZINC数据库中，恰巧为分子的名称。最后，全部命名为ZINC××××××.mol2。
3、AutoDock官方发布版本不能用于并行运算（有的课题组开发了并行版本，但并未发布）。为了提高效率，可以将几十万个分子文件分为几百个部分，每个1000个分子。
将specs数据库的20万个分子文件按照1000个分子为一组（每组为一个目录），分为20个目录
$mkdir specs_0.0_01
$mv specs0.0.mol2 ../specs_0.0_01 利用通配符将前999个分子转移入新文件夹，实现分装目的。
$ls ../specs_0.0_01 | wc [-w] 查看文件个数
4、在AutoDock输入文件的准备中，使用了很多外部程序，比如Sybyl(Tripos Ltd.)。小分子三维结构数据库都有在线免费的。如果想要构建自己的化合物库，可以用ChemDraw（商业）或ISIS/Draw（免费）画好后，存为二维格式（推荐用mol格式），然后用Corina等软件转化成三维结构。
5、AutoDock可能会发生堆栈溢出的情况，错误提示为“Segment Fault”，可以在系统的配置文件（通常为/etc/bashrc)中添加如下行：
unlimit stacksize
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
※※※※※※ Autodock分子对接流程 ※※※※※※
1、受体文件准备：
Atudock运算的受体加极性氢（polar hydrogens, essential_only hydrogens),加Kollman United点电荷。首先把PDB蛋白文件中除蛋白以外的所有物质全部删除（可以在Sybyl中操作），保存未macro.pdb。加氢加电荷通过自由软件pmol2q可以很好实现：
$pmol2q macro.pdb macro.mol2
选择加氢，加电荷，pmol2q可以输出mol2，pdbq以及pdbqs文件。
在Autodock3中，受体的格式为pdbqs，可以通过如下方法准备：
i）通过AutoDock软件包中程序：
SGI Irix下运行：
%mol2fftopdbq macro.mol2 > macro.pdbq
%addsol macro.pdbq macro.pqdbs
以上两个命令可以通过一个命令实现：
%mol2topdbqs macro.mol2
Linux下运行：
$mol2topdbq macro.mol2 > macro.pdbq
$addsol macro.pdbq macro.pqdbs
ii）通过ADT中python脚本：
$pythonsh $adtpy/prepare_receptor.py -r macro.mol2 [-o macro.pdbqs] -C
注：
-C:保留原有电荷
在Autodock4中，受体的格式为pdbqt，与PDB格式类似，包含点电荷 ('Q')和AutoDock4原子类型('T')。可以通过如下方法准备：
$pythonsh $adtpy/prepare_receptor4.py -r macro.mol2 [-o macro.pdbqs] -C
注：
也可以由python脚本中的pdbqs_to_pdbqt.py来将autodock3的pdbqs文件转化。
-o:不加则自动生成文件名相对应的pdbqt文件。
-C:保留原有电荷，否则添加Gaster电荷
2、配体文件准备：
在Autodock3中，配体的格式为pdbq，可以通过如下方法准备：
i）通过AutoDock软件包中程序：
%deftors lig.mol2
ii）通过ADT中python脚本：
$pythonsh $adtpy/prepare_ligand.py -l lig.mol2 [-o lig.pdbq] -C
注：
-C:保留原有电荷
-o：不加则自动生成文件名相对应的pdbq文件。
在Autodock4中，配体的格式为pdbqt.可以通过如下方法准备：
$pythonsh $adtpy/prepare_ligand4.py -l lig.mol2 [-o lig.pdbqt] -C
注：
也可以由python脚本中的pdbq_to_pdbqt.py来将autodock3的gpf文件转化:
$pythonsh pdbq_to_pdbqt -s lig [-o lig.pdbqt]
3、准备GPF(grid parameter file)文件：
Autodock3：
i）通过AutoDock软件包中程序：
$mkgpf4 lig.pdbq macro.pdbqs
注：
程序包中有mkgpf3，不过mkgpf4是其修正后的版本，建议使用mkgpf4，这个与autodock4没有关系。
ii）通过ADT中python脚本：
$pythonsh $adtpy/prepare_gpf.py -r macro.pdbqs -l lig.pdbq
Autodock4(只能通过ADT中python脚本)：
$pythonsh $adtpy/prepare_gpf4.py -r macro.pdbqt -l lig.pdbqt [-p gridcenter='x y z'] [-o lig.macro.gpf]
注：也可以由python脚本中的gpf3_to_gpf4.py来将autodock3的gpf文件转化。
-o:输出文件名
-p gridcenter='x y z' 注意用单引号
注：
1、准备gpf文件时，最好通过ADT图形界面进行，这是因为：
i）gpf文件中活性位点（center grid）默认根据配体位置定义（auto），需要修改为晶体结构中配体的具体坐标值；当然，在命令行中使用-p gridcenter='x y z'参数也可以定义。
ii）盒子（box）大小在图形界面中可以很方便的调整。
iii）在虚拟筛选中需要添加更多的院子类型，在图形界面中操作非常方便，可以通过“Grid”菜单中的“Directly Set”来输入需要的原子类型。
2、配体的原子类型说明：
AD3和AD4的原子类型写法差异很大，而且，在gpf文件中，表达格式差异也很大。在AD4中，各种原子类型写在一行中，每种类型以空格间隔，但是在AD3中，须对每种原子类型的具体参数（作用半径）进行定义，比较繁琐。下面是两种版本中原子类型的写法：
－－－－－－－－－－－－－－－－－－－－－－－
AD3:
C carbon (aliphatic)
A carbon (aromatic)
N nitrogen
O oxygen
P phosphorus
S sulphur
H hydrogen
f iron
F fluorine
c chlorine
b bromine
I iodine
－－－－－－－－－－－－－－－－－－－－－－－
AD4（共22种，下表中有7中重复，在autodocksuite-4.0.0/src/autodock-4.0.0/AD4_parameters.dat文件）
H Non H-bonding Hydrogen
HD* Donor 1 H-bond Hydrogen
HS Donor S Spherical Hydrogen
C* Non H-bonding Aliphatic Carbon
A* Non H-bonding Aromatic Carbon
N* Non H-bonding Nitrogen
NA* Acceptor 1 H-bond Nitrogen
NS Acceptor S Spherical Nitrogen
OA* Acceptor 2 H-bonds Oxygen
OS Acceptor S Spherical Oxygen
F Non H-bonding Fluorine
Mg Non H-bonding Magnesium
MG Non H-bonding Magnesium
P Non H-bonding Phosphorus
SA* Acceptor 2 H-bonds Sulphur
S Non H-bonding Sulphur
Cl Non H-bonding Chlorine
CL Non H-bonding Chlorine
Ca Non H-bonding Calcium
CA Non H-bonding Calcium
Mn Non H-bonding Manganese
MN Non H-bonding Manganese
Fe Non H-bonding Iron
FE Non H-bonding Iron
Zn Non H-bonding Zinc
ZN Non H-bonding Zinc
Br Non H-bonding Bromine
BR Non H-bonding Bromine
I Non H-bonding Iodine
－－－－－－－－－－－－－－－－－－－－－－－
注：
*默认在gpf中存在的原子类型
$bable -imol2 lig.mol2 -obox box
0.375
4、准备DPF(docking parameter file)文件：
Autodock3：
i）通过AutoDock软件包中程序：
$mkdpf3 lig.pdbq macro.pdbqs
ii）通过ADT中python脚本：
$pythonsh $adtpy/prepare_dpf.py -r macro.pdbqs -l lig.pdbq
Autodock4(只能通过ADT中python脚本)：
$pythonsh $adtpy/prepare_dpf4.py -r macro.pdbqs -l lig.pdbqt -p ga_num_evals=25000000 -p ga_run=100 -p ga_pop_size=300
注：
也可以由python脚本中的dpf3_to_dpf4.py来将autodock3的gpf文件转化。
-p parameter=value：输入参数值
注：
此处生成的dpf需要编辑其中的参数，对遗传算法的精度进行设置。dpf文件中参数设置：
其中，"ga_num_evals"和"ga_num_evaluations"用于限定程序收敛条件，决定运算费时长短。
'ga_' 起始行参数：genetic algorithm
'sw_' 起始行参数：Solis and Wets local search method
ga_num_evals=25000000 Torsions<10，设定为2.5e5 － 2.5e7；Torsions>10，设定为2.5e7（默认2.5e6)
ga_run=100 最好在50以上（默认10）
ga_pop_size=300（默认150）
5、通过AutoGrid计算格点能量：
$autogrid3/4 -p macro.gpf -l macro.glg
注：
AutoGrid默认计算8种原子类型，可以通过修改源代码，使其计算更多院子类型的格点能量。
6、通过AutoDock进行分子对接：
$autodock3/4 -p lig.macro.gpf -l lig.macro.glg
注：
和AutoGrid一样，AutoDock默认计算8种原子类型，可以通过修改源代码，使其计算更多院子类型的格点能量。
7、提取结果：
i）通过AutoDock软件包中程序：
$get-docked lig.macro.dlg
注：
生成lig.macro.dlg.pdb文件
ii）通过ADT中python脚本：
pythonsh $adtpy/summarize_results.py
pythonsh $adtpy/summarize_results4.py
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
※※※※※※ AutoDock结果分析 ※※※※※※
AutoDock产生的构象和评分值（单位kcal/mol，AD4提供了RMSD值）保存为PDB格式，这种格式没有mol2格式操作方便。
1、编写脚本提取运算结果的评分值并对评分值进行排序
adscore3（适用于AutoDock3）：提取PDB文件中的评分值并对评分值进行排序（需要用get-docked命令生成pdb文件）
#####################
egrep 'Run|Final Docked Energy' *.pdb >> adscore3.tmp
#perl -pi -e 's/\nUSER Final/ Final/;' adscore.tmp
#sed "s/USER //g" adscore.tmp | sort -r +8 score.tmp > adscore.log
vi -e adscore3.tmp << EOF
:%s/USER //g
:%s/\nFinal/ Final/g
:%s/\[=(1)+(2)\]//g
:wq
EOF
sort -r +7 adscore3.tmp* > adscore3.log
rm -f *.tmp*
cat adscore3.log
#####################
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
adgrep4（适用于AutoDock4）：提取DLG文件中的评分值
#####################
#!/bin/bash
egrep 'DOCKED: USER Run|Cluster Rank|RMSD from reference structure|USER Estimated' *.dlg
#####################
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
adcluster4（适用于AutoDock4）：提取DLG文件中的Cluster Rank部分
#####################
#!/bin/bash
egrep 'Cluster Rank' -A5 -B2 *.dlg > adcluster4.tmp
sed "s/\[=(1)+(2)+(3)-(4)\]//g" adcluster4.tmp > adcluster4.log
cat adcluster4.log
#####################
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
adscore4（适用于AutoDock4）：提取DLG文件中的能量值和Ki值并根据能量值进行排序
#####################
#!/bin/bash
egrep 'USER Run|USER Estimated' *.pdb > adscore4.tmp1
egrep '.dpf' -v adscore4.tmp1 > adscore4.tmp2
sed "s/USER //g" adscore4.tmp2 | egrep 'USER' -v > adscore4.tmp3
vi -e adscore4.tmp3 << EOF
:%s/\nEstimated Free Energy of Binding/, Estimated Free Energy of Binding/g
:%s/ \[=(1)+(2)+(3)-(4)\]\nEstimated Inhibition Constant,/, /g
:%s/\[Temperature = 298.15 K\]//g
:wq
EOF
sort -r +9 adscore4.tmp3 > adscore4.log
rm -f *tmp*
cat adscore4.log
#####################
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
2、用于分析AutoDock运算结果中受体－配体相互作用的图形化软件中，理论上最好用的应该是ADT（我没有用过，在Linux下非常慢，可能是由于从工作站上调用，显示延迟的缘故，而在Windows机器上总是无故退出）。我用过的最好用的是Vida（Openeye Ltd.），以前提供2个月的试用版，我还申请过1年的学术用户版，但一年后License过期，就没的用了。当然，把系统时间调整回到License有效期内，还是可以使用（Openeye公司的软件都是精品呵）。Vida可以同时读入受体、参考分子（这里是晶体结构中的配体）和对接得到的构象，显示非常好看的ribbon形态和氢键，作图效果一流，不过操作起来不是很人性化。
AutoDock的手册示例中使用Sybyl软件对结果进行分析，不过Sybyl每次只能从PDB文件中读入1个分子，还是不很方便。
没有尝试UCSF chimera是不是可以，不过，看其对DOCK处理的得心应手，应该差不多。
－－－－－－－－－－
非常感谢
有时间一定对照着学习学习. 太好了，对这个我做了几次，但是不是太熟悉，准备照着学习学习，谢谢fwangj非常感谢
兄弟，太感谢啦。我在巴黎这边正在学习这个东西，但是因为半路出家，还是比较茫然的。兄弟的这篇文章帮我大忙了。真是太感谢了。:) 太好了，照着学习学习 有没有简单易学的教材啊，多谢 谢谢分享。不错的教程回复 1# fwangj 的帖子
谢谢分享,太有用了 我也是半路出家 挺茫然的 照着做做页:
[1]