Thursday, August 20, 2015

Curating the PAINS filters

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 some bugs, 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.

In [1]:
from rdkit import Chem
import time,random
from rdkit.Chem.Draw import IPythonConsole
IPythonConsole.molSize = (450,250)
from rdkit import rdBase
from __future__ import print_function
%load_ext sql
print(rdBase.rdkitVersion)
print(time.asctime())

2015.09.1.dev1
Sun Aug  9 15:04:15 2015

H atoms in SMARTS

In order for any of this to make sense, it's first important to understand what the RDKit does with H atoms in SMARTS.

Here's one of the SMARTS from the PAINS set:

"[#8]=[#16](=[#8])-[#6](-[#6]#[#7])=[#7]-[#7]-[#1]","<regId=cyano_imine_C(12)>"

Here's a rendering of the pattern:

In [2]:
sma = '[#8]=[#16](=[#8])-[#6](-[#6]#[#7])=[#7]-[#7]-[#1]'
pains8 = Chem.MolFromSmarts(sma)
pains8
Out[2]:

And a molecule that should match:

In [3]:
mol8 = Chem.MolFromSmiles(r'COC(=O)c1sc(SC)c(S(=O)(=O)C(C)C)c1N/N=C(\C#N)S(=O)(=O)c1ccccn1') #CHEMBL3211428
mol8
Out[3]:

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():

In [11]:
patt = Chem.MolFromSmarts('[#6]([#1])[#1]',mergeHs=True)
Chem.MolToSmarts(patt)
Out[11]:
'[#6&!H0&!H1]'

The code can handle recursive SMARTS properly:

In [14]:
patt = Chem.MolFromSmarts('[$([#6]-[#7]),$([#6]-[#1]),$([#6])]',mergeHs=True)
Chem.MolToSmarts(patt)
Out[14]:
'[$([#6]-[#7]),$([#6&!H0]),$([#6])]'

But there's nothing it can do about OR queries (this is a solvable problem but the logic is complex, here's the bug):

In [15]:
patt = Chem.MolFromSmarts('[#6]-[#1,#6]',mergeHs=True)
Chem.MolToSmarts(patt)
Out[15]:
'[#6]-[#1,#6]'

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])]

Back to the PAINS

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:

print("   reading patts")
smas = [x[0] for x in csv.reader(open('./wehi_pains.csv'))]
opatts = [Chem.MolFromSmarts(x,mergeHs=False) for x in smas]
patts = [Chem.MolFromSmarts(x,mergeHs=True) for x in smas]

print("   reading mols")
smis = [x[0] for x in csv.reader(open('./test_data/wehi_mols.csv'))]
ms = [Chem.MolFromSmiles(x) for x in smis]
mhs = [Chem.AddHs(x) for x in ms]

print("   filtering")
matches=[]
found=0
for i,(m,mh) in enumerate(zip(ms,mhs)):
    for j,(patt,opatt) in enumerate(zip(patts,opatts)):
        t1 = m.HasSubstructMatch(patt)
        t2 = mh.HasSubstructMatch(opatt)
        if t1:
            found+=1        
        if t1^t2:
            matches.append((i,j,smis[i],smas[j]))
            print(i,j,smis[i],smas[j])
    if not (i+1)%100:
        print("    Done: ",i+1," matches: ",len(matches)," found: ",found)
print("    Done: ",i+1," matches: ",len(matches)," found: ",found)

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.

Producing a test set

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:

import urllib,requests,json
base = 'http://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/fastsubstructure/smarts/%s/property/CanonicalSmiles/json'
pubchem={}
for key in missing_keys2:
    sma = smas[key]
    url = base%urllib.quote(sma)
    print(sma)
    try:
        r = requests.post(url)
        pubchem[key]=json.loads(r.text)
    except KeyboardInterrupt:
        break
    except:
        import traceback
        traceback.print_exc()

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:

