Quite a while ago I did a blog post showing how to use the Python function AllChem.ConstrainedEmbed() to generate conformations where the positions of a set of atoms are constrained to match the coordinates of a template molecule. More recently a question came up on the mailing list about how to use the core embedding functionality that lies underneath AllChem.ConstrainedEmbed(); this short post will look at that.
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem.Draw import rdDepictor
rdDepictor.SetPreferCoordGen(True)
IPythonConsole.molSize = (250,200)
from rdkit.Chem import AllChem
import math
import rdkit
import time
print(rdkit.__version__)
print(time.asctime())
We're going to work with cyclosporine here, since it's delightfully complicated, which makes it fun. :-) In this case we're going to ignore atomic stereochemistry in order to speed the conformation generation up (the RDKit tends to take a while to generate conformations for molecules with a large number of stereocenters).
We'll also define the macrocylce as the core; this is what we're going to use to provide constraints.
m = Chem.MolFromSmiles('C/C=C/CC(C)C(O)C1C(=O)NC(CC)C(=O)N(C)CC(=O)N(C)C(CC(C)C)C(=O)NC(C(C)C)C(=O)N(C)C(CC(C)C)C(=O)NC(C)C(=O)NC(C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(CC(C)C)C(=O)N(C)C(C(C)C)C(=O)N1C')
rdDepictor.Compute2DCoords(m)
core = Chem.MolFromSmiles('C1C(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)NCC(=O)N1')
m.GetSubstructMatch(core)
m
Start by generating a conformation for cyclosporine itself.
mh = Chem.AddHs(m)
AllChem.EmbedMolecule(mh,randomSeed=0xf00d)
For the other molecule, we'll take a molecule I pulled from ChEMBL that includes the core. I also removed stereochemistry here.
smi = 'CCC1NC(=O)C(C(O)C(C)CCCC(=O)OC)N(C)C(=O)C(C(C)C)N(C)C(=O)C(CC(C)C)N(C)C(=O)C(CC(C)C)N(C)C(=O)C(C)NC(=O)C(C)NC(=O)C(CC(C)C)N(C)C(=O)C(C(C)C)NC(=O)C(CC(C)C)N(C)C(=O)CN(C)C1=O'
newm = Chem.MolFromSmiles(smi)
rdDepictor.Compute2DCoords(newm)
newm.GetSubstructMatch(core)
newm
Generate an unconstrained conformation for that:
newmh = Chem.AddHs(newm)
AllChem.EmbedMolecule(newmh,randomSeed=0xf00d)
Now calculate the RMSD between the core atoms in cylosporine and in our test molecule.
m_match = mh.GetSubstructMatch(core)
newm_match = newmh.GetSubstructMatch(core)
delta2 = 0.0
for mi,newmi in zip(m_match,newm_match):
d = (mh.GetConformer().GetAtomPosition(mi) - newmh.GetConformer().GetAtomPosition(newmi)).LengthSq()
delta2 += d
print('core RMSD:',math.sqrt(delta2/len(m_match)))
There's no reason to expect this RMSD to be anything other than huge: we're using different conformations of a flexible core and haven't aligned them to each other.
We can go ahead and do an alignmeent and see how that affects the RMSD:
AllChem.AlignMol(newmh,mh,atomMap = list(zip(newm_match,m_match)))
delta2 = 0.0
for mi,newmi in zip(m_match,newm_match):
d = (mh.GetConformer().GetAtomPosition(mi) - newmh.GetConformer().GetAtomPosition(newmi)).LengthSq()
delta2 += d
print('aligned core RMSD:',math.sqrt(delta2/len(m_match)))
Again, it's not surprising that this is a large RMSD: the core is quite flexible and we haven't constrained it at all.
Adding constraints is the point of this blog post though, so let's move onto that.
AllChem.EmbedMolecule() can be given constraints for the positions of certain atoms using the coordMap argument, which expects a dictionary that provides a Point3D for each atom that should have a fixed position.
cmap = {newm_match[i]:mh.GetConformer().GetAtomPosition(m_match[i]) for i in range(len(m_match))}
AllChem.EmbedMolecule(newmh,randomSeed=0xf00d,coordMap=cmap)
delta2 = 0.0
for mi,newmi in zip(m_match,newm_match):
d = (mh.GetConformer().GetAtomPosition(mi) - newmh.GetConformer().GetAtomPosition(newmi)).LengthSq()
delta2 += d
print('core RMSD:',math.sqrt(delta2/len(m_match)))
Wait... what happened here? Shouldn't this number be smaller since we introduced constraints? This was the question on the mailing list.
The function AllChem.EmbedMolecule() uses the coordinates provided in coordMap to set elements of the distance bounds matrix that is used to generate conformations (details about the RDKit's distance-geometry-based conformation generator are in the documentation). This results in conformations where the d
Note that, because the coordinates are being constrained using the distances between them, you should expect rigid shifts of the core atoms relative to the constraints. This is solveable by aligning the core of the test molecule to the core of the reference:
AllChem.AlignMol(newmh,mh,atomMap = list(zip(newm_match,m_match)))
delta2 = 0.0
for mi,newmi in zip(m_match,newm_match):
d = (mh.GetConformer().GetAtomPosition(mi) - newmh.GetConformer().GetAtomPosition(newmi)).LengthSq()
delta2 += d
print('aligned core RMSD:',math.sqrt(delta2/len(m_match)))
Note that the output coordinates don't match the constraint coordinates exactly. This will almost always be the case; they should be close, but some differences are, unforunately, expected due to the nature of the algorithm.
Aside 1: while thinking about this blog post, I came up with a relatively straightforward way to force the coordinates to match exactly (and to remove the need for the alignment step) in some circumstances. Assuming it works as I expect it to, this will hopefully end up in the 2019.03 RDKit release.
Aside 2: if it's really important to you that the constrained coordinates exactly match the constraints, take a look at AllChem.ConstrainedEmbed and that old blog post mentioned above. AllChem.ConstrainedEmbed does a couple of extra steps to ensure that the output coordinates of the constrained atoms match the constraint positions very closely.
It's worth taking a look at the conformations of our molecules. We'll do that using py3Dmol. There's some more info on using this with the RDKit in another old blog post.
import py3Dmol
def drawit(ms, p=None, confId=-1):
if p is None:
p = py3Dmol.view(width=400, height=400)
p.removeAllModels()
for m in ms:
mb = Chem.MolToMolBlock(m, confId=confId)
p.addModel(mb, 'sdf')
p.setStyle({'stick':{}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
return p.show()
drawit((mh,newmh))
