Using vagrant with the RDKit

This is more of a technical post that may be of interest to people doing RDKit development work. It's about getting a vagrant box setup to build and test the RDKit.

Background

There are a lot of people using the RDKit via our KNIME integration. The KNIME nodes are an easy way to access a subset of the RDKit functionality in a nice graphical workflow environment. The nodes themselves, like KNIME, are written in Java and make use the RDKit Java wrappers, so whenever we do a new release of the RDKit, I need to build new versions of those wrappers so that we can update the KNIME nodes. Since we have KNIME users who are stuck using older versions of Linux, I try to maintain compatibility back as far as RHEL5 (released 8.5 years ago, but still being somewhat supported by RedHat). This means that I need to build the Linux version of the Java wrappers with an ancient version of gcc: v4.1. Ubuntu (the Linux distrib that I normally use for RDKit development) hasn't supported gcc v4.1 since the 10.04 release (not still being supported), so I also need an old Linux install. I used to keep a couple of vmware images, one 32bit and one 64bit, around for the builds, but this became a pain, so I wanted to try something else.

Enter Vagrant

Vagrant is a tool for creating reproducible development environments. There's a lot more information, including some excellent documentation, on the other side of that link.

I started by installing Virtual Box on my current development machine (running Ubuntu 15.04). I tracked down a couple of vagrant boxes for Ubuntu 10.04 at https://atlas.hashicorp.com/boxes/search (mrgcastle/ubuntu-lucid32 and f500/ubuntu-lucid64). For each of these I created a local box. Here's the process for the 64bit box:

mkdir lucid64
cd lucid64
vagrant box add f500/ubuntu-lucid64
vagrant init f500/ubuntu-lucid64

This downloads the box image, and configures everything. After editing the resulting Vagrantfile to mount my local RDKit pull (this isn't necessary, I could have just as easily pulled the code from github):

# Share an additional folder to the guest VM. The first argument is
# the path on the host to the actual folder. The second argument is
# the path on the guest to mount the folder. And the optional third
# argument is a set of non-required options.
config.vm.synced_folder "/scratch/RDKit_git/", "/home/vagrant/RDKit_git"

I created a file to provision the machine so that it had everything required to build the RDKit Java wrappers and then brought the machine up:

vagrant up

A couple minutes later I could connect to the machine (with vagrant ssh) and build the Java wrappers.

Between that Vagrantfile and bootstrap.sh, I have everything I need to easily and reproducibly create an VM that is configured and ready to go to build the RDKit. Nice!

Given that such an ancient configuration probably isn't that useful to others, here's a Vagrantfile and bootstrap.sh for a 64bit Ubuntu 12.04 box. This can be used on another Linux machine, a Mac, or Windows (I haven't actually tried Windows).

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

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 [ ]: