Sunday, August 10, 2014

Picking diverse compounds from large sets

Draft post This post will probably be edited over the next day or so.

This post explores two different approaches for picking diverse small subsets of compounds from large sets: random selection and the MaxMin algorithm (Ashton, M. et. al., Quant. Struct.-Act. Relat., 21 (2002), 598-604).

In [1]:

import numpy as np
from __future__ import print_function

from rdkit import Chem
from rdkit.Chem import Draw,rdMolDescriptors,AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit import SimDivFilters,DataStructs

import gzip
%pylab inline

from rdkit import rdBase
print(rdBase.rdkitVersion)
import time
print(time.asctime())

Populating the interactive namespace from numpy and matplotlib
2014.09.1pre
Mon Aug 11 07:46:10 2014

Start by reading in a set of molecules from the Zinc Biogenic Compounds (formerly known as "Zinc Natural Products") subset and generating Morgan2 fingerprints for the diversity calculation. The compound structures are available here: http://zinc.docking.org/subsets/zbc

In [2]:

fps = []
inf = gzip.open('../data/znp.sdf.gz')
suppl = Chem.ForwardSDMolSupplier(inf)
for m in suppl:
    if m is None:
        continue
    fp=rdMolDescriptors.GetMorganFingerprintAsBitVect(m,2)
    fps.append(fp)
    if not len(fps)%5000: print("Done: %d"%len(fps))

print("Got %d fingerprints"%len(fps))

Done: 5000
Done: 10000
Done: 15000
Done: 20000
Done: 25000
Done: 30000
Done: 35000
Done: 40000
Done: 45000
Done: 50000
Done: 55000
Done: 60000
Done: 65000
Done: 70000
Done: 75000
Done: 80000
Done: 85000
Done: 90000
Done: 95000
Done: 100000
Done: 105000
Done: 110000
Done: 115000
Done: 120000
Done: 125000
Done: 130000
Done: 135000
Done: 140000
Got 144069 fingerprints

In [5]:

import time
mmp =SimDivFilters.MaxMinPicker()
t1=time.time()
bv_ids = mmp.LazyBitVectorPick(fps,len(fps),50)
t2=time.time()
print("That took %.2f seconds"%(t2-t1))

That took 43.99 seconds

In [10]:

dist_hist=[]
for i in range(len(bv_ids)):
    for j in range(i+1,len(bv_ids)):
        dist_hist.append(DataStructs.TanimotoSimilarity(fps[bv_ids[i]],fps[bv_ids[j]]))
_=hist(dist_hist,bins=20)
_=title("MaxMin Picks")
_=xlabel("Similarity")

Those are certainly diverse

What about simply picking random compounds?

In [11]:

import random
random.seed(0xF00D)
rr = list(range(len(fps)))
random.shuffle(rr)
r_ids = rr[:50]
dist_hist2=[]
for i in range(len(r_ids)):
    for j in range(i+1,len(r_ids)):
        dist_hist2.append(DataStructs.TanimotoSimilarity(fps[r_ids[i]],fps[r_ids[j]]))
_=hist(dist_hist2,bins=20)
_=title("Random Picks")
_=xlabel("Similarity")

At first glance, it looks like the random picks are much less diverse than the picks with the MaxMin algorithm: the peak in the first histogram is significantly shifted to higher similarity.

