Friday, 12 September 2014

Not to scale

The latest release of the CDK (1.5.8) includes a new generator for rendering structure diagrams. A detailed introduction to configuring the new generator is available on the CDK wiki[1].

The new generator can be used as a drop in replacement in existing code. However, one aspect of rendering that I've struggled with previously was getting good sized depictions with the CDK - most notably with vector graphic output. This post will look at how we can size depictions and will provide code in an example project as a reference.

ChEBI's current entity of the month - maytansine [CHEBI:6701] will be used to demonstrate the sizing.

Parameters

Three parameters that are important in the overall sizing of depictions. These are the BondLength, Scale, and Zoom which are all registered as BasicSceneGenerator parameters. The Zoom is not needed if we allow our diagram to be fitted automatically.

The BondLength can be set by the user and has a default value of '40' whilst the Scale is set during diagram generation. BondLength units are arbitrary - for now we'll consider this as '40 px'.

Scaling

The Scale parameter is used to render molecules with different coordinate systems consistently[2,3]. The value is determined using the BondLength parameter and the bond length in the molecule. For maytansine [CHEBI:6701] the median bond length is ~0.82. Again, the units are arbitrary - this could be Angstroms (it isn't).

The Scale is therefore the ratio of the measured bond length (0.82) to the desired bond length (40 px). For this example, the scale is 48.48. The coordinates must be scaled by ~4800% such that each bond is drawn as 40 px long.

Bounds

Now we know our scale (~48.48), how big is our depiction going to be? It depends how we measure it. One way would be to check the bounding box that contains all atom coordinates (using GeometryUtil.getMinMax()). However, this does not consider the positions of adjuncts and would lead to parts of the diagram being cut off[4].

Fortunately the new generator provides a Bounds rendering element allowing us to determine the true diagrams bounds of 8.46x8.03. Since the scale is ~48.48 the final size of the depiction would be ~410x390 px. A margin is also added.


Current Rendering API

Now we have the size of our diagram we can render raster images. Unfortunately the current rendering API makes this a little tricky as the diagram is generated after the desired image size is provided by the user. To get the correct size we need to generate the diagram twice (to get the bounds) or use an intermediate representation (we'll see this later).

Code 1 - Current rendering API
// structure with coordinates
IAtomContainer container = ...; 

// create the renderer - we don't use a font manager
List<IGenerator<IAtomContainer>> generators = 
        Arrays.asList(new BasicSceneGenerator(),
                      new StandardGenerator(new Font("Verdana", Plain, 18));
AtomContainerRenderer renderer = new AtomContainerRenderer(generators, null); 

Graphics2D g2 = ...; // Graphics2D to draw raster / vector graphics
Rectangle2D bounds = ...; // need the bounds here!
renderer.paint(new AWTDrawVisitor(g2),
               bounds);

Vector graphics

To render scalable graphics we can use the VectorGraphics2D[5] implementations of the Java Graphics2D class. Vector graphics output can use varied units (e.g. pt, mm, px) - the VectorGraphics2D uses mm.

Without adjusting our scaling the render of maytansine [CHEBI:6701] would be displayed with bond lengths of 40 mm and a total size of ~410x390 mm. The output can be rescaled after rendering but the default width of 41 cm is a bit large. We therefore need to change our desired bond length.

The bond length of published structure diagrams varies between journals. A common and recommended style for wikipedia [6] is 'ACS 1996' - the style has a bond length of '5.08 mm'. Although setting the BondLength parameter to '5.08' would work, other parameters would need adjusting such as Font size (which is provided in pt!).

To render the diagram with the same proportions as the raster image we can instead resize the bounds and fit the diagram to this. Since the desired bond length is '5.08 mm' instead of '40 mm' we need rescale the diagram by 12.7 %. Our final diagram size is then ~52x50 mm. The border for ACS 1996 is '0.56 mm' which can be added to the diagram size.

Example code

To help demonstrate the above rendering I've put together a quick GitHub project johnmay/efficient-bits/scaled-renders. The code provides a convenient API and a command line utility for generating images.

Code 2 - Intermediate object
// structure with coordinates
IAtomContainer container = ...; 

// create the depiction generator
Font font = new Font("Verdana", Plain, 18);
DepictionGenerator generator = new DepictionGenerator(new BasicSceneGenerator(),
                                                      new StandardGenerator(font));

// generate the intermediate 'depiction'
Depiction depiction = generator.generate(container);

// holds on to the rendering primitives as well as the size
double w = depiction.width();
double h = depiction.height(); 

// draw at 'default' size
depiction.draw(g2, w, h);

// generate a PDF (or SVG)
String pdfContent = depiction.toPdf(); // default size
String pdfContent = depiction.toPdf(1.5); // 1.5 * default size
String pdfContent = depiction.toPdf(0.508, 0.056); // bond length, margin

The command line utility provides several options to play with and can load from molfile, SMILES, InChI, or name (using OPSIN[7]).

Code 3 - Command line
# In the project root set the following alias
$: alias render='mvn exec:java -Dexec.mainClass=Main'

# Using OPSIN to load porphyrin and generate a PDF
$: render -Dexec.args="-name porphyrin -pdf"

# Highlight one of the pyrrole in porphyrin
$: render -Dexec.args="-name porphyrin -pdf -sma n1cccc1"

# Show atom numbers
$: render -Dexec.args="-name porphyrin -pdf -atom-numbers"

# Show CIP labels
$: render -Dexec.args="-name '(2R)-butan-2-ol' -pdf -cip-labels"

# Generate a PDF / SVG for ethanol SMILES
$: render -Dexec.args="-smi CCO -pdf ethanol.pdf -svg ethanol.svg"

# Load a molfile
$: render -Dexec.args="-mol ChEBI_6701.mol -pdf chebi-6701.pdf"

You can even play with the font

$: render -Dexec.args="-name 'caffeine' -svg cc-caffeine.svg -font-family 'Cinnamon Cake' -stroke-scale 0.6 -kekule"

Links/References

  1. https://github.com/cdk/cdk/wiki/Standard-Generator
  2. http://gilleain.blogspot.co.uk/2010/09/scaling-and-text.html
  3. http://gilleain.blogspot.co.uk/2010/09/consistent-zoom-with-models-of.html
  4. http://onlinelibrary.wiley.com/doi/10.1002/minf.201200171/abstract
  5. http://trac.erichseifert.de/vectorgraphics2d/
  6. http://en.wikipedia.org/wiki/Wikipedia:Manual_of_Style/Chemistry/Structure_drawing
  7. http://opsin.ch.cam.ac.uk/

Thursday, 11 September 2014

CDK Release 1.5.8

10.5281/zenodo.11681

CDK 1.5.8 has been released and is available from sourceforge (download here) and the EBI maven repo (XML 1).

The release notes (1.5.8-Release-Notes) summarise and detail the changes.

XML 1 - Maven POM configuration
<repository>
  <url>http://www.ebi.ac.uk/intact/maven/nexus/content/repositories/ebi-repo/</url>
</repository>
...
<dependency>
  <groupId>org.openscience.cdk</groupId>
  <artifactId>cdk-bundle</artifactId>
  <version>1.5.8</version>
</dependency>

Tuesday, 22 July 2014

Polish-ed SMARTS parsing

As introduced in previous posts, SMARTS is a concise notation for describing chemical substructure queries. There are several aspects to a SMARTS implementation: subgraph graph matching, parsing, generating, and even optimisation[1,2].

In this post I'll show a way of parsing the binary atom expressions that I found quite neat.

Preliminaries

Conceptually, a SMARTS atom expression is composed of primitives and operators (binary and unary). The primitives test whether some property of a atom (e.g. element, charge, valence, etc) has a certain value[3]. The operators invert and combine these primitives through conjunction (and), disjunction (or), and negation (not).

Some examples of atom expressions are:

[O&X1]
[!C&!N]
[C,c;X3&v4]
[N&!H0&X3]
[!#6&X4]
[O,S,#7,#15]
[C&X3;$([H2]),$([H1][#6]),$(C([#6])[#6])]

The operators in these expressions ordered by their precedence are:

! unary not
& binary and (high)
, binary or
; binary and (low)

The default operator is '&' and can often be omitted such that the first pattern would read [OX1]. There are two 'and' operators with difference precedence allowing logical expressions like:

[C,N&X1]  C or (N and X1)
[C,N;X1]  (C or N) and X1

More complex expression trees can be accomplished with recursive SMARTS.

A formal grammar for SMARTS that respects precedence looks something like this (lifted from the CDK javacc implementation):

SMARTS EBNF grammar
AtomExpression    ::= "[" <LowAndExpression> "]"
LowAndExpression  ::= <OrExpression> [ ";" <LowAndExpression> ]
OrExpression      ::= <HighAndExpression> [ "," <OrExpression> ]
HighAndExpression ::= <NotExpression> [ '&' <HighAndExpression> ]
NotExpression     ::= [ "!" ] <AtomPrimitive>

Notice this is a recursive procedure where I ascend up the precedence hierarchy while descending into the grammar. The small number of operators in SMARTS means this is generally good enough. However there is a non-recursive alternative.

Reverse Polish notation

Reverse Polish notation (RPN) is a notation where the operator follows the operands of an expression[4]. Some simple mathematical expressions are written as follows:

5 + 1              5 1 +
3 + 4 * 2          3 4 2 * +
(3 + 4) * 2        3 4 + 2 *

RPN is extremely useful and simple for implementing and performing operations on stack-based machines[5]. An excellent property is that the operators are applied as soon as they are encountered. Notice that I don't need parentheses to change the multiply and addition order. Also notice that a lookahead check for operator validity isn't needed, when an operator is applied the primitives have already been parsed.

SMARTS operators are infix but let us see what RPN SMARTS might look like:

[O&X1]             O X1 &
[!C&!N]            C ! N ! &
[C,c;X3&v4]        C c , X3 v4 & ;
[N&!H0&X3]         N H0 ! X3 & &
[!#6&X4]           #6 ! X4 &
[O,S,#7,#15]       O S #7 #15 , , ,

RPN SMARTS is much simpler to write a parser for that respect precedence. All that is needed is a way to convert from infix to postfix. The Shunting-yard algorithm[6] does just that.

Implementation

The Shunting-yard algorithm is explained well in many other webpages so I'll neglect that here. I will be converting from infix to postfix and build the expression at the same time. To do this, two stacks are needed, one for atom primitives and one for operators. The atom primitive stack is essentially the output of the Shunting-yard but I apply the operators instead of appending them to the postfix string.

Code 1 consumes characters from the input and either shunts an operator or parses a primitive. Once all the input is consumed the remaining operators are applied. The created query atom is on top of the stack and is returned. A low precedence no-operator is pushed on the stack to make thinks simpler and buffer the shunting.

To handle the implicit '&' between primitives a little more work is needed. Essentially one would optionally invoke shunt(atoms, operators, '&'); as needed at each iteration.

Code 1 - Primary loop
IQueryAtom parse(CharBuffer buffer) throws IOException {

    // stacks of atom primitives and operators
    Deque<IQueryAtom> atoms     = new ArrayDeque<>();
    Deque<Character>  operators = new ArrayDeque<>();
    operators.push(Character.MAX_VALUE); // a pseudo low precedence op

    while (buffer.hasRemaining()) {
        char c = buffer.get();
        if (isOperator(c)) // c == '!' or '&' or ',' or ';'? 
           shunt(atoms, operators, c);
        else                
           atoms.push(parsePrimitive(buffer));
    }

    // apply remaining operators
    while (!operators.isEmpty())
        apply(operators.pop(), atoms);

    return atoms.pop();
}

Code 2 shows the creation of query atom primitives, here they are delegated to several self explanatory utility methods. For compactness only a subset of primitives are read.

Code 2 - Parsing selected primitives
IQueryAtom parsePrimitive(CharBuffer buffer) throws IOException {
    switch (buffer.get(buffer.position() - 1)) {
        case 'A': return newAliphaticQryAtm();
        case 'C': return newAliphaticQryAtm(6);
        case 'N': return newAliphaticQryAtm(7);
        case 'O': return newAliphaticQryAtm(8);
        case 'P': return newAliphaticQryAtm(15);
        case 'S': return newAliphaticQryAtm(16);

        case 'a': return newAromaticQryAtm();
        case 'c': return newAromaticQryAtm(6);
        case 'n': return newAromaticQryAtm(7);
        case 'o': return newAromaticQryAtm(8);
        case 'p': return newAromaticQryAtm(15);
        case 's': return newAromaticQryAtm(16);

        case '#': return newNumberQryAtm(parseNum(buffer));
        case 'X': return newConnectivityQryAtm(parseNum(buffer));
        case 'H': return newHydrogenCountQryAtm(parseNum(buffer));
        case 'R': return newRingMembershipQryAtom(parseNum(buffer));
        case 'v': return newValenceQryAtom(parseNum(buffer));
    }
    throw new IOException("Primitive not handled");
}

To apply an operator, take the operands (primitives) off the top of the atom stack, create a new query atom, and push it back on to the stack (Code 3). If there aren't enough operands, the expression is invalid (not shown).

Code 3 - Applying an operator
void apply(char op, Deque<IQueryAtom> atoms) {
    if (op == '&' || op == ';')
        atoms.push(and(atoms.pop(), atoms.pop()));
    else if (op == ',')
        atoms.push(or(atoms.pop(), atoms.pop()));
    else if (op == '!')
        atoms.push(not(atoms.pop()));
}

Finally, to handle the operator (Code 4), check if the operator currently on top of the stack has precedence over the new operator. If so, pop it from the stack and apply it. The new operator is then added to the stack. Conveniently the code point of the operator character can be used as the precedence.

Code 4 - Handling operator precedence
void shunt(Deque<IQueryAtom> atoms, Deque<Characters> operators, char op) {
    while (precedence(operators.peek()) < precedence(op))
        apply(operators.pop(), atoms);
    operators.push(op);
}

static int precedence(char c) {
    return c; // in ASCII, '!' < '&' < ',' < ';' 
}

With the exception of a few utility methods these four snippets are essentially the whole implementation. You can find the fully functional code on the GitHub project[7].

Not only is the code is incredibly compact and elegant but it can easily be expanded. Several convenience extensions to SMARTS have been made in the past – for example, #X for !#1!#6. A common requirement in general expressions and the Shunting-yard is to handle parenthesis. These need special treatment but it is only a simple modification to the shunting and the precedence value (Code 5).

Code 5 - Handling parenthesis
void shunt(Deque atoms, Deque operators, char op) {
    if (op == ')') {
        while ((op = operators.pop()) != '(')
            apply(op, atoms);
    } else {
        if (op != '(') {
            while (precedence(operators.peek()) < precedence(op))
                apply(operators.pop(), atoms);
        }
        operators.push(op);               
    }
}

int precedence(char c) {
    switch (c) {
        case '!': return 1;
        case '&': return 2;
        case ',': return 3;
        case ';': return 4;
        case '(':
        case ')': return 5;
        default:  return 6;
    }
}

The parser will now correctly handle the following expressions without recursive SMARTS:

[!(C,N,O,P,S)]              C N O P , , , !
[!(C,N,O&X1)]               C N O X1 & , , !
[((C,N)&X3),((O,S)&X2)]     C N , X3 & O S , X2 & ,

All source code is available at github/johnmay/efficient-bits/polished-smarts.

References

  1. PATSY, NextMove Software
  2. SMARTS Optimisation & Compilation: Introduction & Optimisation, Tim Vandermeersch
  3. Daylight theory manual, Daylight CIS
  4. Reverse Polish notation, Wikipedia
  5. Reverse Polish notation and the stack, Computerphile
  6. Shunting-yard algorithm, Wikipedia
  7. github/johnmay/efficient-bits/polished-smarts

Friday, 18 July 2014

CDK Release 1.5.7

CDK 1.5.7 has been released and is available from sourceforge (download here) and the EBI maven repo (XML 1).

The release notes (1.5.7-Release-Notes) summarise and detail the changes. Among the new bug fixes and features, several plugins have been added to the build. The release notes describe how these plugins can be run and what they do so be sure check the notes out if you're a contributor.

XML 1 - Maven POM configuration
<repository>
  <url>http://www.ebi.ac.uk/intact/maven/nexus/content/repositories/ebi-repo/</url>
</repository>
...
<dependency>
  <groupId>org.openscience.cdk</groupId>
  <artifactId>cdk-bundle</artifactId>
  <version>1.5.7</version>
</dependency>

Wednesday, 26 March 2014

Mischievous SMARTS Queries

Last year I extended the CDK SMARTS implementation to match component groupings and stereochemistry. Specifying stereochemistry presents some interesting logical predicate that might be tricky to handle.

Here are some examples that I came up with for testing the correctness of query handling. They start simple before getting a little mischievous. First, recursion and component grouping.

querytargetsnmatchComment
Component grouping (fragment)
(O).(O)O=O0Example from Daylight
OCCO0
O.CCO2
Component grouping (connected)
(O.O)O=O2Example from Daylight
OCCO2
O.CCO0
Recursion, ad infinitum
[$(CC[$(CCO),$(CCN)])]CCCCO1
CCCCN1
CCCCC0
Recursive component grouping
[O;D1;$(([a,A]).([A,a]))][CH]=OOC=O.c1ccccc11Feature/Bug #1312
OC=O0

These next ones are concerned with logic and stereochemistry.
querytargetsnmatchComment
Ensure local stereo matching
*[@](*)(*)(*)O[C@](N)(C)CC12tetrahedrons have 12 rotation symmetries
O[C@@](N)(C)CC12
O[C](N)(C)CC0
Implicit (hydrogen or lone-pair) neighbour
CC[S@](C)=OCC[S@](C)=O1
CC[S@@](C)=O0
CC[S](C)=O0
Either (tetrahedral)
CC[@,@@](C)OCC[C@H](C)O1
CC[C@@H](C)O1
CCC(C)O0
Both (tetrahedral)
CC[@&@@](C)OCC[C@H](C)O0
CC[C@@H](C)O0
CCC(C)O0
Respect logical precedence 1
CC[@,Si@@](C)OCC[C@H](C)O1
CC[C@@H](C)O0
CCC(C)=O0
Respect logical precedence 2
CC[C@,Si@@](C)OCC[C@H](C)O1
CC[C@@H](C)O0
CCC(C)O0
CC[Si@H](C)O0
CC[Si@@H](C)O1
CC[Si](C)O0
Unspecified
CC[@@?](C)OCC[C@H](C)O0
CC[C@@H](C)O1
CCC(C)O1
Negation
CC[!@](C)OCC[C@H](C)O0!@@ is also equivalent to @?
CC[C@@H](C)O1
CCC(C)O1
Neither (tetrahedral) using 'or unspecified'
CC[@?@@?](C)OCC[C@H](C)O0
CC[C@@H](C)O0
CCC(C)O1
Neither (tetrahedral) using negation
CC[!@!@@](C)OCC[C@H](C)O0
CC[C@@H](C)O0
CCC(C)O1
Either (geomeric)
C/C=C/,\CC/C=C/C1
C/C=C\C1
CC=CC0
Neither (geomeric)
C/C=C!/!\CC/C=C/C0
C/C=C\C0
CC=CC1
The last two are quite tricky (and not currently implemented) but once the atom-centric handling is correct it's a simple reduction. It's quite fun to work out so i'll leaf that up to the reader.