Thursday, August 20, 2015

Back in 2010 Rajarshi Guha blogged about converting the PAINS substructure filters from SLN to SMARTS. Given the paucity of tools that effectively support SLN and the usefulness of the PAINS filters, these SMARTS patterns really caught on. Rajarshi ended up also publishing a paper together with J. Baell, one of the original authors of the PAINS paper, that provided KNIME workflows for using the RDKit and Indigo to filter PAINS compounds. A big focus of the paper was differences in the number of matches provided by the two toolkits. The RDKit matched a lot less structures than Indigo did, but neither matched as many as expected.

With the advent of the FilterCatalog functionality in the RDKit, we included the "standard" SMARTS version of PAINS. Axel Pahl quickly pointed out some problems that he had experienced using the PAINS filters and the RDKit. Axel also pointed out that we really didn't have a good test set for the filters.

This all led me to spend some time looking at the SMARTS versions of the PAINS filters and the way the RDKit handles them. It's fortunate that I had a good-sized block of time, because this turned into a larger task than I had anticipated. I ended up making a number of improvements to the ways Hs from SMARTS are handled, fixing somebugs, and modifying a lot of the PAINS SMARTS definitions. This required a bunch of iterative tweaking and testing that I won't describe in detail here, but I think it's worthwhile getting into a bit of what I did and why the changes were necessary.

An aside here: SLN is an interesting format that never really caught on. This set provides a good opportunity for testing and refining the RDKit's SLN support, which has been there for a while at least for standard molecules, but is not extensively tested and would need to be extended to support more query features.

Though we can clearly see that there should be a substructure match, we don't get one:

In [5]:

mol8.HasSubstructMatch(pains8)

Out[5]:

False

The problem is that explicit H in the SMARTS. In order to get a match we need to either add hydrogens:

In [6]:

mol8h=Chem.AddHs(mol8)mol8h.HasSubstructMatch(pains8)

Out[6]:

True

or merge the H atom into the atom it's attached to:

In [7]:

pains8h=Chem.MergeQueryHs(pains8)mol8.HasSubstructMatch(pains8h)

Out[7]:

True

It's important to note that this still matches the molecule with Hs:

In [8]:

mol8h.HasSubstructMatch(pains8h)

Out[8]:

True

So what did this do to the query?

In [9]:

Chem.MolToSmarts(pains8h)

Out[9]:

'[#8]=[#16](=[#8])-[#6](-[#6]#[#7])=[#7]-[#7&!H0]'

MergeQueryHs() finds explicit H atom queries and merges them with the query on the attached atom to add an hydrogen count query. In this case it's specified that the nitrogen has at least one H atom.

The handling of explicit Hs in queries was, I suspect a large part of the reason that the RDKit didn't generate many matches in the KNIME PAINS paper: 391 of the 480 PAINS SMARTS patterns contain an explicit H atom or an explicit H as part of an atom query. Some of those would still generate matches, but many would not.

I'll show a few more examples of how the merging works. Rather than using MergeQueryHs() explicitly in these examples, I will tell the RDKit to perform the merge as part of building the molecule from the SMARTS using the mergeHs argument to MolFromSmartS():

While fixing the PAINS SMARTS definitions (see below) I worked around this shortcoming in the current version of the code by manually editing the affected SMARTS, so something like the above example would become:
[#6;!H0,$([#6]-#6])]

The SLN->SMARTS translation that Rajarshi did resulted in SMARTS that contain explicit Hs. So in order to have the PAINS be useful, I needed to make sure that the H merging was working as well as possible and that

The first step in testing and cleaning up the SMARTS was to find molecules that match when explicit Hs are present, but which don't match when the Hs are implicit. Here's some code I used for doing that:

The idea is simple: find cases where one form of the pattern matches and the other doesn't.

Starting from the SMARTS definitions and 10K test molecules provided as part of the KNIME PAINS paper I did a number of passes of pattern-tweaking and bug-fixing until the above code produced no matches.

The next step, more time consuming, was to repeat the process for a larger set of molecules: the 1.4 million molecules in the ChEMBL20 SDF.

At this point I had a set of SMARTS definitions that produced the same results as on molecules without Hs as the original SMARTS did for molecules with Hs. This had been tested on ChEMBL20 and the 10K test compounds.

Another nice thing to have would be a test set for the PAINS filters: a set of molecules where we know which PAINS filters they should match (and which they should not match!).

Building a good test set is hard and takes more time than I had available, but I wanted to get a start on it by finding at least one molecule that matched each of the PAINS filters.

I started using the 10K test molecules and ChEMBL20. The 10K set produced matches for 144 PAINS, ChEMBL produced matches for 293. Taking duplicates into account I now had examples for 309 of 480 PAINS. Running the set I now had across around 16 million compounds from the ZINC full set turned up another 106 matches, bringing the total to 399; 81 to go!

I figured my best bet to get the remaining 81 was PubChem, but I really didn't feel like pulling down the full set of 40+ million compounds and processing it. So I opted to use the PubChem web services API. Here's a bit of the code I used for that:

This returned matches that led to another round of SMARTS editing, most of which were due to differences in aromaticity models. After finishing the process, I had examples for an additional 63 SMARTS where the RDKit could generate matches. There are 18 left that I still don't have a working example for:

A request: getting this put together was a fair amount of work and it would be very nice to get some feedback on it/credit for it. If you use these, please let me know either via email, posts to the mailing list, a comment here, or a message on github. The data files and tests are under the same BSD license as the rest of the RDKit, so there's no requirement that you do so, but it would be a nice gesture. :-)