A couple of substructure search topics

I had a couple of related conversations at Sheffield about substructure searching which seemed like they could be combined into a short and hopefully interesting blog post.
The first conversation actually started before the meeting: John Mayfield (at NextMove Software) submitted a pull request with some improvements to the RDKit's implementation of the vf2 algorithm, which is what is used to actually do substructure matching. John did a presentation about speeding up substructure matches at the meeting and implemented some of the improvements he described in his talk for the RDKit. If you're interested in details, here's the pull request that has been merged and will be part of the next release. As I'll show below, this leads to a significant improvement in the performance of the RDKit's substructure matcher. Thanks John!
Of course, the best way to speed up substructure matching is to not have to do it at all. This is what the RDKit's pattern fingerprint is for, and the second conversation was about how exactly the pattern fingerprint works. Fortunately I'd recently written this up and added it to the RDKit docs. The new documentation is now part of the "RDKit Book" and will be in the online documentation for the next release. Here's the version of that documentation as of the writing of this blog post:

These fingerprints were designed to be used in substructure screening. These are, as far as I know, unique to the RDKit. The algorithm identifies features in the molecule by doing substructure searches using a small number (12 in the 2019.03 release of the RDKit) of very generic SMARTS patterns - like [*]~[*]~[*](~[*])~[*] or [R]~1[R]~[R]~[R]~1, and then hashing each occurence of a pattern based on the atom and bond types involved. The fact that particular pattern matched the molecule at all is also stored by hashing the pattern ID and size. If a particular feature contains either a query atom or a query bond (e.g. something generated from SMARTS), the only information that is hashed is the fact that the generic pattern matched.
For the 2019.03 release, the atom types use just the atomic number of the atom and the bond types use the bond type, or AROMATIC for aromatic bonds).
NOTE: Because it plays an important role in substructure screenout, the internals of this fingerprint (the generic patterns used and/or the details of the hashing algorithm) may change from one release to the next.

Doing this work made me realize that it's been quite a while since I did any benchmarking to see how effective the pattern fingerprint. The last time I looked at this was a series of blog posts in 2013: here, here, and here. There have been a number of changes to the toolkit since then, so it seemed worthwhile to revisit that benchmarking exercise.
Since benchmarking screenout performance and substructure search speeds seems like a useful thing to be able to do, I created a new script that will be part of the RDKit source distribution in future releases; that's here.

Rather than duplicating a bunch of code from the new benchmarking script here, I'll just show how to run it and then talk about the results.
Note that if you want to follow along with this, you will need to download the datasets that are being used (they are big enough that I didn't want to make them part of the RDKit source distro. URLs for where to find the data are in the script's source. The datasets themselves are described in the first blog post mentioned above (though note that the 25K pairs of molecules are pulled from ChEMBL21, not ChEMBL16 as in the original post).
Running the script in Jupyter is pretty easy:

In [2]:

from rdkit import RDConfig

In [8]:

loc = f"{RDConfig.RDBaseDir}/Regress/Scripts"
!cd $loc;python fingerprint_screenout.py

[07:21:19] INFO: mols from smiles
[07:21:27] INFO: Results1:  7.77 seconds, 50000 mols
[07:21:27] INFO: queries from smiles
[07:21:27] INFO: Results2:  0.16 seconds
[07:21:27] INFO: generating pattern fingerprints for mols
[07:21:43] INFO: Results3:  16.11 seconds
[07:21:43] INFO: generating pattern fingerprints for queries
[07:21:43] INFO: Results4:  0.34 seconds
[07:21:43] INFO: testing frags queries
[07:22:03] INFO: Results5:  19.90 seconds. 6753 tested (0.0003 of total), 3989 found,  0.59 accuracy. 0 errors.
[07:22:03] INFO: testing leads queries
[07:22:23] INFO: Results6:  19.77 seconds. 1586 tested (0.0001 of total), 1067 found,  0.67 accuracy. 0 errors.
[07:22:23] INFO: testing pieces queries
[07:23:19] INFO: Results7:  55.37 seconds. 3333202 tested (0.0810 of total), 1925628 found,  0.58 accuracy. 0 errors.
| 2019.09.1dev1 | 7.8 | 0.2 | 16.1 | 0.3 | 19.9 | 19.8 | 55.4 |

Let's start by looking at how effective the screenout is. This is captured in the lines for Results5 (fragments) Results6 (leads) and Results7 (pieces). Comparing the screenout accuracy (this is the fraction of compounds passing the fingerprint screen that actually had a match) we see comparable values to 2013:

	2013	2019
Fragments	0.59	0.59
Leads	0.72	0.67
Pieces	0.57	0.58

This is what I'd hope to see: we haven't made any substantial changes to the fingerprinter itself since 2013 except to improve the performance on single-atom queries of "strange" elements, so the screenout accuracy shouldn't have changed much. Note that we wouldn't expect the results to be exactly the same since the set of molecules being searched through is different.
It's worth pointing out the effectiveness of the pattern fingerprints in reducing the number of substructure searches that actually have to be done here: for the fragments queries only 6753 of the 25 million possible comparisons, 0.03%, actually need to be done. For the leads it's 0.01%, and for the pieces it's 8.1%. Nice!

Now let's look at the impact of John's vf2 changes to the overall runtime. For that we can just look at summary timing information in the last line of the output above and compare it to what I got when I ran the same script using the most recent RDKit release:

| 2019.03.2     | 7.8 | 0.2 | 24.9 | 0.5 | 20.6 | 20.3 | 65.5 |
| 2019.09.1dev1 | 7.8 | 0.2 | 16.1 | 0.3 | 19.9 | 19.8 | 55.4 |

Given the low number of substructure searches run for the fragments and leads queries, the times in columns 6 and 7 are dominated by the fingerprint comparisons, so there's not much difference. The pieces queries, on the other hand, do show a nice improvement: the overall runtime (including the fingerprint screening) drops from 65.5 seconds to 55.4 seconds. The other nice improvement is in the amount of time required to generate the pattern fingerprints for the 50K molecules (column 4): this drops from 24.9 seconds to 16.1 seconds: it's now running in 65% of the time.
In order to get a better feeling for the speedup from the vf2 changes I ran the fingerprint_screenout.py script with the --validate argument; this performs all substructure matches in order to validate that the pattern fingerprints aren't filtering out any true matches (they aren't). That takes a lot longer to run, so I will just show the results I got when I ran it:

| 2019.03.2     | 8.1 | 0.2 | 24.9 | 0.5 | 356.3 | 372.1 | 432.1 |
| 2019.09.1dev1 | 7.9 | 0.2 | 16.2 | 0.3 | 205.1 | 207.4 | 276.4 |

From these results we can see that the substructure searches now run in 56% to 64% of the time. Very nice!

In [ ]: