Goal: construct a set of molecular pairs that can be used to compare similarity methods to each other.
I want to start with molecules that have some connection to each other, so I will pick pairs that have a baseline similarity: a Tanimoto similarity using count based Morgan0 fingerprints of at least 0.7. This threshold was selected empirically.
Note: this notebook and the data it uses/generates can be found in the github repo: https://github.com/greglandrum/rdkit_blog
I'm going to use ChEMBL as my data source, so I'll start by adding a table with Morgan0 fingerprints to my local copy of ChEMBL 16:
chembl_16=# select molregno,morganbv_fp(m,0) as mfp0 into rdk.tfps from rdk.mols;
SELECT 1291111
chembl_16=# create index fps_mfp0_idx on rdk.tfps using gist(mfp0);
CREATE INDEX
And now create a smaller table that only contains molecules with molwt<600:
chembl_16=# select molregno,mfp0 into table rdk.tfps_smaller from rdk.tfps join rdk.props using (molregno) where amw<600;
SELECT 1182986
chembl_16=# create index sfps_mfp0_idx on rdk.tfps_smaller using gist(mfp0);
CREATE INDEX
And now I'll build the set of pairs using Python. This is definitely doable in SQL, but my SQL-fu isn't that strong.
Start by getting a set of 35K random small molecules with MW<600:
import psycopg2
cn = psycopg2.connect(dbname='chembl_16')
curs = cn.cursor()
curs.execute('select molregno,m from rdk.mols join rdk.props using (molregno) where amw<=600 order by random() limit 35000')
qs = curs.fetchall()
And now find one neighbor for 25K of those from the mfp0 table of smallish molecules:
curs.execute('set rdkit.tanimoto_threshold=0.7')
keep=[]
for i,row in enumerate(qs):
curs.execute('select molregno,m from rdk.mols join (select molregno from rdk.tfps_smaller where mfp0%%morgan_fp(%s,0) and molregno!=%s limit 1) t2 using (molregno)',(row[1],row[0]))
d = curs.fetchone()
if not d: continue
keep.append((row[0],row[1],d[0],d[1]))
if len(keep)==25000: break
if not i%500: print 'Done: %d'%i
Done: 0 Done: 500 Done: 1000 Done: 1500 Done: 2000 Done: 2500 Done: 3000 Done: 3500 Done: 4000 Done: 4500 Done: 5000 Done: 5500 Done: 6000 Done: 6500 Done: 7000 Done: 7500 Done: 8000 Done: 8500 Done: 9000 Done: 9500 Done: 10000 Done: 10500 Done: 11000 Done: 11500 Done: 12000 Done: 12500 Done: 13000 Done: 13500 Done: 14000 Done: 14500 Done: 15000 Done: 15500 Done: 16000 Done: 16500 Done: 17000 Done: 17500 Done: 18000 Done: 18500 Done: 19000 Done: 19500 Done: 20000 Done: 20500 Done: 21000 Done: 21500 Done: 22000 Done: 22500 Done: 23000 Done: 23500 Done: 24000 Done: 24500 Done: 25000
Finally, write those out to a file so that we can use them elsewhere:
import gzip
outf = gzip.open('../data/chembl16_25K.pairs.txt.gz','wb+')
for idx1,smi1,idx2,smi2 in keep: outf.write('%d %s %d %s\n'%(idx1,smi1,idx2,smi2))
outf=None
Start by loading the pairs from the file we saved and creating RDKit molecules from them
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
import gzip
rows=[]
for row in gzip.open('../data/chembl16_25K.pairs.txt.gz').readlines():
row = row.split()
row[1] = Chem.MolFromSmiles(row[1])
row[3] = Chem.MolFromSmiles(row[3])
rows.append(row)
Look at some pairs:
t = []
for x in rows[:5]:
t.append(x[1])
t.append(x[3])
Draw.MolsToGridImage(t,molsPerRow=2)
Each plot below contains two histograms. The one in blue is for the first set of molecules, the one in green is for the neighbor molecules.
from rdkit.Chem import Descriptors
mws = [(Descriptors.MolWt(x[1]),Descriptors.MolWt(x[3])) for x in rows]
nrots = [(Descriptors.NumRotatableBonds(x[1]),Descriptors.NumRotatableBonds(x[3])) for x in rows]
logps = [(Descriptors.MolLogP(x[1]),Descriptors.MolLogP(x[3])) for x in rows]
_=hist(([x for x,y in mws],[y for x,y in mws]),bins=20,histtype='bar')
xlabel('AMW')
<matplotlib.text.Text at 0x8d6e6890>
_=hist(([x for x,y in logps],[y for x,y in logps]),bins=20,histtype='bar')
xlabel('mollogp')
<matplotlib.text.Text at 0x47288310>
_=hist(([x for x,y in nrots],[y for x,y in nrots]),bins=20,histtype='bar')
xlabel('num rotatable bonds')
<matplotlib.text.Text at 0x56ef06d0>
and a histogram of the similarities we used to construct the set
from rdkit import DataStructs
from rdkit.Chem import rdMolDescriptors
sims = [DataStructs.TanimotoSimilarity(rdMolDescriptors.GetMorganFingerprint(x[1],0),rdMolDescriptors.GetMorganFingerprint(x[3],0)) for x in rows]
_=hist(sims,bins=20)
xlabel('MFP0 sims within pairs')
<matplotlib.text.Text at 0x472ae0d0>
compare to MFP2 similarity (more on this in a later post)
sims2 = [DataStructs.TanimotoSimilarity(rdMolDescriptors.GetMorganFingerprint(x[1],2),rdMolDescriptors.GetMorganFingerprint(x[3],2)) for x in rows]
_=scatter(sims,sims2,marker='o',edgecolors='none')
xlabel('MFP0 sim')
ylabel('MFP2 sim')
<matplotlib.text.Text at 0x7d634250>
Look at the distribution of MFP0 similarities in random molecule pairs (more on this in a later post)
import random
idxs = list(range(len(rows)))
random.shuffle(idxs)
ms1 = [x[1] for x in rows]
ms2 = [rows[x][3] for x in idxs]
sims = [DataStructs.TanimotoSimilarity(rdMolDescriptors.GetMorganFingerprint(x,0),rdMolDescriptors.GetMorganFingerprint(y,0)) for x,y in zip(ms1,ms2)]
_=hist(sims,bins=20)
xlabel('MFP0 sim in random pairs')
<matplotlib.text.Text at 0x452b390>