"[#8]-[#6](=[#8])-[#6](-[#1])(-[#1])-[#16;X2]-[#6](=[#7]-[#6]#[#7])-[#7](-[#1])-c:1:c:c:c:c:c:1","<regId=cyanamide_A(1)>"
"c:1-3:c(:c:c:c:c:1)-[#16]-[#6](=[#7]-[#7]=[#6]-2-[#6]=[#6]-[#6]=[#6]-[#6]=[#6]-2)-[#7]-3-[#6](-[#1])-[#1]","<regId=colchicine_het(1)>"
"c:1(:c(:c:2:c(:n:c:1-[#7](-[#1])-[#1]):c:c:c(:c:2-[#7](-[#1])-[#1])-[#6]#[#7])-[#6]#[#7])-[#6]#[#7]","<regId=cyano_amino_het_A(1)>"
"[#6](-[#1])-[#6]:2:[#7]:[#7](-c:1:c:c:c:c:c:1):[#16]:3:[!#6&!#1]:[!#1]:[#6]:[#6]:2:3","<regId=het_thio_N_55(5)>"
"[#6]-2(=[#16])-[#7]-1-[#6]:[#6]-[#7]=[#7]-[#6]-1=[#7]-[#7]-2-[#1]","<regId=thio_urea_K(2)>"
"[#7](-[#1])(-[#1])-c:1:c(:c(:c(:c(:c:1-[#7](-[#1])-[#16](=[#8])=[#8])-[#1])-[#7](-[#1])-[#6](-[#1])-[#1])-[F,Cl,Br,I])-[#1]","<regId=anil_NH_no_alk_B(1)>"
"[#7]-4(-c:1:c:c:c:c:c:1)-[#6](=[#7+](-c:2:c:c:c:c:c:2)-[#6](=[#7]-c:3:c:c:c:c:c:3)-[#7]-4)-[#1]","<regId=het_5_inium(1)>"
"c:1:3:c(:c:c:c:c:1)-[#7]-2-[#6](=[#8])-[#6](=[#6](-[F,Cl,Br,I])-[#6]-2=[#8])-[#7](-[#1])-[#6]:[#6]:[#6]:[#6](-[#8]-[#6](-[#1])-[#1]):[#6]:[#6]:3","<regId=anil_OC_alk_B(3)>"
"[#6]-1(=[#6](-!@[#6]=[#7])-[#16]-[#6](-[#7]-1)=[#8])-[$([F,Cl,Br,I]),$([#7+](:[#6]):[#6])]","<regId=thiaz_ene_C(11)>"
"s:1:c(:c(-[#1]):c(:c:1-[#6]-3=[#7]-c:2:c:c:c:c:c:2-[#6](=[#7]-[#7]-3-[#1])-c:4:c:c:n:c:c:4)-[#1])-[#1]","<regId=het_76_A(1)>"
"[#7]=[#6]-1-[#7](-[#1])-[#6](=[#6](-[#7]-[#1])-[#7]=[#7]-1)-[#7]-[#1]","<regId=het_6_imidate_A(4)>"
"[#6]-2(=[#7]-c1c(c(nn1-[#6](-[#6]-2(-[#1])-[#1])=[#8])-[#7](-[#1])-[#1])-[#7](-[#1])-[#1])-[#6]","<regId=het_65_G(1)>"
"c:1:2(:c(:c(:c(:o:1)-[#6])-[#1])-[#1])-[#6](=[#8])-[#7](-[#1])-[#6]:[#6](-[#1]):[#6](-[#1]):[#6](-[#1]):[#6](-[#1]):[#6]:2-[#6](=[#8])-[#8]-[#1]","<regId=anthranil_acid_I(1)>"
"[#6](-[#1])(-c:1:c(:c(:c(:c(:c:1-[#1])-[#1])-[Cl])-[#1])-[#1])(-c:2:c(:c(:c(:c(:c:2-[#1])-[#1])-[Cl])-[#1])-[#1])-[#8]-[#6](-[#1])(-[#1])-[#6](-[#1])(-[#1])-[#6](-[#1])(-[#1])-c3nc(c(n3-[#6](-[#1])(-[#1])-[#1])-[#1])-[#1]","<regId=misc_imidazole(1)>"
"c2(c-1n(-[#6](-[#6]=[#6]-[#7]-1)=[#8])nc2-c3cccn3)-[#6]#[#7]","<regId=het_65_H(1)>"
"[#7](-[#1])(-c:1:c(:c(:c(:c(:c:1-[#1])-[#1])-[#1])-[#1])-[#8]-[#1])-[#6]-2=[#6](-[#8]-[#6](-[#7]=[#7]-2)=[#7])-[#7](-[#1])-[#1]","<regId=het_6_imidate_B(1)>"
"[#8]=[#6]-3-c:1:c(:c:c:c:c:1)-[#6]-2=[#6](-[#8]-[#1])-[#6](=[#8])-[#7]-c:4:c-2:c-3:c:c:c:4","<regId=quinone_C(2)>"
"c:1:c:c-2:c(:c:c:1)-[#7](-[#6](-[#8]-[#6]-2)(-[#6](=[#8])-[#8]-[#1])-[#6](-[#1])-[#1])-[#6](=[#8])-[#6](-[#1])-[#1]","<regId=misc_aminal_acid(1)>"

