LeDock虚拟高通量筛选结果分析脚本
本帖最后由 fireflying 于 2015-6-2 01:57 编辑分子对接的结合能计算一般来说很不靠谱，因为：1）分子对接能量计算采用格点方式，格点间隔大小对结果精度有影响，譬如AutoDock默认格点间隔为0.375埃，也就是说对接精度的误差至少在0.375埃左右。更为精细的格点间隔设置可以提高一点精度，不过内存需求也会大为增加。2）构象搜索采用随机搜索算法，搜索有可能提前收敛，换句话说构象不一定是（全局）最优。3）能量方程无论是基于物理的AutoDock或者机器学习的Vina或者其它对接程序的方程基本都不太靠谱，是完备力场方程的阉割或者简化版本。4）去溶剂化效果的计算基本都很糟糕。5）对小分子可旋转键的旋转能垒基本没有考虑。6)某些对接程序如Autodock或Vina没有考虑非极性氢。基于以上的原因，分子对接的能量误差通常都在2 Kcal/mol左右，这个误差有多大呢，基本上就是说对接能量差不多没有计算上的参考意义。如果是同一个分子骨架系列衍生物的对接（譬如药化的日常工作），分子对接的打分也许有一定帮助。至于虚拟高通量筛选，化合物数量越大，基于分子对接能量的排序就越糟糕，更多的是用来设定一个阈值，譬如如下文中的用法：
http://www.sciencedirect.com/sci ... i/S0960894X14001279
要比较准确的预测结合能，分子对接后的第一步是要做分子能量最小化（energy minimization），然后采用更为精确的结合能方程基于完备的力场方程、参数等进行计算。基于这个原因，LeDock在设计之初并没有考虑对筛选结果进行分析。考虑到某些朋友不太熟悉脚本，本文附件提供一个分析虚拟高通量对接结果的文件，可方便的得到每个分子的对接能量以及配体结合效率（Ligand Efficiency）。建议大家在虚拟高通量筛选中更多使用Ligand Efficiency来进行排序，具体原因可以参考以下文章（文章中通过虚拟筛选发现的化合物与靶点蛋白的复合物结构见PDB code： 4PCI、4PCE）：
http://www.sciencedirect.com/sci ... i/S0960894X14003539
如果有条件使用CHARMM或者其它程序可以对对接构象进行能量最小化的朋友，建议采用结合LeWater的打分函数，或者采用LeWater作为一个filter，预先去除掉不合理的对接构象（分子）。具体可参考以下文章（PDB code： 4P4C）：
http://journals.plos.org/plosone ... ournal.pone.0019923
虽然笔者不建议用对接打分来进行排序，不过需要指出的是LeDock在虚拟筛选中预测化合物的结合构象是相当准确的，譬如上述工作中，LeDock预测的结合构象跟最终晶体结构中解析出来的配体结合构象是完全一致的。
附件使用方式
1. g++ ledock_anal.cpp -o ledock_anal -O2
2. chmod +x ledock_anal.csh
3. 将编译好的ledock_anal以及ledock_anal.csh添加到路径里面
4. cd到包含所有dok文件的目录，运行ledock_anal.csh。结果在docking_summary.txt文件中。如果要对对接结果进行排序，使用以下命令
5. less docking_summary.txt | sort -k2n 根据结合能排序，或者less docking_summary.txt | sort -k3n根据Ligand efficiency排序。
本帖最后由 fireflying 于 2015-6-2 01:51 编辑
ZhouDeng 发表于 2015-6-1 15:29
我阅读了您的《Hydrogen Bonding Penalty upon Ligand Binding》这篇文章。有几个问题想请教一下：
1）如何 ...
简单介绍一下LeWater可以计算什么，以及原理。LeWater评价在蛋白配体结合界面上极性原子的氢键（采用基于DFT的半经验方程）以及水化状态，综合给出一个指数（Hydrogen bonding penalty），这个指数越小越好，通常小于2，高通量虚拟筛选可以用1作为阈值。Hydrogen bonding penalty指标的计算对对接构象的精度要求很高（基于QM的计算嘛{:soso_e128:}），所以对接构象最好经过能量最小化（energy minimizaiton），否则不建议用。退而求其次，LeWater可以自动检测配体跟蛋白的氢键情况，基于DFT方程评价每个氢键的状态，数值从0到1。0表示没有氢键，1表示氢键最优。具体的参数设置如下：
Pocket3 36 11 39 24 52第二行数字对应xmin、xmax、ymin、ymax、zmin、zmax，如果你使用过LeDock，对这个应该是比较熟悉的。不过在计算Hydrogen bonding penalty时候口袋要定义的比LeDock更大，建议是在LeDock口袋基础上：xmin = xmin（LeDock）-6，xmax = xmax（LeDock）+6。
Donor919 N HN1046 N HN此处定义蛋白上的给体，氢键给体包括氢键给体的重原子（N或O）以及氢原子。第一列数字表示氨基酸序列号，后面是氢键的重原子名称，然后是氢原子名称。这些名称须跟PDB文件一致。可以用pymol查看，点击某个氨基酸残基，上方文本框会显示氨基酸残基编号，切换到编辑状态，点击某个原子则会显示原子名。
Acceptor917 O885 OE2 0此处定义蛋白上的受体。方法同上。第三列的数字是在计算Hydrogen bonding penalty时候的权重，缺省为1。计算氢键并不需要用到这个参数。0意味着该原子的penalty为0。
以上是参数文件的准备，譬如保存为lewater.in。因为氢键需要考虑氢原子，所以蛋白必须是加好氢的，可以用任何加氢软件生成蛋白结构，包括LePro，此处假设蛋白文件为pro.pdb。
接下来计算就很简单了。ls *dock*.pdb > list #该命令生成一个对接构象的列表lewater lewater.in pro.pdb list
大功告成！输出文件名为filter，查看输出：cat filter4ag8_lig_dock001.pdb 1.275P: 0.110L: 1.165Dpro: 0.9 1.0Apro: 1.0 1.0第一列为结合构象的文件名，第二列为Hydrogen bonding penalty，P：0.110对应来自蛋白的penalty数值，L：1.165对应来自于配体的penalty值。Dpro:后面的数值对应参数文件里面在相应位置形成氢键的情况，譬如0.9表示配体与蛋白919 N HN形成氢键，氢键综合评价为0.9（1最高），如果数值为零则表示配体没有与919 N HN形成氢键，以此类推。
作filter的时候用awk命令，譬如我们要求配体必须跟917 O形成氢键并且指数大于0.5，则cat filter | awk '$11>0.5'#注意该氢键对应第11列的1.0上述命令给出所有符合要求的配体结合构象。
接下来介绍一下LeWater的一个极富特色的功能，生成WaterMap。ls lig_dock001.pdb >list #计算WaterMap的时候list里面只能放一个结合构象lewater lewater.in pro.pdb list -omap生成一个文件名为omap.pdb文件，用pymol加载，我们就得到一个在蛋白配体结合界面的WaterMap，如下图。红色球表示晶体水，绿色小球表示预测的WaterMap。图中每个晶体水周围都有绿色小球球，表示我们成功的预测了结晶水的位置，绿色球球的多少表示该结晶水的熵。球球越多，说明该结晶水越Happy。结构优化的一个重要方向就是通过引入基团来取代结晶水，好啦，不多说了。Cheers！
谢谢。对虚拟结果的分析在某些时候是非常重要的。 您好。我在尝试的过程中没有成功。编译出来的ledock_anal双击打开后显示如下
Last login: Sat May9 11:14:56 on ttys001dengzhoudeMacBook-Pro:~ dengzhou$ /Users/dengzhou/ledock_anal ; exit;Segmentation fault: 11logout
[进程已完成]
我将编译好的这个文件连同另外两个文件一起拷贝到了我的dok文件夹下，运行ledock_anal.csh，得到的txt文件里没有结果。
运行窗口显示如下：dengzhoudeMacBook-Pro:result dengzhou$ /Users/dengzhou/DZ/***\ Data/***\ docking/***/result/ledock_anal.csh ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.ledock_anal: Command not found.
ledock_anal: Command not found.
我使用的是OS X 10.10.3。您能帮我看看是什么问题吗？
ZhouDeng 发表于 2015-5-9 11:23
您好。我在尝试的过程中没有成功。编译出来的ledock_anal双击打开后显示如下
Last login: Sat May9 11:14 ...
分析脚本需要在terminal（终端）里面通过命令行来运行。既然你已经编译好了，估计对terminal有一定的了解了。接下来就很简单了。
运行错误信息是系统没有找到ledock_anal这个可执行文件，这是因为1）系统默认路径里面没有ledock_anal；2）当前目录没有加到默认的搜索路径里面。
解决方法：1）把ledock_anal复制到系统默认的路径里面或者把当前目录添加到path里面；或者2）修改ledock_anal.csh,如下
./ledock_anal $f >>docking_summary.txt
./代表当前目录。
OS X系统的话，建议用LeDock的命令行版本，可以方便的做高通量筛选，而不需要一个个分子的手动对接。命令行版本的输入文件中需要设定一个待对接的所有小分子的列表文件。譬如有1000个分子需要对接，用 ls *mol2 > ligands。然后运行ledock dock.in &。一个命令即可，ledock会自动把1000个分子依次对接好。
具体问题可以加入QQ群：325901617 我阅读了您的《Hydrogen Bonding Penalty upon Ligand Binding》这篇文章。有几个问题想请教一下：
1）如何选取蛋白中donor和acceptor，是利用原始的结晶结构，还是lepro处理之后的结构文件？
2）我已经利用lewater做出了一个filter，如何利用这个filter去筛选我的结合构象呢？ :loveliness::loveliness:高大上
页:
[1]