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.
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.html
As 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.
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.
These sets were discussed in this thread on the mailing list: http://www.mail-archive.com/rdkit-discuss@lists.sourceforge.net/msg02066.html and this presentation: http://www.hinxton.wellcome.ac.uk/advancedcourses/MIOSS%20Greg%20Landrum.pdf
Look at the size distributions:
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.
importing Chem from rdkit on engine(s)
importing DataStructs from rdkit on engine(s)
importing pyAvalonTools from rdkit.Avalon on engine(s)
Verify that it works:
The above runtimes were collected with a cluster of 4 workers, so it looks like the parallization is working.
The test harness
Tests
----------------------------
Doing Avalon-2K
Leads
Found 886 matches in 2044 searches with 0 failures. Accuracy: 0.433
Frags
Found 4479 matches in 13656 searches with 0 failures. Accuracy: 0.328
Summarize that:
Avalon-2K Pieces 0.465 Fragments 0.328 Leads 0.433
Layered Pieces 0.417 Fragments 0.105 Leads 0.187
Pattern-1K Pieces 0.500 Fragments 0.285 Leads 0.169
Pattern-2K Pieces 0.572 Fragments 0.590 Leads 0.715
Pattern-4K Pieces 0.635 Fragments 0.602 Leads 0.729
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:
And now we can get some examples:
Failure: Cc1ccc(C)c(S(=O)(=O)c2nnn3c4ccsc4c(Nc4ccc(C)c(C)c4)nc23)c1 c1cn[nH]n1
Failure: COc1ccc(OCC#Cc2cn([C@H](C)C[C@@H]3CC[C@H]([C@@H](C)C(=O)N(C)Cc4ccccc4)O3)nn2)cc1 c1cn[nH]n1
Failure: Cc1c(C(=O)Nc2ccccc2)nnn1Cc1ccccc1O c1cn[nH]n1
Look at a specific example:
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:
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:
This ends up mattering because the RDKit's substructure matcher ignores the H completely:
Unless, of course, you construct the query from SMARTS:
Other fingerprints
The RDKit fingerprint
Failure: COc1ccc(-n2c(SC(C)=O)nc3sc4c(c3c2=O)CCCC4)cc1 CSC(C)=O
Failure: CCCCCC[C@H]1SC(=O)c2ccccc2[C@@H]1C(=O)O CSC(C)=O
Failure: COc1ccccc1SC(=O)c1cccc(C=O)n1 CSC(C)=O
Out[27]:
(25000000, 1399, 1240, 34)
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:
Found 1274 matches in 125651 searches with 0 failures. Accuracy: 0.010
Out[53]:
(25000000, 125651, 1274, 0)
Using a shorter maxPath length also hurts accuracy:
Found 1240 matches in 2331 searches with 0 failures. Accuracy: 0.532
Out[35]:
(25000000, 2331, 1240, 0)
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:
Failure: CC(C)(C)c1cc2c(cc1SC1=C(O)OC(CCc3ccccc3)(c3ccccc3)CC1=O)CCCC2 C1=COCCC1
Failure: CC(C)(C)c1cc2c(cc1SC1=C(O)OC(CCc3ccccc3)(c3ccccc3)CC1=O)CCCC2 O=C1C=CO[C@H](c2ccccc2)C1
Failure: Cc1ccc(C)c(S(=O)(=O)c2nnn3c4ccsc4c(Nc4ccc(C)c(C)c4)nc23)c1 c1cn[nH]n1
Failure: CCOC(=O)Cc1c(C)n(CC2CCCCC2)c2ccc(OC)cc12 Cc1cc2cc(O)ccc2[nH]1
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:
Failure: CC(C)(C)c1cc2c(cc1SC1=C(O)OC(CCc3ccccc3)(c3ccccc3)CC1=O)CCCC2 C1=COCCC1
Failure: Cc1ccc(C)c(S(=O)(=O)c2nnn3c4ccsc4c(Nc4ccc(C)c(C)c4)nc23)c1 c1cn[nH]n1
Failure: CCOC(=O)C1=COC(C)(CCc2ccccc2)CC1=O C1=COCCC1
Failure: CCOC(=O)Cc1c(C)n(CC2CCCCC2)c2ccc(OC)cc12 Cc1cc2cc(O)ccc2[nH]1
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:
Found 15 matches in 23 searches with 0 failures. Accuracy: 0.652
----------------------------
Doing 1024
Pieces
Found 1935315 matches in 3023990 searches with 0 failures. Accuracy: 0.640
Frags
Found 4659 matches in 13765 searches with 0 failures. Accuracy: 0.338
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.