These, along with having a set of molecules which are known not to match each PAINS, and looking into the aromaticity model changes in some more detail, are things to come back to. In the meantime, the test set, along with code that runs it, is here: https://github.com/rdkit/rdkit/tree/master/Data/Pains/test_data The updated version of the PAINS filters is here: https://github.com/rdkit/rdkit/blob/master/Data/Pains/wehi_pains.csv

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. :-)

Posted by at 9 comments:
Labels: , ,

Saturday, August 15, 2015

Tuning Substructure Queries

This is inspired by a question from a colleague this week.

I'm going to present a method to exert a bit more control over substructure queries without having to do a massive amount of work.

In [1]:
from rdkit import Chem,DataStructs
import time,random
from collections import defaultdict
import psycopg2
from rdkit.Chem import Draw,PandasTools,rdMolDescriptors
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw,rdqueries
from rdkit import rdBase
from __future__ import print_function
import requests
from xml.etree import ElementTree
import pandas as pd
%load_ext sql
print(rdBase.rdkitVersion)
print(time.asctime())

2015.09.1.dev1
Sat Aug 15 17:54:55 2015

Our starting point is a query against ChEMBL for molecules that contain a phenyl ring with a substituent.

In [2]:
data = %sql postgresql://localhost/chembl_20 \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>'c1ccccc1*'::qmol limit 10;

10 rows affected.
In [3]:
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)
Out[3]:

Though some of the structures above (1593124 and 1594761) contain the fragment we are looking for, the others all include additional substituents on the phenyl.

It's straightforward to make the query more specific by hand-editing the SMARTS, here's one example:

In [4]:
sma = '[cH]1[cH][cH][cH][cH]c1*'
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)            

10 rows affected.
Out[4]:

But that would be something of a pain for more complex queries.

Here's another approach, the idea is to modify the query for each non-wildcard atom and specify that its degree should be what is found in the query:

In [5]:
sma = 'c1ccccc1-[*]'
m = Chem.MolFromSmarts(sma)
# start by creating an RWMol:
m = Chem.RWMol(m)
for atom in m.GetAtoms():
    # skip dummies:
    if not atom.GetAtomicNum():
        continue
    atom.ExpandQuery(rdqueries.ExplicitDegreeEqualsQueryAtom(atom.GetDegree()))
sma = Chem.MolToSmarts(m)
print("Result SMARTS:",sma)

Result SMARTS: [c&D2]1:,-[c&D2]:,-[c&D2]:,-[c&D2]:,-[c&D2]:,-[c&D3]:,-1-*

This also produces the results we want:

In [6]:
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)            

10 rows affected.
Out[6]:

This only works for query atoms (i.e. the molecule must be constructed from SMARTS), but we can adapt it to deal with non-query atoms as well:

In [7]:
sma = 'c1ccccc1-[*]'
m = Chem.MolFromSmiles(sma)
# start by creating an RWMol:
m = Chem.RWMol(m)
for atom in m.GetAtoms():
    # skip dummies:
    if not atom.GetAtomicNum():
        continue
    oa = atom
    if not atom.HasQuery():
        needsReplacement=True
        atom = rdqueries.AtomNumEqualsQueryAtom(oa.GetAtomicNum())
        atom.ExpandQuery(rdqueries.IsAromaticQueryAtom(oa.GetIsAromatic()))
        if(oa.GetIsotope()):
            atom.ExpandQuery(rdqueries.IsotopeEqualsQueryAtom(oa.GetIsotope()))
        if(oa.GetFormalCharge()):
            atom.ExpandQuery(rdqueries.FormalChargeEqualsQueryAtom(oa.GetFormalCharge()))
    else:
        needsReplacement=False
    atom.ExpandQuery(rdqueries.ExplicitDegreeEqualsQueryAtom(oa.GetDegree()))
    if needsReplacement:
        m.ReplaceAtom(oa.GetIdx(),atom)
