Hi dear,
I would like to perform covalent docking using the flexible side chain
method as recently described in the paper "Covalent docking using autodock:
Two-point attractor and flexible side chain methods".
I have used the script provided in
http://autodock.scripps.edu/resources/covalentdocking, and
I have followed the procedure in the README text (attached).
But, I got the following openbabel error.
~ /Desktop/adCovalentDockResidue/adcovalent/prepareCovalent.py --ligand
ligand.mol2 \
--ligindices 1,2\
--receptor 3upo_protein.pdb\
--residue B:SER222\
--outputfile ligcovalent.pdb
==============================
*** Open Babel Error in ReadFile
Cannot read from 3upo_protein.pdb
[ ERROR ] Fail to read molecule "3upo_protein.pdb"
I guess openbabel (ver 2.4.1) (pybel) has been properly installed in Mac OS
X 10.10.5, and imported in preparecovalent.py (attached).
What is the reason of this openbabel error, and how can I fix this error?
Any suggestions or comments are highly appreciated.
Thank you.
----------------------------------------------------------------------------------
README.txt
This tutorial describes how to apply the 'flexible side chain' covalent
docking as described in the paper:
"Covalent docking using autodock: Two-point attractor and flexible side
chain methods"
http://www.ncbi.nlm.nih.gov/pubmed/26103917
If you use this method, please cite the paper.
Happy dockings!
==============
REQUIREMENTS
==============
- Python >=2.6
- OpenBabel Python bindings
(On Debian/Ubuntu derivatives, it can be installed with "apt-get install
python-openbabel")
- Basic experience with docking with AutoDock
==============
INSTALLATION
==============
- Copy the directory 'adcovalent' in a location that is accessible from the
user (i.e., '~/local/adcovalent').
==============
USAGE
==============
The preparation code is used to generate the geometry of the ligand
covalently bound to the
residue. Once the ligand initial coordinate is generated, the PDBQT files of
the receptor and
flexible ligand need to be generated using prepare_flexreceptor4.py script
in AutoDockTools.
The directory '3upo_test' contains an example case of a ligand ready to
processed.
The following examples assume a terminal opened in that location.
1. Generate the alignment
--------------------------
The input ligand structure must be modeled with a portion of the alkylated
residue already present,
i.e. for a ligand bound to a serine:
LIGAND-OG-CB-CA-RESIDUE [MOL1]
the ligand structure to process should be:
LIGAND-O-C [MOL2]
where C and O are the two atoms that are going to be used for the alignment.
More in general:
0. Initial configuration
B
/
[L]----A I---J--[R]
/
X
1. Translate [L] to overlap A->I
B
/
[L]----A -> I---J--[R]
/
X
B
/
[L]----AI---J--[R]
/
X
2. Rotate [L] around the normal of plane (B,A,J) to overlap B->J
[L]
\
AI--BJ--[R]
/
X
The alignment could be skipped if the ligand is already attached to the
receptor structure
(i.e., a re-docking experiment from a PDB is going to be perfoemd).
To perform the alignment, a receptor file and a residue name (Chain:ResNum)
must
be specified, then a ligand file and an alignment criterion.
The ligand alignment can be defined in two ways. The first is by specifying
the atom indices in the ligand file (i.e., first and second atoms):
python ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
--ligindices 1,2\
--receptor 3upo_protein.pdb\
--residue B:SER222\
--outputfile ligcovalent.pdb
or, in a more general approach, by specifying a SMARTS pattern and the
indices
of the alignment atoms in the pattern itself:
python ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
--ligindices 4,3\
--ligsmart "C(=O)-O-C"\
--receptor 3upo_protein.pdb\
--residue B:SER222\
--outputfile ligcovalent.pdb
2. Generate PDBQT files
--------------------------
The standard PDBQT files for AutoDock need to be generated for both the
receptor structure. The standard ADT scripts
are going to be used, assuming that $MGLROOT is the path where ADT has been
installed.
*** NOTE: if a covalent ligand is already bound to the residue in the
protein structure, it must be removed! ***
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r
3upo_protein.pdb -A hydrogens
[ this generates the file "3upo_protein.pdbqt" ]
and the covalent ligand that has been aligned:
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r
ligcovalent.pdb
[ this generates the file "ligcovalent.pdbqt" ]
3. Generate flexible/rigid PDBQT files
---------------------------------------
The PDBQT files are going to be used to generate the rigid and flexible
components that are going to be used
for the docking.
First, the receptor is processed to extract the rigid part by specifying
which residue is going to be made flexible:
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r
3upo_protein.pdbqt -s 3upo_protein:B:SER222
[ this generates the files "3upo_protein_rigid.pdbqt" and
"3upo_protein_flex.pdbqt" ]
The same thing is going to be done for the processed ligand:
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r
ligcovalent.pdbqt -s ligcovalent:B:SER222
[ this generates the files "ligcovalent_rigid.pdbqt" and
"ligcovalent_flex.pdbqt" ]
NOTE: if the ligand has not been prepared in step 1, make sure that all
atoms in the ligand have the poper residue id
corresponding to the residue modified (i.e., Ser222A).
4. Generate GPF and DPF files
------------------------------
The parameter files for the actual calculation are going to be generated.
The GPF for AutoGrid is generated with:
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_gpf4.py -r
3upo_protein_rigid.pdbqt\
-x ligcovalent_flex.pdbqt\
-l ligcovalent_flex.pdbqt\
-y -I 20\
-o 3upo_priotein.gpf
The DPF for AutoDock is generated with:
$MGLROOT/bin/pythonsh
$MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_dpf4.py -r
3upo_protein_rigid.pdbqt\
-x ligcovalent_flex.pdbqt\
-l ligcovalent_flex.pdbqt\
-o ligcovalent_3upo_protein.dpf\
-p move='empty'
This command instructs the script to dock the covalent ligand as a flexible
residue and ignore any 'true' ligand ("move='empty'"). To do this,
an empty file must be created, and it can be done with the following Unix
command:
touch empty
Also, in order to include the ligand atoms in the calculation of RMSD and
clustering, the following parameters need to be added:
rmsdatoms all
which will include flexible residue (and ligand) atoms in the RMSD
calculation.
Finally, the DPF file must be manually edited to set the appropriate energy
model so that the docking score corresponds to the interaction between
the flexible residue (the ligand) and the rigid receptor. For this, the
entry:
unbound_model bound
must be replaced with:
unbound_energy 0.0
5. Run AutoGrid and AutoDock
------------------------------
Both programs can be executed using the standard procedure:
autogrid4 -p 3upo_priotein.gpf -l 3upo_priotein.glg
autodock4 -p ligcovalent_3upo_protein.dpf -l
ligcovalent_3upo_protein.dlg
The output generated at each step of the process is available in the
directory '3upo_test/output'.
------------------------------------------------------------------------------------------------------------
prepareCovalent.py
#!/usr/bin/env python
#
# Covalent Docking preparation for AutoDock
#
# v.1.0c Stefano Forli
#
# Copyright 2014, Molecular Graphics Lab
# The Scripps Research Institute
# _
# (,) T h e
# _/
# (.) S c r i p p s
# '\_
# (,) R e s e a r c h
# ./'
# ( ) I n s t i t u t e
# "
#
#################################################################################
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>
#
import sys, os, math
import pybel
ob = pybel.ob
import argparse
import numpy as np
from math import sqrt, sin, cos, acos, degrees
import ResidueProfiler
from CopyMol import MoleculeDuplicator
# TODO improve the method to build the bond directly, calculating
# bond lenght and its placement
#### Debug functions
def makePdb(coord, keyw = "ATOM ", at_index = 1, res_index = 1, atype =
'X', elem = None,
res = "CNT", chain ="Z", bfactor = 10,pcharge = 0.0):
if not elem: elem = atype
# padding bfactor
bfactor = "%2.2f" % bfactor
if len(bfactor.split(".")[0]) == 1:
bfactor = " "+bfactor
if len(atype) == 1:
atype = atype + " "
#atom = "%s%5d %2s %3s %1s%4d %8.3f%8.3f%8.3f 1.00 %02.2f %8.3f
%1s" % (keyw,
# SOURCE: http://cupnet.net/pdb-format/
"%-6s%5d %4s%1s%3s %1s%4d%1s %8.3f%8.3f%8.3f%6.2f%6.2f
%2s%2s"
atom = "%s%5d %2s %3s %1s%4d %8.3f%8.3f%8.3f 1.00 %s %8.3f %1s" %
(keyw,
at_index, elem, res, chain, res_index,
coord[0], coord[1], coord[2], bfactor, pcharge, atype)
return atom
def writeList(filename, inlist, mode = 'w', addNewLine = False):
if addNewLine: nl = "\n"
else: nl = ""
fp = open(filename, mode)
for i in inlist:
fp.write(str(i)+nl)
fp.close()
def vector(p1 , p2 = None, norm = 0): # TODO use Numpy?
if not p2 == None:
vec = np.array([p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]],'f')
else:
vec = np.array([p1[0], p1[1], p1[2] ], 'f' )
if norm:
return normalize(vec)
else:
return vec
def dot(vector1, vector2): # TODO remove and use Numpy
dot_product = 0.
for i in range(0, len(vector1)):
dot_product += (vector1[i] * vector2[i])
return dot_product
def vecAngle(v1, v2, rad=1): # TODO remove and use Numpy?
angle = dot(normalize(v1), normalize(v2))
try:
if rad:
return acos(angle)
else:
return degrees(acos(angle))
except:
print "#vecAngle> CHECK TrottNormalization"
return 0
def norm(A): # TODO use Numpy
"Return vector norm"
return math.sqrt(sum(A*A))
def normalize(A): # TODO use Numpy
"Normalize the Vector"
return A/norm(A)
def calcPlane(p1, p2, p3):
# returns the plane containing the 3 input points
v12 = vector(p1,p2)
v13 = vector(p3,p2)
return normalize(np.cross(v12, v13))
def rotatePoint(pt,m,ax):
"""
Rotate a point applied in m around a pivot ax ?
pt = point that is rotated
ax = vector around wich rotation is performed
"""
# point
x=pt[0]
y=pt[1]
z=pt[2]
# rotation pivot
u=ax[0]
v=ax[1]
w=ax[2]
ux=u*x
uy=u*y
uz=u*z
vx=v*x
vy=v*y
vz=v*z
wx=w*x
wy=w*y
wz=w*z
sa=sin(ax[3])
ca=cos(ax[3])
p0 =(u*(ux+vy+wz)+(x*(v*v+w*w)-u*(vy+wz))*ca+(-wy+vz)*sa)+ m[0]
p1=(v*(ux+vy+wz)+(y*(u*u+w*w)-v*(ux+wz))*ca+(wx-uz)*sa)+ m[1]
p2=(w*(ux+vy+wz)+(z*(u*u+v*v)-w*(ux+vy))*ca+(-vx+uy)*sa)+ m[2]
return np.array([ p0, p1, p2])
class CovalentDockingMaker:
""" This class generate the input for covalent dockings.
Required input are :
- ligand structure (any OB-supported 3D format)
- receptor structure (any OB-supported 3D format)
Ligand and receptor structures should share at least 2 atoms
to be overlapped, e.g.:
Cys-Ca-Cb-S Cb-S-Lig
Alignment can be performed using one of the following criteria:
- indices: two indices pairs for ligand and receptor
- SMARTS: a pattern and two atom indices must be defied for the
ligand;
if nothing is defined for the receptor (default), the
residue
type will be used to pick its strongst nucleophyle
atom,
i.e. Cys-SG, Ser/Thr-OG, Lys-NZ.
- bond: a bond for each structure is provided and indices
calculated from them; bonds must be terminal bonds
(not counting H's)
(not working yet)
"""
def __init__(self,
lig, # ligand openbabel molecule
rec, # receptor openbabel molecule
resName,# residue id to be modifled (i.e. A:THR276)
smartsLig=None, smartsRec=None, # use SMARTS patterns
# patterns should be aligned
smartsIdxLig=(), smartsIdxRec=(), # align using atomIdx from
SMARTS patterns
bondIdxLig=None, bondIdxRec=None, # use bonds
indicesLig=None, indicesRec=None, # use atom indices in both
genFullRec = False, # add the covalent ligand to
the receptor
# as modified residue; this is
*NOT*
# recommended, because it can
lead to
# spurious bonds between the
mod res
# and the rest of the receptor
genLigand = False,
auto=True, # start automatically
debug = False,
verbose=False):
self.verbose = verbose
self.debug = debug
# LIGAND
dup = MoleculeDuplicator(lig)
self.lig = dup.getCopy()
# parameters(REQUIRED)
self.smartsLig = smartsLig # smarts
self.smartsIdxLig = smartsIdxLig # smarts indices
self.bondIdxLig = bondIdxLig # bonds
self.indicesLig = indicesLig # indices (direct)
# RECEPTOR
self.fullRec = rec
self.resName = resName.upper()
# parameters (OPTIONAL)
self.smartsRec = smartsRec # smarts
self.smartsIdxRec = smartsIdxRec # smarts indices
self.bondIdxRec = bondIdxRec # bonds
self.indicesRec = indicesRec # indices (direct)
# collect receptor information
self.genFullRec = genFullRec
self.setMode()
self.initResidue()
if auto:
self.process()
def vprint(self, msg):
""" """
if not self.verbose:
return
print("VERBOSE " + msg)
def setMode(self):
""" """
# ligand mode
if self.indicesLig:
self._ligmode = 'idx'
elif self.smartsLig:
self._ligmode = 'smarts'
elif self.bondIdxLig:
self._ligmode = 'bond'
self.vprint("[setMode] ligand alignment mode: %s" % self._ligmode)
# receptor mode
if self.indicesRec:
self._recmode = 'idx'
elif self.smartsRec:
self._recmode = 'smarts'
elif self.bondIdxRec:
self._recmode = 'bond'
else:
self._recmode = 'smarts'
self.vprint("[setMode] no mode specified-> switching to
default")
self.vprint("[setMode] receptor alignment mode: %s" % self._recmode)
def initResidue(self):
""" initialize residue information"""
string = self.resName
chain, res = string.split(":")
self.resInfo = { 'chain' : chain,
'name' : res[0:3],
'num': int(res[3:]),
}
if self.debug: print "PROCESSING RESIDUE", self.resName
self.residueProfile()
def process(self):
""" do the stuff"""
if self.verbose: print "Aligning ligand on residue...", self.resName
self.align()
self.cleanup()
def residueProfile(self):
""" tries to recognize the receptor residue
and rename atoms to facilitate Calpha-Cbeta identification
"""
self.receptorProfiler =
ResidueProfiler.AminoAcidProfiler(self.fullRec,
resId=self.resName, setLabels=True,
auto=True, debug=self.debug)
self.resType = self.receptorProfiler.residues[self.resName]['type']
if self.verbose:
print "Selected residue recognized as:", self.resType.upper()
dup = MoleculeDuplicator(self.fullRec, resList=[self.resName])
self.modifiedRes = dup.getCopy()
if self._recmode == 'idx':
org = self.indicesRec[:]
self.indicesRec = [ self.modifiedRes._numbering[x] for x in
self.indicesRec]
self.vprint( ('converting full receptor indices [%s] to '
'residue copy indices [%s]' % (org, self.indicesRec)))
self.setResidueParms()
def setResidueParms(self):
""" populate automatically SMARTS and indices for
the residue if nothing was specified already
"""
if not self._recmode == 'smarts':
msg = "[setResidueParm] receptor alignment mode is (%s),
skipping"
self.vprint (msg % self._recmode)
return
if self.smartsRec:
msg = ("[setResidueParms] receptor SMARTS already defined [%s],
skipping")
self.vprint(msg % self.smartsRec)
return
self._residueParms = {'cys' : { 'patternIdx' : [3,2], # when
using SMARTS patterns from ResProfiler
'atNames': ['SG', 'CB'], # when looking
for [setMode] atom names in the structure
},
'ser' : { 'patternIdx' : [3,2],
'atNames' : ['OG', 'CB'],
},
'thr' : { 'patternIdx' : [4,2],
'atNames' : ['OG', 'CB'],
},
'lys' : { 'patternIdx': [6,5],
'atNames': ['NZ','CE'],
},
}
resInfo = self._residueParms[self.resType]
self.smartsRec = self.receptorProfiler.getResPattern(self.resType)
self.smartsIdxRec = resInfo['patternIdx']
if self.verbose:
print "Auto-parameters for residue (%s)" % (self.resType)
print " SMARTS pattern : %s" % self.smartsRec
print " atom indices : %s" % self.smartsIdxRec
print " atom names : %s" % resInfo['atNames']
def align(self):
""" where the magic happens... perform the alignment"""
self.populateIndices()
if not self.ready:
print "Molecule not ready... Aborting."
return
aligner = VectorMolAligner(lig=self.lig, rec=self.modifiedRes,
ligIndices=self.indicesLig, recIndices=self.indicesRec,
verbose=self.verbose, debug=self.debug)
def populateIndices(self):
""" populate indices pairs accordingly to the modes"""
# XXX add a check for missing matches
# ligand indices:
self.ready = True
if self._ligmode == 'smarts':
smarter = IndicesFromSMARTS(mol = self.lig,
pattern = self.smartsLig, indices = self.smartsIdxLig,
firstOnly = False, debug = self.debug)
if not smarter.matches:
print "*** WARNING *** Pattern not found in the ligand"
self.ready = False
return
self.indicesLig = smarter.matches[0]
elif self._ligmode == 'bond':
pass
else:
self.vprint("[align] Using provided indices for ligand")
# receptor indices
if self._recmode == 'smarts':
smarter = IndicesFromSMARTS(mol = self.modifiedRes,
pattern = self.smartsRec, indices = self.smartsIdxRec,
debug = self.debug)
if not smarter.matches:
print "*** WARNING *** Pattern not found in the receptor"
self.ready = False
return
self.indicesRec = smarter.matches[0]
elif self._recmode == 'bond':
pass
else:
self.vprint("[align] Using provided indices for receptor")
def cleanup(self):
""" Lig: remove from A on...
Rec: remove from I on...
merge from CA on to the residue
"""
if not self.ready:
return
self.cleanLigand()
self.cleanReceptor()
self.cleanModRes()
if self.genFullRec:
self.generateFullRec()
def cleanLigand(self):
""" remove atoms A and B from the ligand"""
# delete ligand atoms used to overlap
for a in [ self.lig.GetAtom(i) for i in self.indicesLig ]:
self.lig.DeleteAtom(a)
# delete hydrogens
self.lig.DeleteHydrogens()
# create new molecule (XXX maybe it should be moved at the
beginning?
# before deleting atoms
dup = MoleculeDuplicator(self.lig, name='Covalent docking of...')
self.covalentLigand = dup.getCopy()
# XXX BUG this triggers the creation of a residue "LIG" for ligands
# with no residues
self.covalentLigand.GetAtom(1).GetResidue()
covalentRes = self.covalentLigand.GetResidue(0)
#if covalentRes:
covalentRes.SetName(self.resType.upper())
# renaming ligand atoms 'C -> CX' to avoid issues later (ADT...)
for atom in ob.OBResidueAtomIter(covalentRes):
oldIdFull = covalentRes.GetAtomID(atom)
# keep the element letter, throw the rest...
oldId = oldIdFull[0]
covalentRes.SetAtomID(atom, oldId+'X')
if self.debug:
pybel.Molecule(self.lig).write('pdb', 'DEBUG_cleanLig.pdb',
overwrite=1)
def cleanReceptor(self):
""" use PDB atom labels to decide atoms to be kept
"""
# XXX CORRECT XXX
#accepted = ['CA', 'N'] +
self.receptorProfiler.getResLabels(self.resType,
# backbone=False)
accepted = ['CA'] + self.receptorProfiler.getResLabels(self.resType,
backbone=True)
if self.debug:
print "Accepted residue atom names:" , accepted
covalentRes = self.covalentLigand.GetResidue(0)
for atom in ob.OBMolAtomIter(self.modifiedRes):
aName = self.getAtomName(atom, True)
if self.debug: print "Scanning atom [%s]" % aName
if aName in accepted:
if self.debug: print "[ accepted ]"
newName = ' {0:3}'.format(aName)
newAtom = self.covalentLigand.NewAtom()
newAtom.Duplicate(atom)
covalentRes.AddAtom(newAtom)
# residue atoms will be ATOM
covalentRes.SetAtomID(newAtom, newName)
if self.debug:
mol = pybel.Molecule(self.covalentLigand)
mol.write('pdb', 'final.pdb', overwrite=1)
def cleanModRes(self, setHet=True, connect=True, addH=True):
""" perform last clean-up operations:
conect : autodetect bonds
addH : add hydrogens
setHet : set all atom records as HETATM
"""
# autodetect bonds
if connect:
self.covalentLigand.ConnectTheDots()
# add hydrogens
if addH:
self.addHydrogens()
# set all atoms to be HETATM records
covalentRes = [ x for x in ob.OBResidueIter(self.covalentLigand)][0]
if setHet:
for a in ob.OBMolAtomIter(self.covalentLigand):
covalentRes.SetHetAtom(a, True)
covalentRes.SetName(self.resType.upper())
chain, res = self.resName.split(':',1)
print "CHAIN", chain
rnum = res[3:]
covalentRes.SetNum(int(rnum))
covalentRes.SetChain(chain)
def generateFullRec(self):
""" replace the original residue with the
modified one
"""
rName = self.resInfo['name']
rChain = self.resInfo['chain']
rNum = self.resInfo['num']
for res in ob.OBResidueIter(self.rec):
if (res.GetChain() == rChain):
if (res.GetNum() == rNum):
self.rec.DeleteResidue(res)
if self.verbose:
print "deleted residue [%s:%s%s]" % (rChain,
rName, rNum)
break
# add modified residue
#covalentRes = [ x for x in
ob.OBResidueIter(self.covalentLigand)][0]
#self.rec.AddResidue(covalentRes)
def addHydrogens(self):
""" manages all the gymnastics for adding hydrogens at pH 7.4
(DEFAULT)"""
#self._temporaryPHFix()
self.covalentLigand.DeleteHydrogens()
self.covalentLigand.UnsetFlag(ob.OB_PH_CORRECTED_MOL)
for a in ob.OBMolAtomIter(self.covalentLigand):
a.SetFormalCharge(0)
self.covalentLigand.SetAutomaticFormalCharge(True)
#self.covalentLigand.PerceiveBondOrders()
self.covalentLigand.AddHydrogens(False, True )#, pH)
#self.covalentLigand.AddHydrogens(False,True)
def _temporaryPHFix(self):
""" this is a nasty workaround for an even nastier bug in OB:
http://sourceforge.net/p/openbabel/bugs/710/
"""
t = 'mol2'
n = '.tmpmol'
m = pybel.Molecule(self.covalentLigand)
m.write(t, n, overwrite=1)
self.covalentLigand = pybel.readfile(t, n).next().OBMol
os.remove(n)
def getAtomName(self, atom, strip=False):
""" return the PDB atom name """
res = atom.GetResidue()
if strip:
return res.GetAtomID(atom).strip()
return res.GetAtomID(atom)
def writeCovalent(self, outFname, _format='pdb'):
""" _format could be any OB format
NOTE: test mol2 to avoid spurious bonds?
"""
if self.genFullRec:
mol = self.rec
else:
mol = self.covalentLigand
if _format == 'pdbqt':
self._writePdbqt(outFname+'.pdbqt')
else:
#fname = "%s.%s" % ( outFname, _format.lower() )
fname = outFname
conv = ob.OBConversion()
conv.SetOutFormat(_format)
conv.WriteFile(mol, fname)
if self.verbose: print "Written: %s" % (fname)
def _writePdbqt(self, outFname):
""" not necessary, for now, we will use prepare_flexres code
later"""
if self.genFullRec:
mol = self.rec
else:
mol = self.covalentLigand
class IndicesFromSMARTS:
""" Given a SMARTS pattern, return the atom indices matching
the specified atoms in the pattern
"""
def __init__(self, mol, pattern, indices=(), firstOnly = True,
debug=False):
""" mol : OBMol()
pattern = SMARTS string
[ indices: indices of the SMARTS atoms that need to be
identified ]
[ firstOnly: by default only first match is returned; multiple
are possible]
"""
self.debug = debug
self.mol = mol
if self.debug:
print "IndicesFromSMARTS> MOLECULE NAME", self.mol.GetTitle()
print "IndicesFromSMARTS> MOLECULE ", self.mol
print "PATTERN", pattern
self.pattern = pattern
self.indices = indices
if self.indices == ():
self.indices = (0,1)
print "pointMatch> WARNING! missing indices, using default
(0,1)"
self.matches = []
self.pointMatch()
def pointMatch(self):
""" correct the indexing from the SMARTS pattern indexing to
the 0-based used to retrieve atoms
"""
result = self.findSmarts()
if not len(result):
print "IndexFromSMARTS: NO MATCH FOUND! troubles ahead..."
return
if self.debug:
print "[IndicesFromSMARTS] >>> raw out", result
if len(result) > 1:
print "TMP: MULTIPLE MATCHES FOUND!", len(result)
i,j = self.indices
if self.debug:
print "[IndicesFromSMARTS]> extracting indices [%d,%d] from
smarts[%s]" % (i,j,result)
self.matches = [ (r[i],r[j]) for r in result ]
if self.debug:
print "FOUND!", self.matches
for a in ob.OBMolAtomIter(self.mol):
if a.GetIdx() in self.matches[0]:
res = a.GetResidue()
print "RESNAME", res.GetName()
return self.matches
def findSmarts(self):
""" OB SMARTS matcher """
obpat = ob.OBSmartsPattern()
obpat.Init(self.pattern)
obpat.Match(self.mol)
if self.debug:
print "PROCESSING MOLECULE", self.getSmi()
return [ x for x in obpat.GetUMapList() ]
def getSmi(self):
""" """
m = pybel.Molecule(self.mol)
return m.__str__()
# XXX INACCURATE! PROBABLY TO BE DROPPED?
class IndicesFromBond:
""" identify atom indices involved in a bond"""
def __init__(self, mol, bond):
""" The requirement
XXX the requirement is that there must be a unique atom type that's
common? -X-S S-Y-
Some smartness could be used to guess what is the 'terminal'
direction to infer the proper orientation of the molecules
"""
pass
class VectorMolAligner:
""" """
def __init__(self, lig, rec, ligIndices=(), recIndices=(), auto=True,
verbose=False, debug=False):
"""
Given two OB Molecules and two pairs of atom indices
(A,B) and (I,J) perform the covalent alignment
0. Initial configuration
B
/
[L]----A I---J--[R]
/
X
1. Translate [L] to overlap A->I
B
/
[L]----A -> I---J--[R]
/
X
B
/
[L]----AI---J--[R]
/
X
2. Rotate [L] around the normal of plane (B,A,J) to overlap B->J
[L]
\
AI--BJ--[R]
/
X
Indexing provided must be as following:
ligIndices = (B, A)
recIndices = (I, J)
"""
self.debug = debug
self.verbose = verbose
self.lig = lig
self.rec = rec
self.ligIndices = ligIndices
self.recIndices = recIndices
if auto:
self.align()
def align(self):
""" perform the alignment based on vectors
"""
self.getVectors()
self.translate()
self.rotate()
def getVectors(self):
""" convert all lig and rec indices to vectors
"""
# XXX ADAPT THIS TO WORK WITH MULTIPLE MATCHES
if self.debug: print "GETTING LIGAND VECTORS..."
self._ligB = self._idxToVec(self.lig, self.ligIndices[0])
self._ligA = self._idxToVec(self.lig, self.ligIndices[1])
if self.debug: print "GETTING RECEPTOR VECTORS..."
self._recI = self._idxToVec(self.rec, self.recIndices[0],
translate=0) # , reverse=)
self._recJ = self._idxToVec(self.rec, self.recIndices[1],
translate=0) # , reverse=True)
#self._recI = self._idxToVec(self.rec, self.recIndices[0],
translate=True, reverse=True)
#self._recJ = self._idxToVec(self.rec, self.recIndices[1],
translate=True, reverse=True)
#self._recI = self._idxToVec(self.rec, self.recIndices[0],
translate=True) # , reverse=)
#self._recJ = self._idxToVec(self.rec, self.recIndices[1],
translate=True) # , reverse=True)
if self.debug:
c = 1
buff = ['MODEL']
buff.append( makePdb(coord=self._ligB, at_index=c, atype='P'))
c+=1
buff.append( makePdb(coord=self._ligA, at_index=c, atype='P'))
c+=1
buff.append( makePdb(coord=self._recI, at_index=c, atype='S'))
c+=1
buff.append( makePdb(coord=self._recJ, at_index=c, atype='S'))
buff.append('ENDMDL')
writeList('DEBUG_getVectors.pdb', buff, addNewLine=1, mode='a')
def _idxToVec(self, mol, atIdx, convert=True,translate=False):
""" return the vector of the atom at index atIdx
if convert requested, return Numpy array instead of OB.Vector3
"""
#print "CALLED BY", mol, atIdx, convert, translate
atIdx = int(atIdx)
if translate:
# DuplicateMolecule obj
atom = mol.GetAtom(atIdx, translate=True)
else:
# OBMol object
atom = mol.GetAtom(atIdx)
try:
vec = atom.GetVector()
except:
print "*** ERROR:", sys.exc_info()[1]
print "*** REQUESTED", atIdx
print "*** FOUND",
for a in ob.OBMolAtomIter(mol):
print a.GetIdx(),
print " "
return
if not convert:
return vec
return self.obVecToNumpy( vec )
def translate(self):
""" perform the 1. translation step """
if self.verbose: print "Translation...",
translation = vector( self._ligA, self._recI)
self._translateMol(self.lig, translation)
# update vector coords after translation
self.getVectors()
if self.debug:
pybel.Molecule(self.lig).write('pdb',
'DEBUG_translated.pdb', overwrite=True)
if self.verbose: print "[ DONE ]"
def rotate(self):
""" perform the 2. rotation step"""
import math
if self.verbose: print "Rotation...",
## rotation
p1 = self._ligB
p2 = self._ligA
p3 = self._recJ
p4 = self._recI
v1 = vector(p1,p2)
v2 = vector(p3,p2)
angle = vecAngle(v1, v2) # rad
center = vector(p2)
normVec = calcPlane( p1, p2, p3)
if angle < 0.09:
if self.verbose:
print ( '[ skipping rotation, angle near-zero (%1.3f ), '
'ligand possibly already in place' % (angle) )
#return
pass
if self.debug:
print "Rotation angle:", math.degrees(angle)
pdb1 = makePdb(coord=p1, at_index = 1, atype ='S')
pdb2 = makePdb(coord=p2, at_index = 2, atype ='S')
pdb3 = makePdb(coord=p3, at_index = 3, atype ='S')
pdb4 = makePdb(coord=normVec+p2, at_index = 4, atype ='N')
self._testrot(p1, p2, normVec) # writes rotor.pdb
writeList('DEBUG_plane.pdb', [pdb1, pdb2, pdb3, pdb4],
addNewLine=1)
rotor = ( normVec[0], normVec[1], normVec[2], angle )
self._rotateMol( self.lig, p2, rotor)
if self.debug:
pybel.Molecule(self.lig).write('pdb', 'DEBUG_rotated.pdb',
overwrite=True)
if self.verbose: print "[ DONE ]"
return
def _translateMol(self, mol, v):
""" translate the molecule using vector3"""
for i in range(mol.NumAtoms()):
atom = mol.GetAtom(i+1)
pre = (atom.GetX(), atom.GetY(), atom.GetZ())
post = (pre[0]+v[0], pre[1]+v[1], pre[2]+v[2])
obvec = ob.vector3()
obvec.Set(post[0], post[1], post[2])
atom.SetVector(obvec)
return
def _rotateMol(self, mol, center, rotor):
""" rotate molecule atoms applying vector around center"""
for a in ob.OBMolAtomIter(mol):
coord = self.obVecToNumpy( a.GetVector() )
vec = vector(center, coord)
newCoord = rotatePoint( vec, center, rotor)
x,y,z = self.vecToDouble(newCoord)
vec = ob.vector3()
vec.Set(x,y,z)
a.SetVector(vec)
def _testrot(self, point, center, norm):
""" debug function to depict the rotor applied in the mol
rotation"""
buff = []
buff = [ makePdb((0.,0.,0.), at_index =1, atype='S') ]
buff.append( makePdb((1.8,0.,0.), at_index =2, atype='S') )
buff.append( makePdb((0.,1.8,0.), at_index =3, atype='S') )
buff.append( makePdb((0.,0.,1.8), at_index =4, atype='S') )
buff.append( makePdb(point, at_index =5, atype='S') )
c = 6
for a in range(0, 360, 30):
a = math.radians(a)
rotor = ( norm[0], norm[1], norm[2], a)
coord = vector(point,center)
new = rotatePoint( coord, center, rotor)
atom = makePdb(new, at_index =c)
buff.append(atom)
writeList('DEBUG_rotor.pdb', buff, addNewLine=1)
def obVecToNumpy(self, ObVec):
""" """
return np.array( ( ObVec.GetX(), ObVec.GetY(), ObVec.GetZ() ), 'f')
def vecToDouble(self, vector):
""" """
return float(vector[0]), float(vector[1]), float(vector[2])
class CovalentReaction:
def __init__(self):
""" simmeglio """
""" IT SHOULD be used to attach the next group to the reaction
to be used for the alignment """
pass
# ########################
class CovalentDockingMaster:
""" manage settings for the CovalentDocking object"""
def __init__(self, debug=False, verbose=False):
""" """
self.debug = debug
self.verbose = verbose
self.initOpts()
self.parseOpts()
self.start()
def getFnameExt(self,string):
""" """
name, ext = os.path.splitext(string)
ext = ext[1:]
return name, ext
def initOpts(self):
""" """
self.initDefaults()
self.initDocs()
self.opts = {
'--ligand' : {'help': 'ligand file (OPTIONAL)',
'action':'store', #'required':True,
'type':str},
'--receptor' : {'help':'receptor file (REQUIRED)',
'action':'store', 'required':True, 'type':str},
'--residue' : {'help': ('residue to modify; the format is
"chain:resNUM (es. "B:THR276")'
'; multiple residues can be specified by
repeating --residue (note: '
'case insensitive)'), 'action':'append',
'default':[], 'type':str},
### '--residuelist' : {'help' : 'modify residues read from file (one
per line)',
### 'action':'store', 'metavar': 'FILE',
'default':None},
###
### '--covresidue' : {'help':'prepare modified residue for the
covalent residue mode (default)',
### 'default': True, 'action': 'store_true'},
###
### '--covmap' : {'help': 'calculate coordinates for Z-map covalent
method',
### 'default': False, 'action': 'store_true'},
###
###
'--ligsmarts' : {'help':('ligand SMARTS pattern; by default, the
first two atoms of '
'the pattern are used; use --ligindices to specify
different atoms'),
'action':'store', 'type':str, 'metavar':'SMARTS',
'default':None},
'--ligindices' : {'help':('indices of ligand atoms to use for
alignment (i.e X-Y-lig)'
'starting from 1; by default, indices refer to atoms in the
ligand file; if --ligsmarts '
'is used, indices specify atoms in the SMARTS pattern'),
'action':'store', 'type':str, 'metavar':'X,Y'},
'--recsmarts' : {'help':('receptor SMARTS pattern; by default, the
first two atoms of '
'the pattern are used; use --recindices to specify
different atoms'),
'action':'store', 'type':str, 'metavar':'SMARTS',
'default':None},
'--recindices' : {'help':('indices of receptor atoms to use for
alignment (i.e. I-J-rec) '
'starting from 1; by default, indices refer to atoms in the
ligand file; if --recsmarts '
'is used, indices specify atoms in the SMARTS pattern'),
'action':'store', 'type':str, 'metavar':'I,J'},
'--genfullrec' : {'help':('generate combined receptor-ligand
structure; '
'by default, only the modified residue (including the
ligand) is written '
'(default : %s)' % self.default_genfullrec ) ,
'action':'store_true' },
'--outputfile' : {'help':('save the output in the file specified'
'(default : use input file'),
'action':'store', 'type':str, 'metavar':'FILENAME'},
### '--generatelig' : {'help' : ('generate structure of flexible
covalent ligand for AutoDock'),
### 'action':'store_true',
'default':self.default_flex},
###
'--generaterec' : {'help' : ('generate structure')},
'--verbose':{'help': ('enable verbose mode'), 'action':'store_true',
'default':False},
'--log' : {'help' : ('write log file'), 'action':'store', 'metavar'
: 'LOGFILE'},
'--help_advanced':{'help': ('show advanced help and documentation'),
'action': 'store_true', 'default':False}
}
self.groups_order = ['INPUT STRUCTURES', 'RESIDUE SPECIFICATION',#
'MODE',
'LIGAND ALIGNMENT DATA', 'RECEPTOR ALIGNMENT DATA (OPTIONAL)',
'OUTPUT', 'LOGGING/HELP']
self.groups = { 'INPUT STRUCTURES' : ['--receptor', '--ligand'],
'RESIDUE SPECIFICATION' : ['--residue' ],
#'MODE' : ['--covresidue', '--covmap'],
'LIGAND ALIGNMENT DATA': ['--ligindices',
'--ligsmarts'],
'RECEPTOR ALIGNMENT DATA (OPTIONAL)':
['--recindices', '--recsmarts'],
'OUTPUT' : [ '--generaterec', '--genfullrec',
'--outputfile'],
'LOGGING/HELP':['--log', '--verbose',
'--help_advanced'],
}
self.usage = '%s --ligand xxxx.xxx --receptor xxxx.xxx [ alignment
mode ]'
self.usage = None
def initDocs(self):
""" """
self.description = (
"""Raccoon PrepareCovalent: typical usages
Create input ligand/receptor structures for covalent docking:
%(name)s --receptor receptor.pdb --ligand ligand.pdb --residue A:CYS113
--ligsmarts CSCC --
Create input ligand/receptor structures for covalent docking:
%(name)s --receptor receptor.pdb --ligand ligand.pdb --residue A:CYS113
--ligsmarts CSCC
"""% { 'name': sys.argv[0] } )
self.epilog = "(C) 2014 Stefano Forli, MGL, TSRI\n Please cite..."
def initDefaults(self):
""" """
self.default_ligPatternIdx = (0,1)
self.default_recPatternIdx = None #(0,1)
self.default_genfullrec = False
self.default_ligpatternidx = (0,1)
self.default_recpatternidx = None #(0,1)
self.default_outputfile = None
self.default_flex = False
def showAdvancedHelp(self):
""" """
advHelp = ('[ADVANCED HELP GOES HERE')
print advHelp
sys.exit(0)
def parseForeignOpts(self):
""" remove foreign options from the command line"""
foreign = {'prepare': []}
# debug
if '--debug' in sys.argv[1:]:
print "\n**** Debug mode activated ***\n\n"
self.debug = True
sys.argv.remove('--debug')
# prepare
for k in foreign.keys():
remove = []
for opt in sys.argv[1:]:
if opt[2:].startswith(k):
foreign[k] = opt
remove.append(opt)
for opt in remove:
sys.argv.remove(opt)
def parseOpts(self):
""" parse command-line options"""
self.parseForeignOpts()
self.parser = argparse.ArgumentParser(description =
self.description,
usage=self.usage, epilog=self.epilog,
formatter_class = argparse.RawDescriptionHelpFormatter)
for g in self.groups_order:
group = self.parser.add_argument_group(g)
for i in self.groups[g]:
opts = self.opts[i]
group.add_argument(i, **opts)
self.args = self.parser.parse_args()
if self.args.help_advanced:
self.showAdvancedHelp()
self.alignerArgs = {}
self.verbose = self.alignerArgs['verbose'] = self.args.verbose
self.validateOpts()
if hasattr(self.args, 'outputfile'):
self._output = self.args.outputfile
else:
self._output = None
def validateOpts(self):
""" validate all options and set the modes for the actual
calculation"""
self.validateMode()
self.setLigandMode()
self.setReceptorMode()
def setLigandMode(self):
""" guess the ligand mode basing on the spefication of indices or
SMARTS
--ligindices
--ligsmarts [ --ligindices ]
"""
args = self.alignerArgs #= {}
self._ligmode = 'idx'
idxError = ('Ligand indices (--ligindices) must be specified as
"number,number" (starting from 1)')
# ligand mode
if self.args.ligsmarts:
self._ligmode = 'smarts'
args['smartsLig'] = self.args.ligsmarts
if self.args.ligindices:
try:
# smarts indices are 0-based internally
args['smartsIdxLig'] = [int(x)-1 for x in
self.args.ligindices.split(',')]
except:
self.showError(idxError)
else:
args['smartsIdxLig'] = self.default_ligPatternIdx
else:
if not self.args.ligindices:
msg = ('Ligand indices (--ligindices) are required, or use
SMARTS (--ligsmarts)')
self.showError(msg)
else:
try:
args['indicesLig'] = [int(x) for x in
self.args.ligindices.split(',')]
except:
self.showError(idxError)
self.vprint("LIGAND MODE:"+self._ligmode)
def setReceptorMode(self):
""" set receptor alignment mode (OPTIONAL)"""
args = self.alignerArgs #= {}
self._recmode = None
idxError = ('Receptor indices (--recindices) must be specified as
"number,number" (starting from 1)')
# receptor mode
if self.args.recsmarts:
self._recmode = 'smarts'
args['smartsRec'] = self.args.recsmarts
if self.args.recindices:
try:
# smarts indices are 0-based internally
args['smartsIdxRec'] = [int(x)-1 for x in
self.args.recindices.split(',')]
except:
self.showError(idxError)
else:
args['smartsIdxRec'] = self.default_recPatternIdx
else:
if self.args.recindices:
self._recmode = 'idx'
try:
args['indicesRec'] = [int(x) for x in
self.args.recindices.split(',')]
except:
self.showError(idxError)
# check that specified residue is managed, then no SMARTS/indices is
required
if (self._recmode == None):
unsupported = self.checkManagedResidue()
if len(unsupported):
msg = ('The specified residue types are not supported
automatically: %s. Use either '
'--recsmarts or --recindices.' %
(','.join(unsupported)))
self.showError(msg)
args['genFullRec'] = self.args.genfullrec
def validateMode(self):
""" set session mode to be:
- normal (lig + receptor + residue|residuetype)
"""
self.mode = 'normal'
if not self.args.residue:
msg = ('No residue specification. Use at least --residue to
specify which residue to modify')
self.showError(msg)
if self.args.residue and not self.args.ligand:
msg = ('option --residue requires --ligand to be specified')
self.showError(msg)
def checkManagedResidue(self):
""" check that requested residue type is
among managed residues: CYS, SER, THR, LYS
"""
supported = ['cys', 'ser', 'thr', 'lys']
unsupported = []
for r in self.args.residue:
if ':' in r:
r = r.split(":")[1][:3].lower()
if not r in supported:
unsupported.append(r)
return set([ x.upper() for x in unsupported])
def parseLigAlignmentOpts(self):
""" validate parameters for ligand alignment """
self._ligmode = None
# SMART INDICES
self.ligsmarts = self.args.ligsmarts
if self.ligsmarts:
self._ligmode = 'smarts'
self.ligsmartsidx = self.args.ligsmartsidx
try:
if not self.ligsmartsidx == None:
self.ligsmartsidx = [ int(x) for x in
self.ligsmartsidx.split(',') ]
except:
print "ERROR (--ligsmartsidx): indices must be specified with
number pairs, e.g. 1,3"
sys.exit(1)
self._ligmode = 'smarts'
# ATOM INDICES
self.ligatomidx = self.args.ligatomidx
if self.ligatomidx and self._ligmode == 'smarts':
msg = ('Multiple alignment modes selected: use either SMARTS
(\'--ligsmarts\')'
'or atom indices(\'--ligatomidx\')')
self.showError(msg)
# INDICES
try:
if not self.ligatomidx == None:
self.ligatomidx = [ int(x) for x in
self.ligatomidx.split(',') ]
except:
print "ERROR (--ligatomidx): indices must be specified with
number pairs, e.g. 1,3"
sys.exit(1)
def parseRecAlignmentOpts(self):
""" """
# validate receptor SMARTS indices
try:
self.recsmartsidx = self.args.recsmartsidx
if not self.recsmartsidx == None:
self.recsmartsidx = [ int(x) for x in
self.recsmartsidx.split(',') ]
except:
print "ERROR (--recsmartsidx): indices must be specified with
number pairs, e.g. 1,3"
sys.exit(1)
# validate receptor atom indices
try:
self.recatomidx = self.args.recatomidx
if not self.recatomidx == None:
self.recatomidx = [ int(x) for x in
self.recatomidx.split(',') ]
except:
print "ERROR (--recatomidx): indices must be specified with
number pairs, e.g. 1,3"
sys.exit(1)
def showError(self, msg, showhelp=True):
""" error reporting facility (modify to STDERR)"""
tag = '[ ERROR ] '
msg = msg.replace('\n', ('\n%s' % (" " * len(tag) ) ) )
print("%s%s" % (tag, msg))
if showhelp:
print ""
self.parser.print_help()
sys.exit(1)
def vprint(self, msg):
""" verbose printer"""
if not self.verbose:
return
print("VERBOSE: "+msg)
def getMolecule(self, fname, perceiveChains=False):
""" molecule parser"""
name, ext = self.getFnameExt(fname)
self.molParser.SetInFormat(ext)
mol = ob.OBMol()
result = False
try:
result = self.molParser.ReadFile(mol, fname)
except:
msg = ('Fail to read molecule "%s" : %s' % (fname,
sys.exc_info()[1]))
self.showError(msg, showhelp=False)
if not result:
msg = ('Fail to read molecule "%s"' % (fname))
self.showError(msg, showhelp=False)
if perceiveChains:
self.chainParser.PerceiveChains(mol)
return mol, name
def loadMolecules(self):
""" import ligand and receptor structures"""
args = self.alignerArgs
args['rec'], self.recName = self.getMolecule(self.args.receptor)
args['lig'], self.ligName = self.getMolecule(self.args.ligand, True)
def getResidues(self, log=None):
""" find which residue(s) will be modified:
- specific residues
- residue types (surface-accessible residues)
"""
residue = self.args.residue
if not residue:
msg = ('No residues specified. Use either to process!')
self.showError(msg)
self.residuePool = residue
chain = None
if len(self.residuePool) == 0:
msg = ('No residues to process!')
self.showError(msg)
self.residuePool.sort()
if log:
fp = open(log, 'w')
for r in self.residuePool:
#fp.write( '%s:%s\n' % (self.recName,r))
fp.write( '%s\n' % (r))
fp.close()
#print "========================================================"
#print "Total number of residues to process [%d]" %
len(self.residuePool)
def initParsers(self):
""" initialize molecule parser and chain perception"""
self.molParser = ob.OBConversion()
self.chainParser = ob.OBChainsParser()
def start(self):
""" """
# XXX validate here that the indices are within the lenght of
matches?
self.initParsers()
self.loadMolecules()
self.getResidues('residues.log')
self.processResidues()
def processResidues(self):
""" perform full covalent protocol"""
outType = 'pdb'
outType = 'mol2'
for count, residue in enumerate(self.residuePool):
#print "Processing residue %s\t (%d/%d)" % (residue,
# count+1, len(self.residuePool)),
print "Processing residue %s" % residue
if not self._output == None:
outname = self._output
outType = os.path.splitext(self._output)[1][1:]
else:
outname = '%s_%s_cov_%s' % (os.path.basename(self.ligName),
os.path.basename(self.recName),
residue.replace(':','').upper())
outname = outname + "." + outType
self.vprint("[start] output filename is: "+ outname)
print ("[start] output filename is: "+ outname)
self.alignerArgs['resName'] = residue
if self.verbose:
print "==============="
print "[start] Calling CovalentDockingMaker with settings:"
for k,v in self.alignerArgs.items():
if isinstance(v, ob.OBMol):
print "KW[%s]: %d atoms" % (k, v.NumAtoms())
else:
print "KW[%s]:" % k, v
print "==============="
#print "\n==============\n", self.alignerArgs
if self.debug:
self.alignerArgs['debug'] = True
aligner = self.x = CovalentDockingMaker(**self.alignerArgs)
if not aligner.ready:
print "\t\t\tSKIPPING?"
continue
#print "\t=>", outname + '.pdb'
# XXX LOCAL FILE HERE
fname = os.path.basename(outname) # + '.pdb'
#aligner.writeCovalent(outname, 'pdb')
aligner.writeCovalent(fname, outType)
print "Writing output filename: %s" % fname
#aligner.writeCovalent(outname, 'pdbqt')
# WORKING PATTERNS
# LIG REC
# -------------------------------
# CSCC '[SH]C'
# '[$(SH1);$(SC)]'
# SCC '[$([SH1]);$(SC)]'
if __name__ == '__main__':
x = CovalentDockingMaster()
--
View this message in context:
http://forums.openbabel.org/Openbabel-error-in-ReadFile-tp4660198.html
Sent from the General discussion mailing list archive at Nabble.com.
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
OpenBabel-discuss mailing list
OpenBabel-discuss@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/openbabel-discuss
