What's up dock? - Calculate the RMSD between docked and crystal poses

In cheminformatics, there are two reasons why one might want to calculate the RMSD between two conformers. The first is to check whether two conformers are very close in structure - e.g. for the purpose of generating a diverse set of conformers. This problem is solved using (1) a least squares alignment followed by (2) calculation of the RMSD.

The other situation is comparing two sets of 3D coordinates to see whether a prediction method has accurately reproduced experimental coordinates (e.g. docking). This just requires step (2) above.

The situation is complicated a little bit by the fact that only "heavy" atoms (i.e. non-H atoms in this context) are typically used to calculate the RMSD. A much greater complication is that automorphisms (well, isomorphisms of two molecules which are identical, to be exact) must be taken into account in both cases above. For example, consider the case where two para-substituted benzene rings must be compared; the RMSD calculation must take into account the fact that a 180 degree flip of the ring might yield a smaller RMSD.

Anyhoo, here's some Pybel code that will calculate the RMSD between a crystal pose and a set of docked poses. The code also illustrates how to access the isomorphisms. You should modify the code for your specific purpose:

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

view raw dockedpose.py hosted with ❤ by GitHub