sma = Chem.MolToSmarts(m)
print("Result SMARTS:",sma)

Result SMARTS: [c&D2]1:[c&D2]:[c&D2]:[c&D2]:[c&D2]:[c&D3]:1-[*]

There's a bit too much book-keeping here, it's kind of screaming for a QueryAtom constructor that takes a standard Atom and does the right thing, but that's not there yet.

This, this already produces something useful:

In [8]:
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)

10 rows affected.
Out[8]:

Let's wrap that in a function and add a couple of additional parameters for even more control

In [9]:
def adjustQuery(m,ringsOnly=True,ignoreDummies=True):
    qm =Chem.RWMol(m)
    if ringsOnly:           
        ri = qm.GetRingInfo()
        try:
            ri.NumRings()
        except RuntimeError:
            ri=None
            Chem.FastFindRings(qm)
            ri = qm.GetRingInfo()
    for atom in qm.GetAtoms():
        if ignoreDummies and not atom.GetAtomicNum():
            continue
        if ringsOnly and not ri.NumAtomRings(atom.GetIdx()):
            continue

        oa = atom
        if not atom.HasQuery():
            needsReplacement=True
            atom = rdqueries.AtomNumEqualsQueryAtom(oa.GetAtomicNum())
            atom.ExpandQuery(rdqueries.IsAromaticQueryAtom(oa.GetIsAromatic()))
            if(oa.GetIsotope()):
                atom.ExpandQuery(rdqueries.IsotopeEqualsQueryAtom(oa.GetIsotope()))
            if(oa.GetFormalCharge()):
                atom.ExpandQuery(rdqueries.FormalChargeEqualsQueryAtom(oa.GetFormalCharge()))
        else:
            needsReplacement=False
        atom.ExpandQuery(rdqueries.ExplicitDegreeEqualsQueryAtom(oa.GetDegree()))
        if needsReplacement:
            qm.ReplaceAtom(oa.GetIdx(),atom)            
    return qm

Try it out:

In [10]:
qm = adjustQuery(Chem.MolFromSmarts('c1ccccc1*'))
sma = Chem.MolToSmarts(qm)
print("Result SMARTS:",sma)
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)

Result SMARTS: [c&D2]1:,-[c&D2]:,-[c&D2]:,-[c&D2]:,-[c&D2]:,-[c&D3]:,-1-,:*
10 rows affected.
Out[10]:

Another use case: a phenyl ring that should have an oxygen substituent meta to another substituent:

In [11]:
sma = 'Oc1cc(*)ccc1'
qm = adjustQuery(Chem.MolFromSmarts(sma))
sma = Chem.MolToSmarts(qm)
print("Result SMARTS:",sma)
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)

Result SMARTS: O-,:[c&D3]1:,-[c&D2]:,-[c&D3](-,:*):,-[c&D2]:,-[c&D2]:,-[c&D2]:,-1
10 rows affected.
Out[11]:

Notice that we're getting ethers and esters as well as alcohols. We can fix that by either adding an H to the O query or by telling adjustQuery() to not ignore non-ring atoms:

In [12]:
sma = 'Oc1cc(*)ccc1'
qm = adjustQuery(Chem.MolFromSmarts(sma),ringsOnly=False)
sma = Chem.MolToSmarts(qm)
print("Result SMARTS:",sma)
data = %sql  \
    select molregno,molfile from rdk.mols join compound_structures using (molregno) where m@>:sma ::qmol limit 10;
mols = [Chem.MolFromMolBlock(y) for x,y in data]
Draw.MolsToGridImage(mols,legends=[str(x) for x,y in data],molsPerRow=4)

Result SMARTS: [O&D1]-,:[c&D3]1:,-[c&D2]:,-[c&D3](-,:*):,-[c&D2]:,-[c&D2]:,-[c&D2]:,-1
10 rows affected.
Out[12]:

The adjustQuery() function is, in some form, going to make it into the RDKit C++ code for the next release so that it's available from the cartridge as well as Python.

I still need to figure out what the API should be (which options should be avaiable). Suggestions are welcome, either as comments here or on the github ticket: https://github.com/rdkit/rdkit/issues/567

In [ ]:
Posted by at No comments:
Labels:
Subscribe to: Posts (Atom)