Fingerprint-based substructure screening 1
[Edited 10/11/2013 due to the bug fix for https://github.com/rdkit/rdkit/issues/151]Some form of fingerprint is often used to make substructure searching more efficient. The idea is to use a fingerprinting algorithm with the property:
FP(query) & FP(mol) = FP(query) if query is a substructure of mol. In words: if
query is a substructure of mol then every bit set in the fingerprint of query should also be set in mol.A bunch of different approaches to this have been developed, I'm not going to review them all here. :-)
Andrew Dalke has done some writing on the topic: http://dalkescientific.com/writings/diary/archive/2012/06/11/optimizing_substructure_keys.html
One of the best-known approaches is the Daylight algorithm: http://www.daylight.com/dayhtml/doc/theory/theory.finger.html
As Andrew mentions in the post I link to above, one of the big problems here is the lack of a reasonable query dataset for benchmarking. The approach I've taken is to use collections of small molecules from ZINC as well as some fragments of PubChem molecules. The database I query is a set of ChEMBL molecules from an earlier post.
Here's the atom-number count distributions of the queries and molecules:
It's not perfect, but it's a start.
TL;DR summary
Using the RDKit pattern fingerprinter (the default used in the postgresql cartridge, more on that in another post), a screenout accuracy of around 60% is achievable with these three query sets. This means that 60% of molecules passing the fingerprint screen actually have a substructure match.To put that in perspective: with the leads query set and the pattern fingerprint, 4650 substructure matches are found in a total of 7601 searches, that's 61% accuracy for the fingerprint. If the fingerprint had not been used to pre-screen, 25000000 (500*50000) searches would have needed to be done, so we've reduced the number of substructure searches by 99.97%. By this logic, even a comparatively poor fingerprint with an accuracy of 5% would help enormously: reducing the number of searches by 99.7%. This will become important when actually using the fingerprints in a database index, but that's a different post.
On the technical side: there's some info below about using IPython's ability to run code in parallel.
In [1]:
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors
from rdkit.Avalon import pyAvalonTools
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
from rdkit import rdBase
from rdkit import DataStructs
import cPickle,random,gzip,time
from __future__ import print_function
print(rdBase.rdkitVersion)
Read in the data
Database molecules
Here we will use the 50K molecules that make up the set of 25K reference pairs generated in an earlier post: http://rdkit.blogspot.ch/2013/10/building-similarity-comparison-set-goal.htmlAs a quick reminder: these are pairs of molecules taken from ChEMBL with MW<600 and a count-based MFP0 similarity of at least 0.7 to each other.
In [2]:
ind = [x.split() for x in gzip.open('../data/chembl16_25K.pairs.txt.gz')]
mols = []
for i,row in enumerate(ind):
m1 = Chem.MolFromSmiles(row[1])
mols.append(m1)
m2 = Chem.MolFromSmiles(row[3])
mols.append(m2)
Query molecules
We'll use three sets:- Fragments: 500 diverse molecules taken from the ZINC Fragments set
- Leads: 500 diverse molecules taken from the ZINC Lead-like set
- Pieces: 823 pieces of molecules obtained by doing a BRICS fragmentation of some molecules from the pubchem screening set.
In [3]:
frags = [Chem.MolFromSmiles(x.split()[0]) for x in file('../data/zinc.frags.500.q.smi')]
leads = [Chem.MolFromSmiles(x.split()[0]) for x in file('../data/zinc.leads.500.q.smi')]
pieces = [Chem.MolFromSmiles(x) for x in file('../data/fragqueries.q.txt')]
Look at the size distributions:
In [97]:
figure(figsize=(16,4))
subplot('141')
hist([x.GetNumAtoms() for x in pieces],bins=20)
title('Pieces')
xlabel('NumAtoms')
subplot('142')
hist([x.GetNumAtoms() for x in frags],bins=20)
title('Fragments')
xlabel('NumAtoms')
subplot('143')
hist([x.GetNumAtoms() for x in leads],bins=20)
title('Leads')
_=xlabel('NumAtoms')
subplot('144')
hist([x.GetNumAtoms() for x in mols],bins=20)
title('Mols')
_=xlabel('NumAtoms')
Setting up the parallelization
I'm going to be doing a fair amount of embarassingly parallel computation, so I'll take advantage of IPython's support for parallel computing. For this to work, you need to have a cluster of workers set up. This is easy in the Notebook (there's a tab in the dashboard), and it's doable from the command line (http://ipython.org/ipython-doc/stable/parallel/parallel_process.html)It took a bit of time to figure out how to do this, but once over that hurdle, this thing is super easy and useful.
In [4]:
from IPython.parallel import Client
rc = Client()
dview = rc[:]
In [5]:
with dview.sync_imports():
from rdkit import Chem
from rdkit import DataStructs
from rdkit.Avalon import pyAvalonTools
Verify that it works:
In [28]:
t1=time.time()
_ = [Chem.RDKFingerprint(m) for m in mols[:10000]]
t2=time.time()
print('Serial: %.2f'%(t2-t1))
In [70]:
fn = lambda x:Chem.RDKFingerprint(x)
t1=time.time()
_ = dview.map_sync(fn,mols[:10000])
t2=time.time()
print('Parallel: %.2f'%(t2-t1))
The above runtimes were collected with a cluster of 4 workers, so it looks like the parallization is working.
The test harness
In [6]:
def subTest(mol,queries,qfps,fpf,verify,silent):
nFound=0
nScreened=0
nTot=0
nFailed = 0
molfp = fpf(mol,False)
for j,queryfp in enumerate(qfps):
nTot += 1
if DataStructs.AllProbeBitsMatch(queryfp,molfp):
nScreened += 1
if mol.HasSubstructMatch(queries[j]):
nFound += 1
elif verify:
if mol.HasSubstructMatch(queries[j]):
nFailed += 1
if not silent:
print('Failure: %s %s'%(Chem.MolToSmiles(mol,True),Chem.MolToSmiles(queries[j],True)))
return nFound,nScreened,nTot,nFailed
def testScreenout(mols,queries,fpf,dview,verify=False,silent=False):
nFound=0
nScreened=0
nTot=0
nFailed = 0
#if not silent: print("Building query fingerprints")
qfps = [fpf(m,True) for m in queries]
#if not silent: print("Running Queries")
t = lambda x,qfps=qfps,queries=queries,fpf=fpf,verify=verify,silent=silent,subTest=subTest:subTest(x,queries,qfps,fpf,verify,silent)
res=dview.map_async(t,mols)
for entry in res:
nFound+=entry[0]
nScreened+=entry[1]
nTot+=entry[2]
nFailed+=entry[3]
if not silent:
accuracy = float(nFound)/nScreened
print("Found %d matches in %d searches with %d failures. Accuracy: %.3f"%(nFound,nScreened,nFailed,accuracy))
return nTot,nScreened,nFound,nFailed
Tests
In [11]:
methods = {
'Avalon-2K':lambda x,y: pyAvalonTools.GetAvalonFP(x,nBits=2048,isQuery=y,bitFlags=pyAvalonTools.avalonSSSBits),
'Pattern-1K':lambda x,y: Chem.PatternFingerprint(x,fpSize=1024),
'Pattern-2K':lambda x,y: Chem.PatternFingerprint(x,fpSize=2048),
'Pattern-4K':lambda x,y: Chem.PatternFingerprint(x,fpSize=4096),
'Layered': lambda x,y:Chem.LayeredFingerprint(x,layerFlags=Chem.LayeredFingerprint_substructLayers),
}
In [12]:
leadResults={}
fragResults={}
pieceResults={}
for method,func in methods.iteritems():
print("----------------------------")
print("Doing %s"%method)
print("Leads")
leadResults[method]=testScreenout(mols,leads,func,dview)
print("Frags")
fragResults[method]=testScreenout(mols,frags,func,dview)
print("Pieces")
pieceResults[method]=testScreenout(mols,pieces,func,dview)
Summarize that:
In [14]:
mns = sorted(methods.keys())
for mn in mns:
print(mn,'Pieces','%.3f'%(float(pieceResults[mn][2])/pieceResults[mn][1]),'Fragments','%.3f'%(float(fragResults[mn][2])/fragResults[mn][1]),
'Leads','%.3f'%(float(leadResults[mn][2])/leadResults[mn][1]))
A caveat
The results from the Avalon fingerprint shown above are not directly comparable to the others since structures are being filtered out that shouldn't be.To demonstrate this, we need a serial form of the screenout test:
In [8]:
def testScreenoutSerial(mols,queries,fpf,verify=False,silent=False):
nFound=0
nScreened=0
nTot=0
nFailed = 0
#if not silent: print("Building query fingerprints")
qfps = [fpf(m,True) for m in queries]
#if not silent: print("Running Queries")
t = lambda x,qfps=qfps,queries=queries,fpf=fpf,verify=verify,silent=silent,subTest=subTest:subTest(x,queries,qfps,fpf,verify,silent)
res=map(t,mols)
for entry in res:
nFound+=entry[0]
nScreened+=entry[1]
nTot+=entry[2]
nFailed+=entry[3]
if not silent:
accuracy = float(nFound)/nScreened
print("Found %d matches in %d searches with %d failures. Accuracy: %.3f"%(nFound,nScreened,nFailed,accuracy))
return nTot,nScreened,nFound,nFailed
And now we can get some examples:
In [16]:
testScreenoutSerial(mols[:500],leads,methods['Avalon-2K'],verify=True)
Out[16]:
Look at a specific example:
In [17]:
mfp = pyAvalonTools.GetAvalonFP('Cc1ccc(C)c(S(=O)(=O)c2nnn3c4ccsc4c(Nc4ccc(C)c(C)c4)nc23)c1',True,2048,isQuery=False)
qfp = pyAvalonTools.GetAvalonFP('c1cn[nH]n1',True,2048,isQuery=True)
DataStructs.AllProbeBitsMatch(qfp,mfp)
Out[17]:
In [18]:
Chem.MolFromSmiles('Cc1ccc(C)c(S(=O)(=O)c2nnn3c4ccsc4c(Nc4ccc(C)c(C)c4)nc23)c1')
Out[18]:
The normal explanation for this would be differences in the aromaticity model. In this case, however, that turns out not be exactly it.
Here's the molecule we're querying with:
Here's the molecule we're querying with:
In [19]:
Chem.MolFromSmiles('c1cn[nH]n1')
Out[19]:
In [20]:
qfp2 = pyAvalonTools.GetAvalonFP('c1cnn[nH]1',True,2048,isQuery=True)
qfp==qfp2
Out[20]:
The deciding factor is the tautomer. The Avalon fingerprinter considers this significant, and that difference explains the different screenout.
If we use the fingerprint from the second tautomer we see the expected screenout behavior:
If we use the fingerprint from the second tautomer we see the expected screenout behavior:
In [102]:
DataStructs.AllProbeBitsMatch(qfp2,mfp)
Out[102]:
This ends up mattering because the RDKit's substructure matcher ignores the H completely:
In [103]:
Chem.MolFromSmiles('c1cnn[nH]1').HasSubstructMatch(Chem.MolFromSmiles('c1cn[nH]n1'))
Out[103]:
In [104]:
Chem.MolFromSmiles('c1cnnn1C').HasSubstructMatch(Chem.MolFromSmiles('c1cn[nH]n1'))
Out[104]:
Unless, of course, you construct the query from SMARTS:
In [105]:
Chem.MolFromSmiles('c1cnnn1C').HasSubstructMatch(Chem.MolFromSmarts('c1cn[nH]n1'))
Out[105]:
Other fingerprints
The RDKit fingerprint
In [27]:
fn = lambda x,v:Chem.RDKFingerprint(x,minPath=1,maxPath=5,nBitsPerHash=1,useHs=False)
testScreenoutSerial(mols,leads,fn,verify=True)
Out[27]:
The failures happen here because the RDKit fingerprinter uses information about aromaticity in the calculation of atom invariants and while hashing bonds. That extra information is a big part of why the accuracy is so high.
Turning off all bond order information (finer grained control is unfortunately not possible) and using atomic number as the atom invariants gives no errors, but a much lower screenout rate:
Turning off all bond order information (finer grained control is unfortunately not possible) and using atomic number as the atom invariants gives no errors, but a much lower screenout rate:
In [53]:
def func(mol,v):
invs = [x.GetAtomicNum() for x in mol.GetAtoms()]
return Chem.RDKFingerprint(mol,minPath=1,maxPath=5,nBitsPerHash=1,useHs=False,useBondOrder=False,atomInvariants=invs)
testScreenoutSerial(mols,leads,func,verify=True)
Out[53]:
Using a shorter maxPath length also hurts accuracy:
In [35]:
fn = lambda x,v:Chem.RDKFingerprint(x,minPath=1,maxPath=4,nBitsPerHash=1,useHs=False)
testScreenout(mols,leads,fn,dview)
Out[35]:
It is most likely luck that there are no failures here.
Chemfp
Andrew Dalke's chemfp package (http://chemfp.com/) includes a pattern-based fingerprinter which is based on the PubChem/CACTVS substructure keys:
In [31]:
from chemfp import rdkit_patterns
dalke_fp=rdkit_patterns.SubstructRDKitFingerprinter_v1.make_fingerprinter()
fn = lambda x,v,fp=dalke_fp:DataStructs.CreateFromBinaryText(fp(x))
testScreenoutSerial(mols[:500],leads,fn,verify=True)
Out[31]:
It's not particularly effective as a screening fingerprint and, as seen above, creates some holes.
MACCS keys
Another key-based fingerprint that isn't particularly effective and results in some misses:
In [9]:
fn = lambda x,v:rdMolDescriptors.GetMACCSKeysFingerprint(x)
testScreenoutSerial(mols[:500],leads,fn,verify=True)
Out[9]:
Improving things
A first pass at improving is to combine some selected patterns with an algorithmic fingerprinter like the pattern fingerprinter. Here's a thread describing an experiment:http://www.mail-archive.com/rdkit-discuss@lists.sourceforge.net/msg02078.html
And a reproduction of that with this test harness:
In [15]:
qs="""0 2 times O largest 55458
1 2 times Ccc largest 29602
2 1 times CCN largest 16829
3 1 times cnc largest 11439
4 1 times cN largest 8998
5 1 times C=O largest 7358
6 1 times CCC largest 6250
7 1 times S largest 4760
8 1 times c1ccccc1 largest 4524
9 2 times N largest 2854
10 1 times C=C largest 2162
11 1 times nn largest 1840
12 2 times CO largest 1248
13 1 times Ccn largest 964
14 1 times CCCCC largest 857
15 1 times cc(c)c largest 653
16 3 times O largest 653
17 1 times O largest 466
18 2 times CNC largest 464
19 1 times s largest 457
20 1 times CC(C)C largest 335
21 1 times o largest 334
22 1 times cncnc largest 334
23 1 times C=N largest 321
24 2 times CC=O largest 238
25 4 times Ccc largest 238
26 1 times Cl largest 230
27 4 times O largest 149
28 2 times ccncc largest 149
29 6 times CCCCCC largest 76
30 2 times c1ccccc1 largest 76
31 1 times F largest 75
32 3 times CCOC largest 44
33 3 times N largest 44
34 1 times c(cn)n largest 44
35 1 times N largest 41
36 9 times C largest 41
37 1 times CC=C(C)C largest 33
38 1 times c1ccncc1 largest 26
39 1 times CC(C)N largest 26
40 1 times CC largest 26
41 4 times CCC(C)O largest 25
42 2 times ccc(cc)n largest 21
43 6 times C largest 21
44 1 times C1CCCC1 largest 18
45 1 times C largest 18
46 5 times O largest 18
47 2 times Ccn largest 14
48 1 times CNCN largest 13
49 3 times cncn largest 13
50 1 times CSC largest 13
51 3 times CC=O largest 11
52 1 times CCNCCCN largest 11
53 1 times CccC largest 11
54 3 times ccccc(c)c largest 10"""
def _initPatts():
ssPatts=[]
for q in qs.split('\n'):
q = q.split(' ')
count = int(q[1])
patt = Chem.MolFromSmiles(q[3],sanitize=False)
patt.UpdatePropertyCache(strict=False) # <- github #149, without this we cannot pickle the queries
ssPatts.append((patt,count))
return ssPatts
_ssPatts=None
def GetCombinedSubstructFP(m,ssPatts=None,fpSize=1024,verbose=False):
if ssPatts is None:
ssPatts = _initPatts()
sz = len(ssPatts)
lfp=Chem.PatternFingerprint(m,fpSize=fpSize)
res = DataStructs.ExplicitBitVect(fpSize+sz)
obls = [x+sz for x in lfp.GetOnBits()]
res.SetBitsFromList(obls)
for i,(p,count) in enumerate(ssPatts):
matches = m.GetSubstructMatches(p,uniquify=True)
if len(matches)>=count:
res.SetBit(i)
return res
In [16]:
ssPatts = _initPatts()
In [65]:
fn = lambda x,v,ssPatts=ssPatts:GetCombinedSubstructFP(x,ssPatts=ssPatts,fpSize=2048)
testScreenoutSerial(mols[:500],leads,fn,verify=True)
Out[65]:
In [17]:
leadResults={}
fragResults={}
pieceResults={}
print("----------------------------")
fn = lambda x,v,ssPatts=ssPatts,fpf=GetCombinedSubstructFP:fpf(x,ssPatts=ssPatts,fpSize=1024)
method='1024'
print("Doing %s"%method)
print("Pieces")
pieceResults[method]=testScreenout(mols,pieces,fn,dview)
print("Frags")
fragResults[method]=testScreenout(mols,frags,fn,dview)
print("Leads")
leadResults[method]=testScreenout(mols,leads,fn,dview)
print("----------------------------")
fn = lambda x,v,ssPatts=ssPatts,fpf=GetCombinedSubstructFP:fpf(x,ssPatts=ssPatts,fpSize=2048)
method='2048'
print("Doing %s"%method)
print("Pieces")
pieceResults[method]=testScreenout(mols,pieces,fn,dview)
print("Frags")
fragResults[method]=testScreenout(mols,frags,fn,dview)
print("Leads")
leadResults[method]=testScreenout(mols,leads,fn,dview)
Comparing to the previous pattern fp results:
Pattern-1K Pieces 0.500 Fragments 0.285 Leads 0.169
Combine-1K Pieces 0.640 Fragments 0.338 Leads 0.186
Pattern-2K Pieces 0.572 Fragments 0.590 Leads 0.715
Combine-2K Pieces 0.689 Fragments 0.608 Leads 0.723
Pattern-4K Pieces 0.635 Fragments 0.602 Leads 0.729
These help a bit, particularly with the pieces, but not a whole lot.
Pattern-1K Pieces 0.500 Fragments 0.285 Leads 0.169
Combine-1K Pieces 0.640 Fragments 0.338 Leads 0.186
Pattern-2K Pieces 0.572 Fragments 0.590 Leads 0.715
Combine-2K Pieces 0.689 Fragments 0.608 Leads 0.723
Pattern-4K Pieces 0.635 Fragments 0.602 Leads 0.729
These help a bit, particularly with the pieces, but not a whole lot.
In []:
No comments:
Post a Comment