SureChEMBL iPython Notebook Tutorial¶
An introduction to patent chemoinformatics using SureChEMBL data and the RDKit toolkit¶
George Papadatos, ChEMBL group, EMBL-EBI¶
In this tutorial:
- Read a file that contains all chemistry extracted from the Levitra US patent (US6566360) along with all the other members of the same patent family.
- Filter by different text-mining and chemoinformatics properties to remove noise and enrich the genuinely novel structures claimed in the patent documents.
- Visualise the chemical space using MDS and dimensionality reduction.
- Identify outliers in the chemistry space.
- Fix wrongly extracted structures.
- Find the Murcko and MCS scaffolds of the claimed compounds in the patent family.
- Compare the derived core with the actual Markush structure as reported in the original patent document.
- Identify key compounds using structural information only. This is based on the publications by Hattori et al. and Tyrchan et al..
- Count the number of NNs per compound.
- Perform R-Group decomposition.
import matplotlib.pyplot as plt
import pandas as pd
from IPython.display import display
from IPython.display import display_pretty, display_html, HTML, Javascript, Image
from rdkit.Chem import PandasTools
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import DataStructs
from rdkit.Chem import MCS
from rdkit.Chem import Draw
Draw.DrawingOptions.elemDict[0]=(0.,0.,0.) # draw dummy atoms in black
from itertools import cycle
from sklearn import manifold
from scipy.spatial.distance import *
import mpld3
mpld3.enable_notebook()
import warnings
warnings.filterwarnings('ignore')
from rdkit import rdBase
rdBase.rdkitVersion
This is the base URL of the publicly accessible ChEMBL Beaker server.
base_url = 'www.ebi.ac.uk/chembl/api'
pd.options.display.mpl_style = 'default'
pd.options.display.float_format = '{:.2f}'.format
rcParams['figure.figsize'] = 16,10
1. Read and filter the input file
The file was generated by extracting all chemistry from a list of patent documents in SureChEMBL. The Lucene query used to retrieve the list of relevant patents was: pn:"US6566360".
df = pd.read_csv('document chemistry_20141011_114421_271.csv',sep=',')
Let's check the contents:
df.info()
df.shape
Add a SCHEMBL ID column and sort:
df['SCHEMBL ID'] = df['Chemical ID'].map(lambda x: 'SCHEMBL{0}'.format(x))
df = df.sort('Chemical ID')
#df.head()
First round of filtering: Novel compounds appear in the description or claims section of the document. Alternatively, they are extracted from images or mol files
dff = df[(df['Claims Count'] > 0) | (df['Description Count'] > 0) | (df['Type'] != "TEXT")]
dff.shape
Second round of filtering: Simple physicochemical properties and counts
dff = dff[(dff['Rotatable Bound Count'] < 15) & (dff['Ring Count'] > 0) & (df['Radical'] == 0) & (df['Singleton'] == 0) & (df['Simple'] == 0)]
dff.shape
dff = dff[(dff['Molecular Weight'] >= 300) & (dff['Molecular Weight'] <= 800) & (dff['Log P'] > 0) & (dff['Log P'] < 7)]
dff.shape
Third round of filtering: Using Global Corpus Count
dff = dff[(dff['Chemical Corpus Count'] < 400) & ((dff['Annotation Corpus Count'] < 400) | (dff['Annotation Corpus Count'].isnull()))]
dff.shape
Keep largest fragment and convert SMILES to RDKit molecules
dff['SMILES'] = dff['SMILES'].map(lambda x: max(x.split('.'), key=len))
PandasTools.AddMoleculeColumnToFrame(dff, smilesCol = 'SMILES')
#dff.head(10)
Fourth round of filtering: Remove duplicates based on Chemical IDs and InChI keys
dff['InChI'] = dff['ROMol'].map(Chem.MolToInchi)
dff['InChIKey'] = dff['InChI'].map(Chem.InchiToInchiKey)
#dff.head()
dff = dff.drop_duplicates('Chemical ID')
dff.shape
dff = dff.drop_duplicates('InChIKey')
dff.shape
Wow, that was a lot of duplicates. This is because the same compounds come from different patents in the same patent family. In addition, in US patents a compound may come from 3 different sources: text, image and mol file.
Fifth round of filtering: Remove Boron-containing compounds as they are likely to be reactants.
dff = dff.ix[~(dff['ROMol'] >= Chem.MolFromSmarts('[#5]'))]
dff.shape
dff = dff.set_index('Chemical ID', drop=False)
That was the end of the filters - let's see a summary of the remaining compounds and their patent document source:
dff_counts = dff[['Patent ID','ROMol']].groupby('Patent ID').count()
dff_counts['Link'] = dff_counts.index.map(lambda x: '<a href="https://www.surechembl.org/document/{0}/" target="_blank">{0}</a>'.format(x))
dff_counts = dff_counts.rename(columns={'ROMol':'# Compounds'})
dff_counts
Right, is Levitra included in this set? ChEMBL tells me that its InChiKey is 'SECKRCOLJRRGGV-UHFFFAOYSA-N'.
dff.ix[dff.InChIKey == 'SECKRCOLJRRGGV-UHFFFAOYSA-N',['SCHEMBL ID','ROMol']]
| SCHEMBL ID | ROMol | |
|---|---|---|
| Chemical ID | ||
| 5441 | SCHEMBL5441 |  |
Good - it is...
2. Visualise the patent chemical space - MDS
Preparation of the pairwise distance matrix.
fps = [Chem.GetMorganFingerprintAsBitVect(m,2,nBits=2048) for m in dff['ROMol']]
dist_mat = squareform(pdist(fps,'jaccard'))
mds = manifold.MDS(n_components=2, dissimilarity="precomputed", random_state=3, n_jobs = 2)
results = mds.fit(dist_mat)
coords = results.embedding_
dff['X'] = [c[0] for c in coords]
dff['Y'] = [c[1] for c in coords]
csss = """
table
{
border-collapse: collapse;
}
th
{
color: #ffffff;
background-color: #848482;
}
td
{
background-color: #f2f3f4;
}
table, th, td
{
font-family:Arial, Helvetica, sans-serif;
border: 1px solid black;
text-align: right;
}
"""
import base64
dff['base64smi'] = dff['SMILES'].map(base64.b64encode)
Let's produce a scatter plot of the chemical space - points represent compounds, color-coded by the patent document they were found in. Thanks to mpld3, the scatter plot is interactive with live structure rendering calls to the local myChEMBL Beaker server.
fig, ax = plt.subplots()
fig.set_size_inches(14.0, 12.0)
#colors = cycle('bgrcmykwbgrcmykbgrcmykwbgrcmykw')
colors = cycle(cm.Dark2(np.linspace(0,1,len(dff_counts.index))))
for name, group in dff.groupby('Patent ID'):
labels = []
for i in group.index:
zz = group.ix[[i],['Patent ID','SCHEMBL ID','base64smi']]
zz['mol'] = zz['base64smi'].map(lambda x: '<img src="https://{0}/utils/smiles2image/{1}/300">'.format(base_url,x))
del zz['base64smi']
label = zz.T
del zz
label.columns = ['Row: {}'.format(i)]
labels.append(str(label.to_html()))
points = ax.scatter(group['X'], group['Y'],c=colors.next(), s=80, alpha=0.8)
tooltip = mpld3.plugins.PointHTMLTooltip(points, labels, voffset=10, hoffset=10, css=csss)
mpld3.plugins.connect(fig,tooltip)
Nice, we get Beaker rendering on the fly. However, there is a bunch of compounds on the left that are distinctively different. This is because of an image extraction error.
At least the error is systematic: extraction breaks the imidazole ring. Let's try to 'repair' the ring and rescue these compounds - we'll use an appropriate intramolecular reaction SMARTS.
resmarts = "([R:5][C;!R:1]=[C;H2].[C;H0,H1;!R:4]-[N;H2:2])>>[R:5][*:1]=[*:2]-[*:4]"
Draw.ReactionToImage(Chem.ReactionFromSmarts(resmarts))
def fixRing(mol,resmarts=resmarts):
rxn = Chem.ReactionFromSmarts(resmarts)
ps = rxn.RunReactants((mol,))
if ps:
m = ps[0][0]
Chem.SanitizeMol(m)
return m
else:
return mol
So this...
mm = dff.ix[12823029]['ROMol']
mm
...will become this:
fixRing(mm)
Let's apply the repair to all affected structures:
for i, p, m in zip(dff.index, dff['Patent ID'], dff['ROMol']):
if p == 'EP-2295436-A1':
dff.ix[i,'ROMol'] = fixRing(m)
dff['InChI'] = dff['ROMol'].map(Chem.MolToInchi)
dff['InChIKey'] = dff['InChI'].map(Chem.InchiToInchiKey)
dff = dff.drop_duplicates('InChIKey')
So after the fix, we have even less structures because some of the fixed ones already existed. Hopefully, however, we have less false positives too:
dff.shape
dff['RD_SMILES'] = dff['ROMol'].map(Chem.MolToSmiles)
dff['base64rdsmi'] = dff['RD_SMILES'].map(base64.b64encode)
We're going to re-calculate the distance matrix and plot the MCS just to double-check they were no side-effects from the applications of the reaction.
fps = [Chem.GetMorganFingerprintAsBitVect(m,2,nBits=2048) for m in dff['ROMol']]
dist_mat = squareform(pdist(fps,'jaccard'))
mds = manifold.MDS(n_components=2, dissimilarity="precomputed", random_state=3, n_jobs = 2)
results = mds.fit(dist_mat)
coords = results.embedding_
dff['X'] = [c[0] for c in coords]
dff['Y'] = [c[1] for c in coords]
fig, ax = plt.subplots()
fig.set_size_inches(14.0, 12.0)
#colors = cycle('bgrcmykwbgrcmykbgrcmykwbgrcmykw')
colors = cycle(cm.Dark2(np.linspace(0,1,len(dff_counts.index))))
for name, group in dff.groupby('Patent ID'):
labels = []
for i in group.index:
zz = group.ix[[i],['Patent ID','SCHEMBL ID','base64rdsmi']]
zz['mol'] = zz['base64rdsmi'].map(lambda x: '<img src="https://{0}/utils/smiles2image/{1}/300">'.format(base_url,x))
del zz['base64rdsmi']
label = zz.T
del zz
label.columns = ['Row: {}'.format(i)]
labels.append(str(label.to_html()))
points = ax.scatter(group['X'], group['Y'],c=colors.next(), s=80, alpha=0.8)
tooltip = mpld3.plugins.PointHTMLTooltip(points, labels, voffset=10, hoffset=10, css=csss)
mpld3.plugins.connect(fig,tooltip)
So the "island" of outliers has disappeared. There's only one obvious outlier on the far left. This might be a reactant which wasn't filtered out.
3. Define the MCS for this set of claimed molecules
A helper function, inspired by Greg Landrum's post here.
def MCS_Report(ms,atomCompare='any',**kwargs):
mcs = MCS.FindMCS(ms,atomCompare=atomCompare,timeout=120,**kwargs)
nAts = np.array([x.GetNumAtoms() for x in ms])
print 'Mean nAts: {0}, mcs nAts: {1}'.format(nAts.mean(),mcs.numAtoms)
print 'MCS SMARTS: {0}'.format(mcs.smarts)
mcsM = Chem.MolFromSmarts(mcs.smarts)
# find the common atoms
if atomCompare == 'any':
mcsM2 = Chem.MolFromSmiles(mcs.smarts,sanitize=False)
atNums=[-1]*mcsM.GetNumAtoms()
for m in ms:
match = m.GetSubstructMatch(mcsM)
if not match:
continue
for qidx,midx in enumerate(match):
anum = m.GetAtomWithIdx(midx).GetAtomicNum()
if atNums[qidx]==-1:
atNums[qidx]=anum
elif anum!=atNums[qidx]:
atNums[qidx]=0
for idx,atnum in enumerate(atNums):
if atnum>0:
mcsM.GetAtomWithIdx(idx).SetAtomicNum(atnum)
mcsM.UpdatePropertyCache(False)
Chem.SetHybridization(mcsM)
smi = Chem.MolToSmiles(mcsM,isomericSmiles=True)
print "MCS SMILES: {0}".format(smi)
img=Draw.MolToImage(mcsM,kekulize=False)
return mcs.smarts,smi,img,mcsM
Let's add the Murcko framework to our data frame:
PandasTools.AddMurckoToFrame(dff)
PandasTools.AddMoleculeColumnToFrame(dff,smilesCol='Murcko_SMILES', molCol='MurMol')
PandasTools.AlignToScaffold(dff)
dff[['ROMol','MurMol']].head()
| ROMol | MurMol | |
|---|---|---|
| Chemical ID | ||
| 5441 |  |
 |
