-
-
Save GeauxEric/4a540c8e13321954d2f3 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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))) | |
| def rmsd_between(mol1, mol2): | |
| """Find the RMSD between two molecules | |
| Keyword Arguments: | |
| mol1 -- first pybel molecule | |
| mol2 -- second pybel molecule | |
| """ | |
| mappings = pybel.ob.vvpairUIntUInt() | |
| bitvec = pybel.ob.OBBitVec() | |
| lookup = [] | |
| for i, atom in enumerate(mol1): | |
| if not atom.OBAtom.IsHydrogen(): | |
| bitvec.SetBitOn(i + 1) | |
| lookup.append(i) | |
| _ = pybel.ob.FindAutomorphisms(mol1.OBMol, mappings, bitvec) | |
| xtalcoords = [atom.coords for atom in mol1 if not atom.OBAtom.IsHydrogen()] | |
| posecoords = [atom.coords for atom in mol2 if not atom.OBAtom.IsHydrogen()] | |
| minrmsd = 99999999 | |
| 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 | |
| return minrmsd | |
| 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))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment