Saturday, August 8, 2015

Impact of the threshold on similarity search times

Update: added timing information for chemfp with the larger ZINC set.

A ChEMBL blog post included the (to me) somewhat surprising result that similarity search times while using the RDKit cartridge did not show much of a dependency on the similarity threshold being used. I wanted to investigate this a bit more closely.

I'm using a local install of ChEMBL20 for these tests. Since I don't have the patience to wait for searches to complete with 1000 molecules, I will randomly select just 10:

chembl_20=# select * into temporary table foo from rdk.fps order by random() limit 10;
SELECT 10
chembl_20=# select molregno,m from rdk.mols join foo using (molregno);
 molregno |                               m                                
----------+----------------------------------------------------------------
  1067143 | C[C@@H](O)[C@@H](NC(=O)N1CCN(c2ccc(C#Cc3ccccc3)cc2)CC1)C(=O)NO
   552129 | Cl.Cl.NCCNCCSc1ccncc1
  1034319 | CC(C)c1ccc(CN2CC(C(=O)Nc3ccccc3)CC2=O)cc1
   244180 | NS(=O)(=O)c1ccc(-n2nc(C(F)(F)F)cc2-c2cccc(Cl)c2)cc1CO
   357411 | CCCC1=C2c3ccc4[nH]nnc4c3CC2(CCC)CCC1=O
  1116827 | CCc1ccc(C2=NN(CCC(=O)NC3CC3)C(=O)CC2)cc1
   783512 | CCc1ccc(-c2cc(C)no2)cc1S(=O)(=O)NCc1ccc(OC)cc1
   623138 | Nc1ccc(-c2cn([C@@H]3O[C@H](CO)[C@H](O)[C@H](O)[C@H]3O)nn2)cc1
  1109930 | CCc1ccc(OCC(=O)OCC(=O)NCCNC(=O)COC(=O)COc2ccc(CC)cc2)cc1
   766599 | CC(=O)c1ccc(NC(=O)C(C)OC(=O)COc2ccc(C#N)cc2)cc1
(10 rows)

And now look at timing results for a number of threshold values:

chembl_20=# set rdkit.tanimoto_threshold=0.4;
SET
Time: 0.104 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
  3142
(1 row)

Time: 7022.514 ms
chembl_20=# set rdkit.tanimoto_threshold=0.5;
SET
Time: 0.103 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
   593
(1 row)

Time: 6966.116 ms
chembl_20=# set rdkit.tanimoto_threshold=0.6;
SET
Time: 0.112 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
   201
(1 row)

Time: 6914.092 ms
chembl_20=# set rdkit.tanimoto_threshold=0.7;
SET
Time: 0.116 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    93
(1 row)

Time: 6961.017 ms
chembl_20=# set rdkit.tanimoto_threshold=0.8;
SET
Time: 0.127 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    32
(1 row)

Time: 6952.292 ms
chembl_20=# set rdkit.tanimoto_threshold=0.9;
SET
Time: 0.108 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    12
(1 row)

Time: 6700.883 ms
chembl_20=# set rdkit.tanimoto_threshold=0.95;
SET
Time: 0.106 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    10
(1 row)

Time: 5512.117 ms
chembl_20=# set rdkit.tanimoto_threshold=0.99;
SET
Time: 0.102 ms
chembl_20=# select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    10
(1 row)

Time: 2315.399 ms

It strikes me as somewhat odd that things are more or less flat until very high threshold values are reached; 0.95 is the first one where a real difference is observable.

Let's look at the results of explain analyze to see if that helps understand:

chembl_20=# set rdkit.tanimoto_threshold=0.4;
SET
Time: 0.138 ms
chembl_20=# explain (analyze on, buffers on) select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;                                                                  QUERY PLAN                                                                   
-----------------------------------------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=291843.41..291843.42 rows=1 width=0) (actual time=6975.350..6975.351 rows=1 loops=1)
   Buffers: shared hit=253423, local hit=1
   ->  Nested Loop  (cost=0.43..289586.70 rows=902685 width=0) (actual time=43.666..6975.065 rows=3142 loops=1)
         Buffers: shared hit=253423, local hit=1
         ->  Seq Scan on foo  (cost=0.00..16.20 rows=620 width=32) (actual time=0.003..0.015 rows=10 loops=1)
               Buffers: local hit=1
         ->  Index Scan using fps_mfp2_idx on fps fps1  (cost=0.43..452.49 rows=1456 width=65) (actual time=72.735..697.462 rows=314 loops=10)
               Index Cond: (mfp2 % foo.mfp2)
               Buffers: shared hit=253423
 Planning time: 0.063 ms
 Execution time: 6975.381 ms
(11 rows)

Time: 6975.672 ms
chembl_20=# set rdkit.tanimoto_threshold=0.7;
SET
Time: 0.118 ms
chembl_20=# explain (analyze on, buffers on) select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
                                                                  QUERY PLAN                                                                  
----------------------------------------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=291843.41..291843.42 rows=1 width=0) (actual time=7032.295..7032.296 rows=1 loops=1)
   Buffers: shared hit=250500, local hit=1
   ->  Nested Loop  (cost=0.43..289586.70 rows=902685 width=0) (actual time=240.418..7032.260 rows=93 loops=1)
         Buffers: shared hit=250500, local hit=1
         ->  Seq Scan on foo  (cost=0.00..16.20 rows=620 width=32) (actual time=0.004..0.015 rows=10 loops=1)
               Buffers: local hit=1
         ->  Index Scan using fps_mfp2_idx on fps fps1  (cost=0.43..452.49 rows=1456 width=65) (actual time=238.074..703.216 rows=9 loops=10)
               Index Cond: (mfp2 % foo.mfp2)
               Buffers: shared hit=250500
 Planning time: 0.063 ms
 Execution time: 7032.327 ms
(11 rows)

Time: 7032.634 ms
chembl_20=# set rdkit.tanimoto_threshold=0.9;
SET
Time: 0.139 ms
chembl_20=# explain (analyze on, buffers on) select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
                                                                  QUERY PLAN                                                                  
----------------------------------------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=291843.41..291843.42 rows=1 width=0) (actual time=6768.003..6768.004 rows=1 loops=1)
   Buffers: shared hit=238999, local hit=1
   ->  Nested Loop  (cost=0.43..289586.70 rows=902685 width=0) (actual time=243.949..6767.985 rows=12 loops=1)
         Buffers: shared hit=238999, local hit=1
         ->  Seq Scan on foo  (cost=0.00..16.20 rows=620 width=32) (actual time=0.004..0.018 rows=10 loops=1)
               Buffers: local hit=1
         ->  Index Scan using fps_mfp2_idx on fps fps1  (cost=0.43..452.49 rows=1456 width=65) (actual time=345.900..676.791 rows=1 loops=10)
               Index Cond: (mfp2 % foo.mfp2)
               Buffers: shared hit=238999
 Planning time: 0.062 ms
 Execution time: 6768.034 ms
(11 rows)

Time: 6768.345 ms
chembl_20=# set rdkit.tanimoto_threshold=0.99;
SET
Time: 0.100 ms
chembl_20=# explain (analyze on, buffers on) select count(*) from rdk.fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
                                                                  QUERY PLAN                                                                  
----------------------------------------------------------------------------------------------------------------------------------------------
 Aggregate  (cost=291843.41..291843.42 rows=1 width=0) (actual time=2287.224..2287.224 rows=1 loops=1)
   Buffers: shared hit=76686, local hit=1
   ->  Nested Loop  (cost=0.43..289586.70 rows=902685 width=0) (actual time=87.785..2287.213 rows=10 loops=1)
         Buffers: shared hit=76686, local hit=1
         ->  Seq Scan on foo  (cost=0.00..16.20 rows=620 width=32) (actual time=0.003..0.013 rows=10 loops=1)
               Buffers: local hit=1
         ->  Index Scan using fps_mfp2_idx on fps fps1  (cost=0.43..452.49 rows=1456 width=65) (actual time=127.327..228.715 rows=1 loops=10)
               Index Cond: (mfp2 % foo.mfp2)
               Buffers: shared hit=76686
 Planning time: 0.063 ms
 Execution time: 2287.256 ms
(11 rows)

Time: 2287.559 ms

And that's it, the number of buffers that need to be accessed ends up being essentially the same until getting to the very high thresholds. As to why this is true? That I don't really know; I suspect figuring it out would require digging more deeply into the innards of how the GIST indexes work in PostgreSQL, and that's not something I'm going to do right now.

A larger test: ZINC

To see if the same trends are there for a larger dataset, I grabbed a copy of the ZINC All Clean set. After removing the molecules that the RDKit doesn't like (there are actually only 14 molecules in the full set that the RDKit rejected), this leaves 16.4 million molecules.

For the record, here's how I built the database:

glandrum@Otter:/scratch/RDKit_git/Data/Zinc$ createdb zinc
glandrum@Otter:/scratch/RDKit_git/Data/Zinc$ psql -c 'create extension rdkit' zinc
CREATE EXTENSION
glandrum@Otter:/scratch/RDKit_git/Data/Zinc$ psql -c 'create table raw_data (id SERIAL, smiles text, zinc_id char(12))' zinc
CREATE TABLE
glandrum@Otter:/scratch/RDKit_git/Data/Zinc$ zcat zinc_all_clean.smi.gz | sed '1d; s/\\/\\\\/g' |psql -c "copy raw_data (smiles,zinc_id) from stdin with delimiter ' '" zinc
COPY 16403864
glandrum@Otter:/scratch/RDKit_git/Data/Zinc$ psql zinc
psql (9.4.2)
Type "help" for help.

zinc=# \timing
Timing is on.
zinc=# select * into mols from (select id,mol_from_smiles(smiles::cstring) m from raw_data) tmp where m is not null;
SELECT 16403848
Time: 8591329.609 ms
zinc=# select id,morganbv_fp(m) as mfp2 into fps from mols;
SELECT 16403848
Time: 1308215.512 ms
zinc=# create index fps_mfp2_idx on fps using gist(mfp2);
CREATE INDEX
Time: 1230138.847 ms
zinc=# select * into temporary table foo from fps order by random() limit 10;
SELECT 10
Time: 3640.389 ms

Pick 10 random molecules:

zinc=# select id,m from mols join foo using (id);
    id    |                                   m                                    
----------+------------------------------------------------------------------------
 11373174 | COc1cccc([C@@H](C)CC(=O)NCCCOC2CCCC2)c1
  6068083 | COc1ccc(CC[C@@H](C)[NH+]2CCC[C@H](O)C2)cc1
 10184316 | COc1cccc(O[C@@H](C)CNc2nccn(C)c2=O)c1
   869387 | COc1ccc(NC(=O)CN(Cc2ccc(F)cc2)C(=O)Cn2nnc3ccccc32)cc1
  2996517 | CN1C(=O)[C@H](NC(=O)Nc2ccc(Cl)cc2)N=C(c2ccccc2)c2ccccc21
  2749280 | CC(C)NS(=O)(=O)Cc1ccc(CNC(=O)c2ccc(CSc3nc4ccccc4s3)cc2)cc1
  5014876 | CCCCOC(=O)[C@@H](O)CC
 10653280 | Cc1cc(NC(=O)NC[C@](C)(O)c2ccc(C)o2)n(Cc2ccccn2)n1
  1253806 | CCCC[C@@H](CC)CNC(=O)c1cccc(OC)c1
  4846692 | O=C(c1ccc(N2CC[NH2+]CC2)c(NS(=O)(=O)c2ccccc2)c1)N1CC[NH+](C2CCCCC2)CC1
(10 rows)

Time: 161.055 ms

Now the searches:

zinc=# set rdkit.tanimoto_threshold=0.5;
SET
Time: 0.141 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2; count 
-------
  7443
(1 row)

Time: 80237.225 ms
zinc=# set rdkit.tanimoto_threshold=0.6;
SET
Time: 0.089 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
   875
(1 row)

Time: 79778.422 ms
zinc=# set rdkit.tanimoto_threshold=0.7;
SET
Time: 0.089 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
   144
(1 row)

Time: 79741.891 ms
zinc=# set rdkit.tanimoto_threshold=0.8;
SET
Time: 0.092 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    44
(1 row)

Time: 79128.903 ms
zinc=# set rdkit.tanimoto_threshold=0.9;
SET
Time: 0.087 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    21
(1 row)

Time: 73175.759 ms
zinc=# set rdkit.tanimoto_threshold=0.95;
SET
Time: 0.108 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    20
(1 row)

Time: 49995.441 ms
zinc=# set rdkit.tanimoto_threshold=0.99;
SET
Time: 0.105 ms
zinc=# select count(*) from fps fps1 cross join (select * from foo) tmp where fps1.mfp2%tmp.mfp2;
 count 
-------
    18
(1 row)

Time: 17011.942 ms

This time we start to see a difference at a threshold of 0.9, but things remain more or less constant until then.

An aside: chemfp performance

For raw search speed, I'm not aware of anything that does better than Andrew Dalke's chemfp, so I figured I'd go ahead and see what exceptional performance looks like with the ChEMBL20 dataset.

I start by creating two FPB files with Morgan2 fingerprints:

(py27)glandrum@Otter:~/RDKit_blog/data$ rdkit2fps -o chembl_20.mfp2.fpb --morgan ~/Downloads/chembl_20.sdf.gz (py27)glandrum@Otter:~/RDKit_blog/data$ rkit2fps -o chembl_random_10.fbp --morgan chembl_random_10.smi

Note that I am using the licensed version of chemfp for these tests. The open-source version is also super fast, but the idea here is to see what the best I can do is:

(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.4 -q chembl_random_10.fpb --times -k all -c chembl_20.mfp2.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.4
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=chembl_20.mfp2.fpb
#query_sources=chembl_random_10.smi
#target_sources=/home/glandrum/Downloads/chembl_20.sdf.gz
30  552129
429 1109930
232 623138
433 1034319
125 766599
66  1067143
99  1116827
201 783512
49  357411
179 244180
open 0.00 search 0.12 total 0.12
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.9 -q chembl_random_10.fpb --times -k all -c chembl_20.mfp2.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.9
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=chembl_20.mfp2.fpb
#query_sources=chembl_random_10.smi
#target_sources=/home/glandrum/Downloads/chembl_20.sdf.gz
2   552129
1   1109930
1   623138
1   1034319
1   766599
1   1067143
2   1116827
1   783512
1   357411
1   244180
open 0.00 search 0.05 total 0.05
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.95 -q chembl_random_10.fpb --times -k all -c chembl_20.mfp2.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.95
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=chembl_20.mfp2.fpb
#query_sources=chembl_random_10.smi
#target_sources=/home/glandrum/Downloads/chembl_20.sdf.gz
1   552129
1   1109930
1   623138
1   1034319
1   766599
1   1067143
1   1116827
1   783512
1   357411
1   244180
open 0.00 search 0.03 total 0.03
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.99 -q chembl_random_10.fpb --times -k all -c chembl_20.mfp2.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.99
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=chembl_20.mfp2.fpb
#query_sources=chembl_random_10.smi
#target_sources=/home/glandrum/Downloads/chembl_20.sdf.gz
1   552129
1   1109930
1   623138
1   1034319
1   766599
1   1067143
1   1116827
1   783512
1   357411
1   244180
open 0.00 search 0.02 total 0.02

Pulling the times out: doing the 10 similarity searches with tanimoto cutoffs of 0.4, 0.9, 0.95, and 0.99 took 120ms, 50ms, 30ms, and 20ms (the numbers are only reported with a precision of 10ms). As usual, chemfp is just stupid fast.

Update: Here's equivalent chemfp output for the zinc_clean set of 16.4 million fingerprints:

(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.4 -q chembl_random_10.fpb --times -k all -c zinc_all_clean.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.4
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=zinc_all_clean.fpb
#query_sources=chembl_random_10.smi
#target_sources=/scratch/RDKit_git/Data/Zinc/zinc_all_clean.smi.gz
92  552129
6675    1109930
116 623138
19890   1034319
5142    766599
679 1067143
131 1116827
1463    783512
0   357411
14  244180
open 0.00 search 1.31 total 1.31
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.9 -q chembl_random_10.fpb --times -k all -c zinc_all_clean.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.9
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=zinc_all_clean.fpb
#query_sources=chembl_random_10.smi
#target_sources=/scratch/RDKit_git/Data/Zinc/zinc_all_clean.smi.gz
0   552129
1   1109930
0   623138
4   1034319
2   766599
0   1067143
0   1116827
1   783512
0   357411
0   244180
open 0.00 search 0.58 total 0.58
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.95 -q chembl_random_10.fpb --times -k all -c zinc_all_clean.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.95
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=zinc_all_clean.fpb
#query_sources=chembl_random_10.smi
#target_sources=/scratch/RDKit_git/Data/Zinc/zinc_all_clean.smi.gz
0   552129
1   1109930
0   623138
2   1034319
2   766599
0   1067143
0   1116827
1   783512
0   357411
0   244180
open 0.00 search 0.41 total 0.41
(py27)glandrum@Otter:~/RDKit_blog/data$ simsearch -t 0.99 -q chembl_random_10.fpb --times -k all -c zinc_all_clean.fpb 
#Count/1
#num_bits=2048
#type=Count threshold=0.99
#software=chemfp/2.1
#queries=chembl_random_10.fpb
#targets=zinc_all_clean.fpb
#query_sources=chembl_random_10.smi
#target_sources=/scratch/RDKit_git/Data/Zinc/zinc_all_clean.smi.gz
0   552129
1   1109930
0   623138
2   1034319
2   766599
0   1067143
0   1116827
1   783512
0   357411
0   244180
open 0.00 search 0.22 total 0.22

Great stuff

Posted by at No comments:
Labels: , ,

Friday, May 29, 2015

Generating a reference dataset of chemical "scaffolds"

A set of scaffolds from CheMBL papers

Andrew Dalke has been working on putting together a dataset for testing substructure screenout performance for a while now. Due to concerns about privacy and/or confidentiality it's really difficult to get "real world" data on the kinds of things that scientists (really mostly chemists) look for in large-scale research systems like the data warehouses one finds in pharma companies. We have the data, but because of the nature of the searches we can't share it. Andrew visited us recently and the topic of substructure screening came up again. (Aside: There were three people in that room who cared deeply about substructure screening and who have invested some real thought and effort in it, how often does that happen? Sometimes I love my job!) During the conversation we came up with a way to fake having access to our actual queries. This post is about that.

One of the common use cases we see for substructure searching in the warehouse is to find all compounds that are in a particular chemical series. We obviously can't disclose the series that people are searching our database for, but we can fake it by collecting a bunch of substructure queries defining medchem-relevant chemical scaffolds and using those. This would have the added advantage of including things that have query features (variable atom and bond lists). Now we just need a source of medchem relevant scaffolds.

I did a post a while ago on the composition of the molecule sets we called "Datasets 2" in the model fusion paper. Since these are pulled from individual publications found in ChEMBL they tend to be made up of a set of molecules from one scaffold and then a few other molecules. One of the things I showed in that post was a method for estimating what the scaffold is in these sets using the RDKit's MCS code. This builds on that.

On to the work

Start with the preliminaries

We need to do a bunch of imports and define some functions we'll use later.

Note: This is a Python3 notebook, so the code below may not work in Python2.

In [1]:
import numpy
from rdkit.Chem.Draw import IPythonConsole
import pandas as pd
from rdkit.Chem import PandasTools
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import rdFMCS as MCS
import pickle
from rdkit import rdBase
print(rdBase.rdkitVersion)

import time
print(time.asctime())

2015.09.1pre
Thu May 28 04:42:14 2015
In [41]:
# define a function to convert a query molecule to SVG using the new rendering code.

from IPython.display import SVG
from rdkit.Chem import rdDepictor
from rdkit.Chem.Draw import rdMolDraw2D
def moltosvg(mol,molSize=(450,250),kekulize=True):
    mc = Chem.Mol(mol.ToBinary())
    if kekulize:
        try:
            Chem.Kekulize(mc)
        except:
            mc = Chem.Mol(mol.ToBinary())
    if not mc.GetNumConformers():
        rdDepictor.Compute2DCoords(mc)
    drawer = rdMolDraw2D.MolDraw2DSVG(molSize[0],molSize[1])
    # the MolDraw2D code is not very good at the moment at dealing with atom queries,
    # this is a workaround until that's fixed.
    # The rendering is still not going to be perfect because query bonds are not properly indicated
    opts = drawer.drawOptions()
    for atom in mc.GetAtoms():
        if atom.HasQuery() and atom.DescribeQuery().find('AtomAtomicNum')!=0:
            opts.atomLabels[atom.GetIdx()]=atom.GetSmarts()
    drawer.DrawMolecule(mc)
    drawer.FinishDrawing()
    svg = drawer.GetDrawingText()
    # It seems that the svg renderer used doesn't quite hit the spec.
    # Here are some fixes to make it work in the notebook, although I think
    # the underlying issue needs to be resolved at the generation step
    return svg.replace('svg:','')
In [43]:
# define a function to find the MCS for a set of molecules and return some summary info about that MCS

import time
# we will use a namedtuple to return the results
from collections import namedtuple
MCSRes=namedtuple('MCSRes',('smarts','numAtoms','numMols','avgNumMolAtoms','mcsTime'))

def MCS_Report(ms,printSmarts=True,
               atomCompare=MCS.AtomCompare.CompareAny,
               bondCompare=MCS.BondCompare.CompareAny,
               completeRingsOnly=True,
               **kwargs):
    t1=time.time()
    mcs = MCS.FindMCS(ms,atomCompare=atomCompare,bondCompare=bondCompare,timeout=30,completeRingsOnly=completeRingsOnly,
                      **kwargs)
    t2=time.time()
    nAts = numpy.array([x.GetNumAtoms() for x in ms])
    print('Mean nAts %.1f, mcs nAts: %d'%(nAts.mean(),mcs.numAtoms))
    if printSmarts: 
        print(mcs.smartsString)
    mcsM = Chem.MolFromSmarts(mcs.smartsString)         
    mcsM.UpdatePropertyCache(False)
    Chem.SetHybridization(mcsM)
    svg = moltosvg(mcsM,kekulize=False)
    tpl = MCSRes(mcs.smartsString,mcs.numAtoms,len(ms),nAts,t2-t1)
    return tpl,svg

Grab the benchmarking data

In [2]:
import glob
pkls = glob.glob('/home/glandrum/Code/benchmarking_platform/compounds/ChEMBL_II/*.pkl')
print(len(pkls))

37

Get target information from our local ChEMBL copy:

In [3]:
%load_ext sql
%config SqlMagic.feedback = False
In [4]:
import os
tgts={}
for pkl in pkls:
    fn = os.path.split(pkl)[-1]
    basen = os.path.splitext(fn)[0]
    tgtn=int(basen.split('_')[-1])
    data = %sql postgresql://localhost/chembl_20 \
          select tid,pref_name,organism from target_dictionary where tid=:tgtn ;
    tgts[data[0][0]]=(pkl,data[0][1],data[0][2])
In [5]:
for idx in sorted(tgts.keys()):
    pkl,tgt,species = tgts[idx]
    print(idx,tgt,species)

15 Carbonic anhydrase II Homo sapiens
25 Glucocorticoid receptor Homo sapiens
43 Beta-2 adrenergic receptor Homo sapiens
51 Serotonin 1a (5-HT1a) receptor Homo sapiens
61 Muscarinic acetylcholine receptor M1 Homo sapiens
65 Cytochrome P450 19A1 Homo sapiens
72 Dopamine D2 receptor Homo sapiens
87 Cannabinoid CB1 receptor Homo sapiens
90 Dopamine D4 receptor Homo sapiens
93 Acetylcholinesterase Homo sapiens
100 Norepinephrine transporter Homo sapiens
107 Serotonin 2a (5-HT2a) receptor Homo sapiens
108 Serotonin 2c (5-HT2c) receptor Homo sapiens
114 Adenosine A1 receptor Homo sapiens
121 Serotonin transporter Homo sapiens
126 Cyclooxygenase-2 Homo sapiens
130 Dopamine D3 receptor Homo sapiens
165 HERG Homo sapiens
259 Cannabinoid CB2 receptor Homo sapiens
10188 MAP kinase p38 alpha Homo sapiens
10193 Carbonic anhydrase I Homo sapiens
10260 Vanilloid receptor Homo sapiens
10280 Histamine H3 receptor Homo sapiens
10434 Tyrosine-protein kinase SRC Homo sapiens
10980 Vascular endothelial growth factor receptor 2 Homo sapiens
11140 Dipeptidyl peptidase IV Homo sapiens
11365 Cytochrome P450 2D6 Homo sapiens
11489 11-beta-hydroxysteroid dehydrogenase 1 Homo sapiens
11534 Cathepsin S Homo sapiens
11575 C-C chemokine receptor type 2 Homo sapiens
11631 Sphingosine 1-phosphate receptor Edg-1  Homo sapiens
12209 Carbonic anhydrase XII Homo sapiens
12252 Beta-secretase 1 Homo sapiens
12952 Carbonic anhydrase IX Homo sapiens
13001 Matrix metalloproteinase-2 Homo sapiens
17045 Cytochrome P450 3A4 Homo sapiens
19905 Melanin-concentrating hormone receptor 1 Homo sapiens

The data is organized in a dictionary with one entry per paper.

The MCS_Report() function defined above prints out the mean number of atoms per molecule in the input along with the size of the MCS. This is intended to help assess whether or not the MCS is actually a scaffold.

In [40]:
sets = pickle.load(open(tgts[11631][0],'rb'))
docid,cmpds = next(iter(sets.items()))
ids,smis= zip(*cmpds)
mols = [Chem.MolFromSmiles(x) for x in smis]

tpl,svg = MCS_Report(mols,threshold=0.8,completeRingsOnly=True)
SVG(svg)

Mean nAts 29.4, mcs nAts: 26
[#8,#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#6]1:[#7,#8]:[#6](-[#6]2:[#6]:[#6]:[#6](:[#6]:[#6]:2)-[#6]-[#7]2-[#6]-[#6](-[#6]-2)-[#6](=[#8])-[#8]):[#7,#8]:[#8,#7]:1
Out[40]:
[#8,#6] [#7,#8] N O OH [#7,#8] [#8,#7]

Find scaffolds in bulk

Now let's look at all the papers

In [38]:
def processTarget(tid):
    pkl,target,species=tgts[tid]
    print('-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*')
    print('\t',target,species)
    sets = pickle.load(open(pkl,'rb'))
    print('\t\t %d papers'%(len(sets)))
    alldata=[]
    for docid,cmpds in sets.items():
        data = %sql  \
              select journal,volume,first_page,year from docs where doc_id=:docid ;
        ref=', '.join(str(x) for x in data[0])
        print( '----------------')
        print(ref)
        ids,smis= zip(*cmpds)
        mols = [Chem.MolFromSmiles(x) for x in smis]
        tpl,svg = MCS_Report(mols,printSmarts=True,threshold=0.8,completeRingsOnly=True)
        alldata.append((tid,docid,ref,tpl,mols))
    return alldata
In [9]:
alldata=processTarget(11631)

-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  Sphingosine 1-phosphate receptor Edg-1  Homo sapiens
   7 papers
----------------
J. Med. Chem., 48, 6169, 2005
Mean nAts 29.4, mcs nAts: 26
[#8,#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#6]1:[#7,#8]:[#6](-[#6]2:[#6]:[#6]:[#6](:[#6]:[#6]:2)-[#6]-[#7]2-[#6]-[#6](-[#6]-2)-[#6](=[#8])-[#8]):[#7,#8]:[#8,#7]:1
----------------
Bioorg. Med. Chem., 15, 663, 2007
Mean nAts 24.9, mcs nAts: 16
[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#8,#6,#7]-[#6]-[#6](-[#7,#8])-[#6]-[#6,#8]-[#15](=[#8])(-[#8])-[#8]
----------------
Bioorg. Med. Chem. Lett., 22, 144, 2012
Mean nAts 30.6, mcs nAts: 29
[#6]-[#6]1=[#6](-[#6]-[#7]2-[#6]-[#6](-[#6]-2)-[#6](=[#8])-[#8])-[#6]-[#6]-[#6]2:[#6]-1:[#6]:[#6]:[#6](:[#6]:2)-[#8]-[#6]-[#6]-[#6]-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1
----------------
Bioorg. Med. Chem. Lett., 22, 1779, 2012
Mean nAts 31.2, mcs nAts: 28
[#7,#6]-[#6,#7]-[#6,#7]-[#6]1:[#6]:[#6]:[#6](:[#6](:[#6]:1)-[#9])-[#6]1:[#7]:[#6]2:[#6]:[#6]:[#6](:[#7]:[#6]:2:[#16]:1)-[#6]1(-[#6]2:[#6]:[#6]:[#6]:[#6]:[#6]:2)-[#6]-[#6]-1
----------------
Bioorg. Med. Chem. Lett., 14, 3351, 2004
Mean nAts 24.6, mcs nAts: 23
[#6]-[#6]-[#6]-[#6]-[#6]-[#6]-[#6]-[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#6]-[#6]-[#6](-[#7])-[#6]-[#6,#8]-[#6,#15,#16](=[#8])-,=[#8]
----------------
Bioorg. Med. Chem. Lett., 16, 3564, 2006
Mean nAts 31.4, mcs nAts: 14
[#6,#8]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#6]1-[#6]-[#6](-[#6]-,=[#6,#8])-[#6]-[#7]-1
----------------
Bioorg. Med. Chem. Lett., 21, 1390, 2011
Mean nAts 26.5, mcs nAts: 13
[#6]-[#7,#6]-[#6]-[#6]=,-[#6,#8]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#8,#6]-,=[#6]

Cool, do it for every target. Note that this takes a while; it would have been clever to do this using IPython's clustering, but I ended up just starting it and going off to do something else.

In [13]:
with open('results/papers_mcs_results.pkl','wb+') as outf:
    for tid in tgts:
        data = processTarget(tid)
        pickle.dump(data,outf)

-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  Cytochrome P450 19A1 Homo sapiens
   6 papers
----------------
J. Med. Chem., 47, 6792, 2004
Mean nAts 20.3, mcs nAts: 12
[#6]1:[#6]:[#6](-[#6]-[#7]2:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:2):[#6]:[#6,#8]:[#6,#8]:1
----------------
Bioorg. Med. Chem., 16, 8349, 2008
Mean nAts 26.1, mcs nAts: 16
[#6]-[#6]-[#7,#6](-[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#17,#6])-[#7]1:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:1
----------------
ACS Med. Chem. Lett., 2, 243, 2011
Mean nAts 30.2, mcs nAts: 21
[#7]#[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#7](-[#6]-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1)-[#7]1:[#6]:[#7]:[#7]:[#6]:1
----------------
J. Med. Chem., 50, 3420, 2007
Mean nAts 22.5, mcs nAts: 14
[#8]=[#6](-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1)-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1
----------------
Bioorg. Med. Chem., 16, 1474, 2008
Mean nAts 20.6, mcs nAts: 17
[#8]=[#6]1-,:[#6]-,:[#6](-[#6]2:[#6]:[#6]:[#6]:[#6]:[#6]:2)-,:[#8]-,:[#6]2:[#6]-,:1:[#6]:[#6]:[#6]:[#6]:2
----------------
J. Med. Chem., 53, 5347, 2010
Mean nAts 21.8, mcs nAts: 12
[#6]1:[#6]:[#7](-[#6]-[#6]2:[#6]:[#6]:[#6]:[#6]:[#6]:2):[#6]:[#7]:1
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  Tyrosine-protein kinase SRC Homo sapiens
   8 papers
----------------
Bioorg. Med. Chem. Lett., 16, 3993, 2006
Mean nAts 37.9, mcs nAts: 32
[#6]-[#8]-[#6]1:[#6]:[#6](-[#7]-[#6]2:[#6](-[#6]#[#7]):[#6]:[#7]:[#6]3:[#6]:2:[#6]:[#6](:[#6](:[#6]:3)-[#6]#[#6]-[#6]2:[#6,#7]:[#7,#6]:[#6]:[#6]:[#6,#7]:2)-[#8]-[#6]):[#6]:[#6]:[#6]:1-[#17,#8]
  
  ... SNIP ...

Read the data back in (this allows us to come back later and repeat analysis work without having to re-calculate all the scaffolds).

In [14]:
alldata = []
with open('results/papers_mcs_results.pkl','rb') as inf:
    while 1:
        try:
            row = pickle.load(inf)
        except EOFError:
            break
        alldata.append(row)
In [15]:
len(alldata)
Out[15]:
37
In [16]:
alldata[0]
Out[16]:
[(65,
  Decimal('35328'),
  'J. Med. Chem., 47, 6792, 2004',
  MCSRes(smarts='[#6]1:[#6]:[#6](-[#6]-[#7]2:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:2):[#6]:[#6,#8]:[#6,#8]:1', numAtoms=12, numMols=17, avgNumMolAtoms=array([24, 23, 18, 19, 17, 17, 19, 19, 24, 19, 19, 17, 25, 17, 19, 25, 24]), mcsTime=0.046170949935913086),
  [<rdkit.Chem.rdchem.Mol at 0x7f55246846d8>,
   <rdkit.Chem.rdchem.Mol at 0x7f55246847c8>,
   <rdkit.Chem.rdchem.Mol at 0x7f5524684688>,
   ...,
   <rdkit.Chem.rdchem.Mol at 0x7f5524613778>]),
 (65,
  Decimal('44596'),
  'Bioorg. Med. Chem., 16, 8349, 2008',
  MCSRes(smarts='[#6]-[#6]-[#7,#6](-[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#17,#6])-[#7]1:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:1', numAtoms=16, numMols=37, avgNumMolAtoms=array([27, 27, 23, 27, 30, 28, 28, 26, 27, 26, 27, 28, 23, 27, 28, 28, 25,
       20, 28, 30, 26, 28, 27, 27, 27, 17, 27, 26, 26, 27, 26, 22, 20, 27,
       28, 23, 28]), mcsTime=0.4323115348815918),
   
   ... SNIP ...

Create a Pandas DataFrame with the data so that we can explore it a bit.

Note The molecule renderings in this table are not optimal. I was unable to figure out how to get the SVG to render properly in the notebook. Still, this is hopefully better than nothing.

In [72]:
tbl = []
for row in alldata:
    for r in row:
        tid,docid,ref,tpl,_ = r
        m = Chem.MolFromSmarts(tpl.smarts)
        m.UpdatePropertyCache(False)
        #svg = moltosvg(m,kekulize=False)
        tbl.append((tid,tpl.numAtoms,tpl.avgNumMolAtoms.mean(),m,tpl.smarts,tpl.mcsTime))
df = pd.DataFrame(data=tbl,columns=('tid','numAtoms','avgNumMolAtoms','mol','smarts','time'))
PandasTools.RenderImagesInAllDataFrames(images=True)
df.head(5)
Out[72]:
tid numAtoms avgNumMolAtoms mol smarts time
0 65 12 20.294118 Mol [#6]1:[#6]:[#6](-[#6]-[#7]2:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:2):[#6]:[#6,#8]:[#6,#8]:1 0.046171
1 65 16 26.081081 Mol [#6]-[#6]-[#7,#6](-[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#17,#6])-[#7]1:[#6]:[#6,#7]:[#7,#6]:[#6,#7]:1 0.432312
2 65 21 30.214286 Mol [#7]#[#6]-[#6]1:[#6]:[#6]:[#6](:[#6]:[#6]:1)-[#7](-[#6]-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1)-[#7]1:[#6]:[#7]:[#7]:[#6]:1 0.608892
3 65 14 22.500000 Mol [#8]=[#6](-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1)-[#6]1:[#6]:[#6]:[#6]:[#6]:[#6]:1 0.044816
4 65 17 20.583333 Mol [#8]=[#6]1-,:[#6]-,:[#6](-[#6]2:[#6]:[#6]:[#6]:[#6]:[#6]:2)-,:[#8]-,:[#6]2:[#6]-,:1:[#6]:[#6]:[#6]:[#6]:2 0.160467

How many "scaffolds" do we have in total?

In [18]:
len(df)
Out[18]:
648

How many are reasonably sized?

In [19]:
len(df[df.numAtoms>6])
Out[19]:
640

Not bad! It's still worth looking at the ones that are super small:

In [20]:
df[df.numAtoms<=6]
Out[20]:
tid numAtoms avgNumMolAtoms mol smarts time
128 13001 6 33.538462 Mol [#6]1:[#6]:[#6]:[#6]:[#6]:[#6,#7]:1 0.010566
188 15 5 11.400000 Mol [#6,#7,#8]-[#6,#7,#8]-[#6](-,=[#6,#8])-,=[#6,#8] 0.001489
189 15 6 16.300000 Mol [#7,#6,#8,#16]-,=[#6]1:,-[#7,#16]:,-[#7,#6,#8,#16]:,-[#6,#7]:[#16,#6,#7]:,-1 0.008936
200 15 5 23.608696 Mol [#6,#7,#8]-,=[#6,#16](=,-[#8,#6])-[#7,#6]-,=[#6] 0.005014
206 15 6 21.000000 Mol [#6]1:[#6]:[#6]:[#6]:[#6,#7]:[#6]:1 0.015131
223 10193 6 15.791667 Mol [#6]1:,-[#6,#7]:,-[#6,#8]:,-[#6]:,-[#6]:,-[#6,#7]:,-1 0.005510
238 10193 5 11.900000 Mol [#7,#6]-[#6](-[#6,#7]-,=[#6,#8])-,=[#6,#7,#8] 0.000738
249 10193 5 11.400000 Mol [#6,#7,#8]-[#6,#7,#8]-[#6](-,=[#6,#8])-,=[#6,#8] 0.001440

We can also render those better:

In [68]:
SVG(moltosvg(df[df.numAtoms<=6].iloc[0,3],kekulize=False,molSize=(250,150)))
Out[68]:
[#6,#7]
In [69]:
SVG(moltosvg(df[df.numAtoms<=6].iloc[1,3],kekulize=False,molSize=(250,150)))
Out[69]:
[#6,#7,#8] [#6,#7,#8] [#6,#8] [#6,#8]
In [70]:
SVG(moltosvg(df[df.numAtoms<=6].iloc[2,3],kekulize=False,molSize=(250,150)))
Out[70]:
[#7,#6,#8,#16] [#7,#16] [#7,#6,#8,#16] [#6,#7] [#16,#6,#7]
In [71]:
SVG(moltosvg(df[df.numAtoms<=6].iloc[3,3],kekulize=False,molSize=(250,150)))
Out[71]:
[#6,#7,#8] [#6,#16] [#8,#6] [#7,#6]

Wrapping up

There's a lot more we can do with this set of scaffolds. There are more posts to come...

Posted by at No comments:
Labels: , ,
Subscribe to: Posts (Atom)