I download the adcovalent package to attempt covalent docking with AutoDock Vina, I have Anaconda, and openbabel installed. pybel.py is a wrapper for openbabel, it looks to me like pybel.py is broken, although I have the most recent version.
This is the script I'm trying to run:
#!/usr/bin/env python
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
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
bfactor = "%2.2f" % bfactor
if len(bfactor.split(".")[0]) == 1:
bfactor = " "+bfactor
if len(atype) == 1:
atype = atype + " "
"%-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 which 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:
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)
self.fullRec = rec
self.resName = resName.upper()
self.smartsRec = smartsRec # smarts
self.smartsIdxRec = smartsIdxRec # smarts indices
self.bondIdxRec = bondIdxRec # bonds
self.indicesRec = indicesRec # indices (direct)
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):
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)
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):
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):
if self.verbose: print "Aligning ligand on residue...", self.resName
self.align()
self.cleanup()
def residueProfile(self):
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)))
###BLAH, BLAH... it goes on...
This is the error I am getting:
(vina) gorrie06@DMGLaptop:~/docking/covalent$ ptor /3upo_test/3upo_protein.pdb --residue B:SER222 --outputfile ligcov.pdb
raceback (most recent call last):
File "/home/gorrie06/docking/covalent/adcovalent/prepareCovalent.py", line 36, in <module>
import pybel
File "/home/gorrie06/anaconda3/envs/vina/lib/python2.7/site-packages/pybel.py", line 89, in <module>
informats = _formatstodict(_obconv.GetSupportedInputFormat())
File "/home/gorrie06/anaconda3/envs/vina/lib/python2.7/site-packages/pybel.py", line 68, in _formatstodict
broken = [(x, y.strip()) for x, y in broken]
ValueError: need more than 1 value to unpack
This is the function in pybel that is broken:
def _formatstodict(list):
if sys.platform[:4] == "java":
list = [list.get(i) for i in range(list.size())]
broken = [x.replace("[Read-only]", "").replace("[Write-only]", "").split(
" -- ") for x in list]
broken = [(x, y.strip()) for x, y in broken]
return dict(broken)
It comes from an openbabel wrapper, pybel:
# -*. coding: utf-8 -*-
# Copyright (c) 2008-2012, Noel O'Boyle; 2012, Adrià Cereto-Massagué
# All rights reserved.
#
# This file is part of Cinfony.
# The contents are covered by the terms of the GPL v2 license
# which is included in the file LICENSE_GPLv2.txt.
"""
pybel - A Cinfony module for accessing Open Babel
Global variables:
ob - the underlying SWIG bindings for Open Babel
informats - a dictionary of supported input formats
outformats - a dictionary of supported output formats
descs - a list of supported descriptors
fps - a list of supported fingerprint types
forcefields - a list of supported forcefields
"""
import sys
import os.path
import tempfile
import json
import uuid
import xml.etree.ElementTree as ET
if sys.platform[:4] == "java":
import org.openbabel as ob
import java.lang.System
java.lang.System.loadLibrary("openbabel_java")
_obfuncs = ob.openbabel_java
_obconsts = ob.openbabel_javaConstants
import javax
elif sys.platform[:3] == "cli":
import System
import clr
clr.AddReference('System.Windows.Forms')
clr.AddReference('System.Drawing')
from System.Windows.Forms import Application, DockStyle, Form, PictureBox
from System.Windows.Forms import PictureBoxSizeMode
from System.Drawing import Image, Size
_obdotnet = os.environ["OBDOTNET"]
if _obdotnet[0] == '"': # Remove trailing quotes
_obdotnet = _obdotnet[1:-1]
clr.AddReferenceToFileAndPath(os.path.join(_obdotnet, "OBDotNet.dll"))
import OpenBabel as ob
_obfuncs = ob.openbabel_csharp
_obconsts = ob.openbabel_csharp
else:
import openbabel as ob
_obfuncs = _obconsts = ob
try:
import Tkinter as tk
import Image as PIL
import ImageTk as piltk
except ImportError: # pragma: no cover
tk = None
def _formatstodict(list):
if sys.platform[:4] == "java":
list = [list.get(i) for i in range(list.size())]
broken = [x.replace("[Read-only]", "").replace("[Write-only]", "").split(
" -- ") for x in list]
broken = [(x, y.strip()) for x, y in broken]
return dict(broken)
def _getplugins(findplugin, names):
return dict([(x, findplugin(x)) for x in names if findplugin(x)])
def _getpluginnames(ptype):
if sys.platform[:4] == "cli":
plugins = ob.VectorString()
else:
plugins = ob.vectorString()
ob.OBPlugin.ListAsVector(ptype, None, plugins)
if sys.platform[:4] == "java":
plugins = [plugins.get(i) for i in range(plugins.size())]
return [x.split()[0] for x in plugins]
_obconv = ob.OBConversion()
_builder = ob.OBBuilder()
informats = _formatstodict(_obconv.GetSupportedInputFormat())
"""A dictionary of supported input formats"""
outformats = _formatstodict(_obconv.GetSupportedOutputFormat())
"""A dictionary of supported output formats"""
descs = _getpluginnames("descriptors")
"""A list of supported descriptors"""
_descdict = _getplugins(ob.OBDescriptor.FindType, descs)
fps = [_x.lower() for _x in _getpluginnames("fingerprints")]
"""A list of supported fingerprint types"""
_fingerprinters = _getplugins(ob.OBFingerprint.FindFingerprint, fps)
forcefields = [_x.lower() for _x in _getpluginnames("forcefields")]
"""A list of supported forcefields"""
_forcefields = _getplugins(ob.OBForceField.FindType, forcefields)
charges = [_x.lower() for _x in _getpluginnames("charges")]
"""A list of supported charge models"""
_charges = _getplugins(ob.OBChargeModel.FindType, charges)
operations = _getpluginnames("ops")
"""A list of supported operations"""
_operations = _getplugins(ob.OBOp.FindType, operations)
ipython_3d = False
"""Toggles 2D vs 3D molecule representations in IPython notebook"""
def readfile(format, filename, opt=None):
"""Iterate over the molecules in a file.
Required parameters:
format - see the informats variable for a list of available
input formats
filename
Optional parameters:
opt - a dictionary of format-specific options
For format options with no parameters, specify the
value as None.
You can access the first molecule in a file using the next() method
of the iterator (or the next() keyword in Python 3):
mol = readfile("smi", "myfile.smi").next() # Python 2
mol = next(readfile("smi", "myfile.smi")) # Python 3
You can make a list of the molecules in a file using:
mols = list(readfile("smi", "myfile.smi"))
You can iterate over the molecules in a file as shown in the
following code snippet:
>>> atomtotal = 0
>>> for mol in readfile("sdf", "head.sdf"):
... atomtotal += len(mol.atoms)
...
>>> print atomtotal
43
"""
if opt is None:
opt = {}
obconversion = ob.OBConversion()
formatok = obconversion.SetInFormat(format)
for k, v in opt.items():
if v is None:
obconversion.AddOption(k, obconversion.INOPTIONS)
else:
obconversion.AddOption(k, obconversion.INOPTIONS, str(v))
if not formatok:
raise ValueError("%s is not a recognised Open Babel format" % format)
if not os.path.isfile(filename):
raise IOError("No such file: '%s'" % filename)
def filereader():
obmol = ob.OBMol()
notatend = obconversion.ReadFile(obmol, filename)
while notatend:
yield Molecule(obmol)
obmol = ob.OBMol()
notatend = obconversion.Read(obmol)
return filereader()
def readstring(format, string, opt=None):
"""Read in a molecule from a string.
Required parameters:
format - see the informats variable for a list of available
input formats
string
Optional parameters:
opt - a dictionary of format-specific options
For format options with no parameters, specify the
value as None.
Example:
>>> input = "C1=CC=CS1"
>>> mymol = readstring("smi", input)
>>> len(mymol.atoms)
5
"""
if opt is None:
opt = {}
obmol = ob.OBMol()
obconversion = ob.OBConversion()
###it goes on from here...
Any help is greatly appreciated, Thank you!