Thursday, 18 January 2018

Implementing the Sayle tautomer hash with Open Babel

One of the consequences of last year's overhaul of handling of kekulization, aromaticity and implicit hydrogens is the ability to easily calculate structure hashes such as Roger Sayle's tautomer hash, which I wrote up on the NextMove blog a while ago.

I routinely use this hash to handle tautomers, particularly when dealing with R groups. It doesn't solve all tautomer issues (e.g. ones that involve carbon) but it can quickly bring you from having no support for handling tautomers to getting you 95% of the way. In fact, I've been thinking about adding this (and some of the other hashes that Roger has come up with) as cansmi output options. Anyhoo, here's an implementation in Python:
import pybel
ob = pybel.ob

def tautomerhash(smi):
    """Take a SMILES and return the Sayle tautomer hash:
    mol = pybel.readstring("smi", smi).OBMol
    mol.DeleteHydrogens() # just in case
    formalcharges = 0
    hcount = 0
    for atom in ob.OBMolAtomIter(mol):
        formalcharges += atom.GetFormalCharge()
        if atom.GetAtomicNum() != 6: # non-carbon
            hcount += atom.GetImplicitHCount()
    for bond in ob.OBMolBondIter(mol):
    mol.SetAromaticPerceived() # no point triggering perception
    return "%s_%d" % (pybel.Molecule(mol).write("can").rstrip(), hcount-formalcharges)

if __name__ == "__main__":
    smis = ["*c1c(c(C(=N)O)cc2nc([nH]c12)C(=O)[O-])N(=O)=O",
    for smi in smis:
The two SMILES in the example code above are those from the original blog post. Here they give the following two identical hashes (note to self: 'fix' the bracketed asterisk):

Sunday, 17 December 2017

Faster toolkit, faster! Part II

A while ago I described some work I was doing to improve the overall speed of the Open Babel toolkit. Here I want to focus on a particular use-case and where I've got to.

This use case is SMILES to SMILES conversion. Now, while this particular transformation might not sound very interesting, it does encompass both SMILES reading and SMILES writing in one handy package, and both of these operations are often important when dealing with databases or datasets of chemical structures. It also exercises several areas of the toolkit such as kekulization, handling of aromaticity, and stereo perception (or not, as we'll see). Canonicalization is also relevant here, but I didn't do any work on that (and it needs some).

To begin with, some timings. To convert 100K ChEMBL molecules from smi to smi took 10m7s with OB 2.4.1. With the current development version it takes 31s. One change to the defaults is that the dev version does not reperceive the stereo. If you turn on stereo perception (-aS), it takes 1m13s. You can speed things up if you also avoid reperceiving the aromaticity (-aa) and read it as provided in the input; then the conversion only takes 19s.

[21 days later] So I was originally going to describe the results from the Visual Studio profiler for this conversion. But then I said, hey, I might as well fix that one, and that one there, and, well, you know how it goes - this part is actually quite fun, when you make a small change and see the speed go up. Anyway, the conversion that used to take 19s, now takes 11.0s. If you're interested, the speedups included things like replacing std::endl by "\n", caching option values, avoiding string copies, avoiding use of stringstream, avoiding SSSR calculation, and using reserve() on vectors. It was often surprising what things appeared high on the list in the profiler. I can see a few more things that could be improved, but I'm going to leave it there for the moment.

So, in summary, this particular conversion has gone from slow to fast, with a speedup of 55x. There's always more that could be done, but it's respectable.

Open Babel in a snap II

This is a quick update on a previous post, where I mentioned that I had made Open Babel available as a snap. That snap is the released version, 2.4.1. I have now put in place a semi-automated procedure that builds a snap of the development version. You can only have one snap version installed at a time, but you can switch between them.

To install the stable version use:
sudo snap install openbabel

To install the development version use:
sudo snap install openbabel --channel=edge

You can switch between them with:
sudo snap refresh openbabel --channel=stable # or edge

To see which you have installed, use "snap list", or run "openbabel.obabel" and look at the version number and date.

Notes: I'm using Launchpad to do this. Rather than base it directly off the openbabel master (which would require me to check-in snapcraft specific files), I've set it up so that it runs off a branch (named "snaps") in my own repo. Every so often, I merge master into this and a new snap will be created. To fully automate it, I will need to have a cronjob to do that merge automatically.

Thursday, 2 November 2017

Announcing more closures

It is a good idea to avoid extending an existing format. If additional information needs to be stored, many formats provide a means to store it separately so that software that does not support the extension will still be able to read and process the core information. For example, with SMILES, ChemAxon have famously used the title field to store additional info.

However, this is not always possible if one runs up against a fundamental limitation of a format.

I recently ran into a situation where I needed to deal with the fact that the SMILES syntax only supports at most two digit ring closure numbers. This is the %nn notation. The problem is that %111 does not mean ring closure 111 but rather ring closures 11 and 1. If you need a three digit number you are out of luck - it simply can't be expressed as a SMILES string.

So, is it really necessary to have this ability in the first place? Daylight didn't think so. And they were mainly right - when writing SMILES, if you reuse the numbers and pick a traversal order that closes rings as soon as possible then it may be possible to avoid requiring three digits (though a diamond or graphene substructure may cause problems - a refutable hypothesis).

But the thing is, any particular implementation may not be reusing ring closure numbers, or may be favoring a different traversal order (e.g. nicer looking SMILES are made by following ring substituents first). Indeed, remembering that this syntax can be used to create a bond between any two atoms, you may wish to use the ring closure syntax as a means to construct a molecule, and there you don't want to be limited to just 100 bonds; e.g. by stitching together a large molecule by combining shorter SMILES strings (see SmiLib for example, and the note under Restrictions), or via this SMILES hack of Andrew's.

When Bob Hanson of Jmol had to deal with this issue he settled on the syntax "%(number)" (discussed here in 2012, and more recently in J. Cheminf. 2016, 8, 50). One alternative might be "%%", but again this would be limited, this time to three digits. Might as well come up with a general solution as Bob has done. It uses an extra character for the three-digit case but this is going to be rare in the first place, and these SMILES are going to be so big that having a few more characters is not going to break the bank.

Support for this syntax has recently been added to RDKit and Open Babel. For more info see the associated pull requests (here and here). Hopefully we will start to see support for this more generally in other toolkits.

1. According to the Jmol docs, "%(n)" supports any non-negative number (integer, presumably?). My implementation supports up to 5 digits.

Image credit:
image by InExtremiss on Flickr licensed CC-BY

Sunday, 1 October 2017

Open Babel in a snap

Maintainers of Linux distributions do an amazing job packaging applications for users, so that "apt/dnf install openbabel"  makes software packages available for use within seconds. The only downside is that the version of the software may be older than desired, due to the infrequent release schedules of many distros combined with the infrequent release schedule of Open Babel itself.

To bridge the gap for such users, I have created a snap package for the latest release of Open Babel. For those distros that support snapd, "snap install openbabel" will install it and typing "openbabel." and hitting tab a few times will list the executables, e.g. "openbabel.obgui". Um, this is all in theory. If you actually get this to work on a distro other than Ubuntu, please leave a comment below.

1. Snaps run in some sort of security sandbox. As far as I understand, the executables I have created can only read/write files in the user's HOME directory (and subdirectories).
2. Right now, I don't include the libraries or include files. This may come with time.
3. Similarly, right now I've only packaged the release version. Snaps are a good way to distribute nightly-builds or snapshot releases and this might be something to look into in future.

Saturday, 23 September 2017

How many cheminformaticians does it take to read a SMILES string?

A few days ago I posted the following poll question on Twitter:
How many Hs are on the N in the molecule described by this SMILES string? N(C)(C)(C)C
I provided four possible answers:
  • Zero
  • One
  • Depends
  • Can't say as no such molecule
Please take a moment to choose your own answer. To guide you, here is Inigo Montoya.
The results of the survey were interesting. I selected the answers to highlight common misconceptions, and indeed each was chosen at least as much as the correct answer. In fact, only two people got the correct answer, and one of those people shares the office with me.

First of all, it's not "can't say as no such molecule" (8/25 people). To even say whether any such molecule exists, one needs to read the SMILES string, and so could answer the question. Also, the scope of SMILES strings is any molecule that can be represented by a valence model; that is, it is not restricted to only molecules that exist. I think what people were really indicating here is that they work back from the molecule to answer how many hydrogens it must have; given that they couldn't work out what molecule it is, they weren't able to count the hydrogens.

And what about "depends"? (2 people) I'm guessing that people chose this answer to indicate that it's either 0 or 1 (or something else) but it's not possible to be certain. This might be where people are trying to second guess what the corret molecule might have been, given that the one presented is clearly wrong in some way. Here, the key point that people are missing is that a SMILES string exactly represents a molecule including all of its hydrogens.

But what is that hydrogen count? 0 or 1? Well, 13 people went for zero and 2 for one.

Which highlights the final misconception, that the number of hydrogens on an atom in a SMILES string is the same as the answer to "how many hydrogens would I guesstimate if I didn't know how many were attached?". The answer to that is probably zero here. However, when reading SMILES, the number of hydrogens on the nitrogen is exactly described by a set of rules outlined by Daylight at [1], and somewhat more clearly by OpenSMILES at [2] (a project which was spearheaded by former Daylight developer, Craig James). The rules are simple and have no special cases. For nitrogen, the rule is add hydrogens to round up the valence to 3 (if valence <= 3), or up to 5 (if 3 < valence <= 5), and zero otherwise. In short, one hydrogen should be added.

If at this point you're thinking, "but that doesn't make sense - you shouldn't add a hydrogen", you are missing the point. The structure that the SMILES writer was presented with had a hydrogen there. We can argue about whether it was intended, but it did have it. The SMILES writer then used the SMILES valence rule described above to decide whether or not it could omit square brackets. It worked out that it could, and then decided to omit them. Any SMILES reader that respects the rules would then read it in and perfectly reproduce the original structure.

...Or would they? I've just been testing the ability of a number of cheminformatics toolkits and drawing programs to agree on the number of hydrogens on each atom of a given SMILES string. The goal is not to give anyone a kicking, but to improve SMILES interoperability across the board, and in the last two weeks I've given feedback to (I think) 7 different toolkits with this goal in mind. So far I have test results for the Avalon toolkit, BIOVIA Draw, CACTVS, the CDK, ChemDraw, Indigo, IWToolkit, JChem, OEChem, Open Babel, OpenChemLib, the RDKit, and Weininger's CEX. If you would like to add another program to this list, then I encourage you to get in touch.

[1] (section 3.2.1)
[2] (section 3.1.5)

Monday, 28 August 2017

My ACS talk on Kekulization and aromatic SMILES

Here are the slides for the talk I presented last week at the ACS meeting in Washington. It describes my understanding of the Daylight toolkit as deduced by John:
The fundamental reason for attempting to describe the behaviour of the Daylight toolkit is the assumption that this is the correct way to read/write SMILES, and that any deviation is either wrong or at least should be justified. The background to this is that I recently worked on the handling of kekulization and reading SMILES in Open Babel, but found that many of the details were not present in the Daylight documentation. As I say on the final slide, please let me know if you feel I have made any mistake. You can do this by email or by leaving a comment below.

Over at the NextMove blog, Peter Shenkin brought up the biphenylene case, which (to my mind) illustrates alternative approaches to reading aromatic SMILES. Consider the SMILES string c12ccccc1c3ccccc23. Some toolkits may read this, work out that only the two six-membered rings can be aromatic, and then make sure that the double bonds are not placed in the four-membered ring. I refer to this approach as dearomatisation, an approach that Open Babel used to use. It involves ring detection, 4n+2 counting and so forth. Apart from taking some time, an obvious problem is different aromaticity models may be used by the reader and writer, thus leading the reader to drop aromaticity from a particular ring, typically by setting those bonds to single bonds and adjusting hydrogens, resulting in a different structure than intended.

In any case, this is not the approach used by the Daylight toolkit, which did not consider 4n+2, or even detect cycles. The approach is described in the talk above so I won't repeat it here. For the SMILES above, I believe that it would generate one of two Kekulé forms depending on the atom order; one with the two double bonds in the four-membered ring, and one with two benzenes. It's for this reason that Daylight would never generate that SMILES for biphenylene ("don't generate aromatic SMILES that you can't kekulize"), but always write a single bond symbol for the bonds connecting the phenyl rings (e.g. something like c12-c3c(-c2cccc1)cccc3). When written that way, kekulization always gives the desired form.
This use of a single bond symbol is rather subtle. When writing an aromatic SMILES, the rule is "use a single bond symbol when a ring bond is between aromatic atoms but is not itself aromatic". This corrects for the fact that all default bonds joining aromatic atoms are themselves considered aromatic (except where the bond is not in a ring but the atoms are).

Following up on a comment by Rajarshi, while differences in aromaticity models are a problem for 'dearomatisation' algorithms, they are not a problem for the kekulization algorithm used by Daylight. So long as the structure is kekulizable (and appropriate single-bond symbols are used) then it can read in any structure without loss of information no matter what aromatic model is used.