steeveslab-blog

Multiple conformations of a molecule in RDKit

In [1]:
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions
In [2]:
import cPickle as pickle
In [3]:
from rdkit.Chem import PyMol
In [4]:
ibuH = pickle.load(open('ibuH.pkl','rb'))
In [5]:
ibuH
Out[5]:

Previously we built a single conformation of this molecule. We'd often like to have more than that to start to capture conformational heterogeneity.

In [6]:
cids = AllChem.EmbedMultipleConfs(ibuH, 
                                  clearConfs=True, 
                                  numConfs=100, 
                                  pruneRmsThresh=1)
In [7]:
print len(cids)
97

In [8]:
for cid in cids: AllChem.MMFFOptimizeMolecule(ibuH,confId=cid)
In [9]:
v= PyMol.MolViewer()
In [10]:
v.DeleteAll()
for cid in cids: 
    v.ShowMol(ibuH,confId=cid,name='Conf-%d'%cid,showOnly=False)
In [11]:
v.server.do('set grid_mode, on')
In [12]:
v.server.do('ray')
v.GetPNG()
Out[12]:

It's not clear from this picture if we have diverse conformations or are just rotating the molecule in space. We'd like to align of the molecules in order to be able to make this assessment. Specifically, let's choose to align them all, so the benzene rings are overlapping. We can select the benzene ring using an object called a SMARTS pattern:

In [13]:
patt = Chem.MolFromSmarts('c1ccccc1');patt #this is just a SMILES string, SMILES are valid SMARTS 
Out[13]:

Find which atoms ids in our molecule match the pattern described above:

In [14]:
match = ibuH.GetSubstructMatch(patt)
In [15]:
match
Out[15]:
(4, 5, 6, 7, 8, 9)
In [16]:
DrawingOptions.includeAtomNumbers=True
ibuH
Out[16]:

Now we can align of the conformers on the identified subset of atoms:

In [17]:
AllChem.AlignMolConformers(ibuH,atomIds=match)

And visualize:

In [18]:
v.DeleteAll()
for cid in cids: 
    v.ShowMol(ibuH,confId=cid,name='Conf-%d'%cid,showOnly=False)
In [19]:
v.server.do('ray')
v.GetPNG()
Out[19]:
In [20]:
pickle.dump(ibuH,open('ibuH_confs.pkl','wb'))