| 5517 |  |
|
| 5530 |  |
|
| 5537 |  |
|
| 5606 |  |
 |
And now we'll extract the MCS core:
mols = list(dff.MurMol)
smarts,smi,img,mcsM = MCS_Report(mols,threshold=0.8,ringMatchesRingOnly=True)
img
4. Compare the extracted MCS core with the actual claimed Markush structure.
Given the MCS core, we'll do an R-Group decomposition first to define the R-Groups and their attachement points.
core = Chem.MolFromSmarts(smarts)
pind = []
rgroups = []
for i, m in zip(dff.index, dff.ROMol):
rgrps = Chem.ReplaceCore(m,core,labelByIndex=True)
if rgrps:
s = Chem.MolToSmiles(rgrps, True)
for e in s.split("."):
pind.append(e[1:4].replace('*','').replace(']',''))
rgroups.append(e)
#print "R-Groups: ", i, s
pind = sorted(list(set(pind)),key=int)
pind = map(int,pind)
Draw.MolToImage(core,includeAtomNumbers=True,highlightAtoms=pind)
Is the MCS similar to the claimed Markush structure? Let's fetch the original document in PDF and check in page 4.
Pretty close!
5. Identify key compounds - in this case, vardenafil itself
This is based on the publication by Hattori et al. and Tyrchan et al.. First we'll try the Hattori et al. version of finding the compounds with the highest number of nearest neighbours (NNs) for a given similarity threshold. We will need to calculate the similarity matrix.
sim_mat = []
CUTOFF = 0.8
for i,fp in enumerate(fps):
tt = DataStructs.BulkTanimotoSimilarity(fps[i],fps)
for i,t in enumerate(tt):
if tt[i] < CUTOFF:
tt[i] = None
sim_mat.append(tt)
sim_mat=numpy.array(sim_mat)
sim_df = pd.DataFrame(sim_mat, columns = dff['Chemical ID'], index=dff['Chemical ID'])
sim_df.head()
d = sim_df.count().to_dict()
nn_df = pd.DataFrame(d.items(),columns=['Chemical ID', '# NNs'])
nn_df.set_index('Chemical ID',inplace=True)
nn_df.head()
pd.merge(dff[['Chemical ID','ROMol']],nn_df,left_index=True, right_index=True, sort=False).sort('# NNs',ascending=False).head()
| Chemical ID | ROMol | # NNs | |
|---|---|---|---|
| Chemical ID | |||
| 12823654 | 12823654 |  |
23 |
| 5441 | 5441 | |
22 |
| 12824064 | 12824064 |  |
21 |
| 5537 | 5537 | |
21 |
| 12822873 | 12822873 |  |
18 |
So, the method seems to be working in this case. Vardenafil has the second largest number of NNs with a cutoff of 0.8. Next, we'll try one of the Tyrchan et al. methods (Frequency of Groups, FOG) and count the frequency of the derived R-Groups:
from collections import Counter
rgroups = [r.replace('[15*]','[11*]').replace('[14*]','[12*]') for r in rgroups]
c = Counter(rgroups)
rgroup_counts = pd.DataFrame(c.most_common(20),columns=['RGroups','Frequency'])
rgroup_counts.head(10)
PandasTools.AddMoleculeColumnToFrame(rgroup_counts, smilesCol = 'RGroups')
rgroup_counts.head(10)
| RGroups | Frequency | ROMol | |
|---|---|---|---|
| 0 | [6*]CCC | 349 |  |
| 1 | [11*]OCC | 306 |  |
| 2 | [8*]C | 294 |  |
| 3 | [8*]CC | 63 |  |
| 4 | [11*]OCCC | 28 |  |
| 5 | [11*]O | 16 |  |
| 6 | [6*]C | 10 |  |
| 7 | [12*]S(=O)(=O)N1CCN(C)CC1 | 8 |  |
| 8 | [13*]OC | 8 |  |
| 9 | [12*]S(=O)(=O)N(CC)CCO | 7 |  |
rgroup_counts[[r.startswith('[12*]') for r in rgroup_counts.RGroups]].head(10)
| RGroups | Frequency | ROMol | |
|---|---|---|---|
| 7 | [12*]S(=O)(=O)N1CCN(C)CC1 | 8 | |
| 9 | [12*]S(=O)(=O)N(CC)CCO | 7 | |
| 10 | [12*]S(=O)(=O)Cl | 6 |  |
| 11 | [12*]S(=O)(=O)NCCCN1CCOCC1 | 4 |  |
| 12 | [12*]S(=O)(=O)N1CCN(CCO)CC1 | 4 |  |
| 13 | [12*]S(=O)(=O)N1CCC(O)CC1 | 4 |  |
| 14 | [12*]S(=O)(=O)N(CC)CC | 4 |  |
| 15 | [12*]S(=O)(=O)N1CCN(CC)CC1 | 4 |  |
| 16 | [12*]S(=O)(=O)N1CCC(CO)CC1 | 3 |  |
| 17 | [12*]S(=O)(=O)NCCOC | 3 |  |
Although the most frequent R-groups in positions 6, 8 and 11 correspond to those in Vardenafil (see below), the exact piperazine R-group in position 12 is the 7th most frequent one in that position.
dff.ix[dff.InChIKey == 'SECKRCOLJRRGGV-UHFFFAOYSA-N',['SCHEMBL ID','ROMol']]
| SCHEMBL ID | ROMol | |
|---|---|---|
| Chemical ID | ||
| 5441 | SCHEMBL5441 | |