However: there is an important detail to keep in mind: the "noise level" of the fingerprint. In an earlier post (http://rdkit.blogspot.ch/2013/10/fingerprint-thresholds.html), I looked at similarities between random pairs of "drug-like" compounds from ChEMBL in order to identify thresholds for when a similarity value is at the noise level (defined by finding the value which 99% of random pairs have a lower similarity). For the fingerprint used here, Morgan2, this value is 0.321.

Using that cutoff, we can easily figure out what fraction of the pairs in the set above are "noise":

In [12]:

noise=0.321
count = sum([1 for x in dist_hist2 if x<noise])
print("Count: %d, fraction: %.3f"%(count,1.*count/len(dist_hist2)))

Count: 1210, fraction: 0.988

So almost 99% of the pairs are below the noise level.

Verify that this wasn't just lucky by doing some stats:

In [15]:

noise=0.321
counts=[]
for rep in range(100):
    random.seed(0xF00D+rep)
    rr = list(range(len(fps)))
    random.shuffle(rr)
    r_ids = rr[:50]
    count=0
    nTot=0
    for i in range(len(r_ids)):
        for j in range(i+1,len(r_ids)):
            d=DataStructs.TanimotoSimilarity(fps[r_ids[i]],fps[r_ids[j]])
            nTot+=1
            if d<noise:
                count+=1
    counts.append(count)
counts = numpy.array(counts,numpy.double)
counts /= nTot
mean = numpy.mean(counts)
dev = numpy.std(counts)
print("Mean: %.3f, std dev: %.3f"%(mean,dev))

Mean: 0.991, std dev: 0.004

Fazit: randomly picking small numbers of compounds from the large set gives subsets that are quite diverse. It's also a lot quicker than using the diversity picker: not only can we skip the diversity picking algorithm, but the molecules and fingerprints don't even need to be constructed.

Clearly this is not going to work as well when picking from compounds that are reasonably closely related to each other, but for small subsets of large diverse sets, it's a solid approach.

Optimizing Diversity Picking in the RDKit

The RDKit has an implementation of the MaxMin algorithm (Ashton, M. et. al., Quant. Struct.-Act. Relat., 21 (2002), 598-604) for diversity picking. This provides a computationally efficient method for picking diverse subsets of molecules from a larger set.

This post explores the performance of the different ways of invoking the algorithm and does a bit of looking at the scaling of the runtime with set size.

In [1]:

import numpy as np
from __future__ import print_function

from rdkit import Chem
from rdkit.Chem import Draw,rdMolDescriptors,AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit import SimDivFilters,DataStructs

import gzip
%pylab inline

from rdkit import rdBase
print(rdBase.rdkitVersion)
import time
print(time.asctime())

Populating the interactive namespace from numpy and matplotlib
2014.09.1pre
Mon Aug 11 07:51:19 2014

In [10]:

ms = []
inf = gzip.open('../data/znp.sdf.gz')
suppl = Chem.ForwardSDMolSupplier(inf)
while len(ms)<20000:
    m = next(suppl)
    AllChem.Compute2DCoords(m)
    if m is not None: ms.append(m)

In [4]:

fps = [rdMolDescriptors.GetMorganFingerprintAsBitVect(m,2) for m in ms]

There are two basic ways to use the MaxMinPicker

Generate the distance matrix in advance and do the picking based on those pre-calculated distances
Use a "lazy" form of the picker that generates distances using a caller-supplied function. For reasons of efficiency, this uses a cache internally so that distances don't have to be calculated more than once.

Here's a demonstration of using the first approach to pick 50 diverse instances from the first 1000 compounds in the set:

In [5]:

def dmat_sim(fps,ntopick):
    ds=[]
    for i in range(1,len(fps)):
         ds.extend(DataStructs.BulkTanimotoSimilarity(fps[i],fps[:i],returnDistance=True))
    mmp =SimDivFilters.MaxMinPicker()
    ids=mmp.Pick(numpy.array(ds),len(fps),ntopick)
    return ids

dmat_ids=dmat_sim(fps[:1000],50)

Note that there are some examples of using this approach floating around on the web that calculate the distance matrix in the wrong order. I've done my best to find and either fix or remove them, but there are no doubt still some bad ones out there.

In [7]:

print(list(dmat_ids)[:10])

[374, 636, 790, 497, 724, 366, 418, 587, 433, 537]

In [11]:

Draw.MolsToGridImage([ms[x] for x in dmat_ids[:10]],molsPerRow=4)

Out[11]:

Those are certainly diverse

Try the same thing using the lazy picker and a user-supplied function

In [13]:

def fn(i,j,fps=fps):
    return 1.-DataStructs.TanimotoSimilarity(fps[i],fps[j])
mmp =SimDivFilters.MaxMinPicker()
lazy_ids = mmp.LazyPick(fn,1000,50)

Make sure the results are the same

In [14]:

list(lazy_ids)==list(dmat_ids)

Out[14]:

True

When working with bit vectors, a third option that can be used that includes the distance calculation directly in the C++ call. I'll show below that this makes a difference in performance.

Note the LazyBitVectorPick call is a new addition to the RDKit. At the time of this writing it is only available in the github version of the code. It will appear in the 2014.09 release.

In [15]:

mmp =SimDivFilters.MaxMinPicker()
bv_ids = mmp.LazyBitVectorPick(fps[:1000],1000,50)

In [16]:

list(bv_ids)==list(dmat_ids)

Out[16]:

True

Now let's look at the relative performance of these approaches for a subset of the data.

Timing data was generated on a not-particularly-modern machine: a four-year old MacBook Pro

In [18]:

mmp =SimDivFilters.MaxMinPicker()
%timeit dmat_sim(fps[:2000],50)

1 loops, best of 3: 1.77 s per loop

In [19]:

mmp =SimDivFilters.MaxMinPicker()
%timeit mmp.LazyPick(fn,2000,50)

1 loops, best of 3: 1.01 s per loop

In [23]:

mmp =SimDivFilters.MaxMinPicker()
%timeit mmp.LazyBitVectorPick(fps,2000,50)

1 loops, best of 3: 458 ms per loop

So we have a clear winner: the LazyBitVector picker is considerably faster than the other alternatives.

Impact of the cache

The LazyBitVectorPicker is very runtime efficient, but the memory usage of the cache that it uses (about nPoints * nPicks * 30bytes) could theoretically cause problems for large data sets. The cache can be disabled, at the expense of substantially longer run time:

In [22]:

mmp =SimDivFilters.MaxMinPicker()
%timeit mmp.LazyBitVectorPick(fps,2000,50,useCache=False)

1 loops, best of 3: 1.37 s per loop

Impact of data set size

Let's look at the evolution of the run time with the number of compounds we're picking from

In [36]:

for sz in range(1,11):
    sz*=1000
    print("Doing: %1d"%(sz))
    mmp =SimDivFilters.MaxMinPicker()
    %timeit mmp.LazyBitVectorPick(fps,sz,50)

Doing: 1000
10 loops, best of 3: 186 ms per loop
Doing: 2000
1 loops, best of 3: 417 ms per loop
Doing: 3000
1 loops, best of 3: 688 ms per loop
Doing: 4000
1 loops, best of 3: 955 ms per loop
Doing: 5000
1 loops, best of 3: 1.21 s per loop
Doing: 6000
1 loops, best of 3: 1.49 s per loop
Doing: 7000
1 loops, best of 3: 1.8 s per loop
Doing: 8000
1 loops, best of 3: 2.34 s per loop
Doing: 9000
1 loops, best of 3: 2.22 s per loop
Doing: 10000
1 loops, best of 3: 2.55 s per loop

The runtime increases more or less linearly with the size of the set being picked from

Impact of the number of points being picked

In [37]:

for sz in range(1,11):
    sz*=20
    print("Doing: %1d"%(sz))
    mmp =SimDivFilters.MaxMinPicker()
    %timeit mmp.LazyBitVectorPick(fps,5000,sz)

Doing: 20
1 loops, best of 3: 201 ms per loop
Doing: 40
1 loops, best of 3: 737 ms per loop
Doing: 60
1 loops, best of 3: 1.69 s per loop
Doing: 80
1 loops, best of 3: 3.03 s per loop
Doing: 100
1 loops, best of 3: 4.78 s per loop
Doing: 120
1 loops, best of 3: 6.91 s per loop
Doing: 140
1 loops, best of 3: 9.64 s per loop
Doing: 160
1 loops, best of 3: 12.4 s per loop
Doing: 180
1 loops, best of 3: 15.8 s per loop
Doing: 200
1 loops, best of 3: 20 s per loop

That's more quadratic in the number of picks. Given the algorithm, which needs to compare each potential point to pick against all the points that have so-far been picked, this makes perfect sense.