Back

利用OpenBabel生成多个构像

此文感谢建民兄的指导,才让我学习到这么有用的工具

起因

有许多软件例如Dock 6.8,对于小分子的柔性对接支持并不那么好,这就导致了因为初始构像的不正确而导致了对接结果的不正确和不可靠,为了解决这个问题一个方法就是进行构像的生成,再进行多对接,整理得分排序。当然还有一个应用就是用于3D药效团的使用。网上有许多在线网站可以进行这方面的工作,但是最好用的我感觉还是openbabel神器!

利用openbabel生成多构像

Open Babel提供了两种构像生成的算法代码:

  1. Confab: 一个成体系的构像生成器,有序的生成所有的低能量构像
  2. 遗传算法: 是一个随机构像产生器,能够根据RMSD或者能量的不同达到优化的目的

1. 遗传算法

了解算法的命令信息可以使用:obabel -L conformer查看。虽然查看是说的构像搜寻,但是可以使用--writeconformers对生成的构像进行保存。可以查看如下的生成30个构像的例子:

1
obabel startingConformer.mol -O ga_conformers.sdf --conformer --nconf 30  --score rmsd --writeconformers

其中--score rmsd可以不写,因为默认就是使用这种方式

当然生成的构像在1个文件里,我们有时候需要split,当然首先想到的还是obabel,其进行拆分相当于高射炮打蚊子–相当好用呀(<-原谅我在这里乱用歇后语),示例如下:

1
obabel ga_conformers.sdf -O new.mol -m

我自己的输出结果如下:

1
2
32 molecules converted
32 files output. The first is molecule/new1.mol2

2. Confab

Confab 生成分子的低能量构像分子集合。设置命令为--confab,同样可以使用obabel -L confab进行查看了解。 使用简单的设置操作如下:

1
obabel <inputfile> -O <outputfile> --confab [confab options]

其还有一个额外的应用就是confabreport format,如下:

1
obabel <inputfile> [-O <outputfile>] -o confabreport -xf <reference_file> [-xr <rmsd>]

主要是为了和初始构像进行比较。

下面是其中一个例子,首先生成100K个构像:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
> obabel bostrom.sdf -O confs.sdf --confab --conf 100000

**Starting Confab 1.1.0
**To support, cite Journal of Cheminformatics, 2011, 3, 8.
..Input format = sdf
..Output format = sdf
..RMSD cutoff = 0.5
..Energy cutoff = 50
..Conformer cutoff = 1000000
..Write input conformation? False
..Verbose? False

**Molecule 1
..title = 1a28_STR_1_A_1__C__
..number of rotatable bonds = 1
..tot conformations = 12
..tot confs tested = 12
..below energy threshold = 10
..generated 3 conformers

... etc, etc

0 molecules converted

统计多少个构像和初始构像相比RMSD小于1.0A

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
> obabel confs.sdf -oconfabreport -xf bostrom.sdf -xr 1.0

**Generating Confab Report
..Reference file = bostrom.sdf
..Conformer file = confs.sdf

..Molecule 1
..title = 1a28_STR_1_A_1__C__
..number of confs = 3
..minimum rmsd = 0.0644801
..confs less than cutoffs: 0.2 0.5 1 1.5 2 3 4 100
..1 1 3 3 3 3 3 3
..cutoff (1) passed =  Yes

... etc, etc

**Summary
..number of molecules = 36
..less than cutoff(1) = 35
52271 molecules converted

我做的一个的简单应用(批量读取文件对接),不提供部分代码:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#-*-coding:utf-8-*-

from ucsfdock.docksh import PreDock
from ucsfdock.dockfile import Dock
import os

#计时器工具
import time
time_start=time.time()

DOCKHOME='/home/kangsgo/install/dock6/bin/'
#run=PreDock(DOCKHOME=DOCKHOME)
#run.onerun()
dock=Dock(DOCKHOME=DOCKHOME)
#遍历molecule目录下所有的文件
for (dirpath,dirnames,filenames) in os.walk('molecule'):
    for i in filenames:
        print('开始执行'+i+'文件')
        dock.dock(mode=2,ligand_atom_file='molecule/'+i,ligand_outfile_prefix='out1/'+i[:-5])
#dock.dock(mode=2,ligand_atom_file='2.mol2')
print("done!")
time_end=time.time()
print("耗时为"+str(time_end-time_start)+" 秒")

参考资料:

  1. Generate multiple conformers
  2. obabel and babel