Instantly share code, notes, and snippets.

# baoilleach/dockedpose.py Created May 16, 2011

What would you like to do?
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 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 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 commented Oct 12, 2017

 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)))```