Skip to content

Instantly share code, notes, and snippets.

@baoilleach
Created May 16, 2011 13:53
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 5 You must be signed in to fork a gist
  • Save baoilleach/974477 to your computer and use it in GitHub Desktop.
Save baoilleach/974477 to your computer and use it in GitHub Desktop.
Using Open Babel to calculate the symmetry-corrected RMSD of a docked pose from a crystal structure
import math
import pybel
def squared_distance(coordsA, coordsB):
"""Find the squared distance between two 3-tuples"""
sqrdist = sum( (a-b)**2 for a, b in zip(coordsA, coordsB) )
return sqrdist
def rmsd(allcoordsA, allcoordsB):
"""Find the RMSD between two lists of 3-tuples"""
deviation = sum(squared_distance(atomA, atomB) for
(atomA, atomB) in zip(allcoordsA, allcoordsB))
return math.sqrt(deviation / float(len(allcoordsA)))
if __name__ == "__main__":
# Read crystal pose
crystal = next(pybel.readfile("mol", "crystalpose.mol"))
# Find automorphisms involving only non-H atoms
mappings = pybel.ob.vvpairUIntUInt()
bitvec = pybel.ob.OBBitVec()
lookup = []
for i, atom in enumerate(crystal):
if not atom.OBAtom.IsHydrogen():
bitvec.SetBitOn(i+1)
lookup.append(i)
success = pybel.ob.FindAutomorphisms(crystal.OBMol, mappings, bitvec)
# Find the RMSD between the crystal pose and each docked pose
xtalcoords = [atom.coords for atom in crystal if not atom.OBAtom.IsHydrogen()]
for i, dockedpose in enumerate(pybel.readfile("sdf", "dockedposes.sdf")):
posecoords = [atom.coords for atom in dockedpose if not atom.OBAtom.IsHydrogen()]
minrmsd = 999999999999
for mapping in mappings:
automorph_coords = [None] * len(xtalcoords)
for x, y in mapping:
automorph_coords[lookup.index(x)] = xtalcoords[lookup.index(y)]
mapping_rmsd = rmsd(posecoords, automorph_coords)
if mapping_rmsd < minrmsd:
minrmsd = mapping_rmsd
print("The RMSD is %.1f for pose %d" % (minrmsd, (i+1)))
@boyu
Copy link

boyu commented Jul 28, 2011

Hi, I run your code on my SuSE11.1 machine with open babel 2.3.0, and I got the following error:
Traceback (most recent call last):
File "dockedpose.py", line 21, in
mappings = pybel.ob.vvpairUIntUInt()
AttributeError: 'module' object has no attribute 'vvpairUIntUInt'

Could you please let me know which ob version you are suing? Thanks a lot!
Boyu

@mirix
Copy link

mirix commented Sep 13, 2011

I have compiled the latest development version of OpenBabel with the python bindings. However, I get the following error:

python dockedpose.py
Traceback (most recent call last):
File "dockedpose.py", line 2, in
import pybel
File "/usr/local/lib/python2.6/dist-packages/pybel.py", line 16, in
import openbabel as ob
File "/usr/local/lib/python2.6/dist-packages/openbabel.py", line 79, in
_openbabel = swig_import_helper()
File "/usr/local/lib/python2.6/dist-packages/openbabel.py", line 75, in swig_import_helper
_mod = imp.load_module('_openbabel', fp, pathname, description)
ImportError: dynamic module does not define init function (init_openbabel)

@filipsPL
Copy link

As IsHydrogen() is no longer in the API, this is the modified version of the script:

#!/usr/bin/env python

import math
import pybel

def squared_distance(coordsA, coordsB):
    """Find the squared distance between two 3-tuples"""
    sqrdist = sum( (a-b)**2 for a, b in zip(coordsA, coordsB) )
    return sqrdist

def rmsd(allcoordsA, allcoordsB):
    """Find the RMSD between two lists of 3-tuples"""
    deviation = sum(squared_distance(atomA, atomB) for
                    (atomA, atomB) in zip(allcoordsA, allcoordsB))
    return math.sqrt(deviation / float(len(allcoordsA)))

if __name__ == "__main__":
    print "Read crystal pose"
    crystal = next(pybel.readfile("sdf", "firstDocked.sdf"))

    # Find automorphisms involving only non-H atoms
    mappings = pybel.ob.vvpairUIntUInt()
    bitvec = pybel.ob.OBBitVec()
    lookup = []
    for i, atom in enumerate(crystal):
        if not atom.OBAtom.GetAtomicNum()==1:
            bitvec.SetBitOn(i+1)
            lookup.append(i)
    success = pybel.ob.FindAutomorphisms(crystal.OBMol, mappings, bitvec)

    print "Find the RMSD between the crystal pose and each docked pose"
    xtalcoords = [atom.coords for atom in crystal if not atom.OBAtom.GetAtomicNum()==1]
    for i, dockedpose in enumerate(pybel.readfile("sdf", "docking10.sdf")):
        posecoords = [atom.coords for atom in dockedpose if not atom.OBAtom.GetAtomicNum()==1]
        minrmsd = 999999999999
        for mapping in mappings:
            automorph_coords = [None] * len(xtalcoords)
            for x, y in mapping:
                automorph_coords[lookup.index(x)] = xtalcoords[lookup.index(y)]
            mapping_rmsd = rmsd(posecoords, automorph_coords)
            if mapping_rmsd < minrmsd:
                minrmsd = mapping_rmsd
        print("The RMSD is %.1f for pose %d" % (minrmsd, (i+1)))

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment