Since it was written the UniversalIsomorphismTester (UIT) has been the goto class in the Chemistry Development Kit (CDK) to check structure identity, substructure and overlap matching. The utility works really well but derives the substructure from the more expensive MCS computation. There has been several efforts to bring faster structure matching through Ullmann [1] and the Vento-Foggia (VF) [2] algorithms. Most notably VF was used for structure matching in OrChem [3] and SMSD [4]. The original SMSD code was added to CDK but since then SMSD has been completely rewritten and optimised. I recently tried reintegrating the latest version but have had several failed attempts.
Substructure and identity checking is fundamental and I felt it needed a more centralised home in the code base. This weekend I wrote the first parts of the API and implementations which I will discuss here.
Lazy Matching
For brevity I won't go over how substructure matching works, Rich Apodaca has two excellent introductions on One of These Things is Not Like The Others and Substructure Search From Scratch which gives the key ideas. In the first post it is mentioned that recursion is at the heart of any implementation. Generally recursive implementations are one-or-all, that is you can get the first match or all matches.
As recursion is simply storing frames on the call stack for matching you can easily replace the recursion with a non-recursive stack which only stores each candidate matchings. Providing front end access to this allows a lazy substructure matcher which improves memory and has another benefit we'll see later.
As recursion is simply storing frames on the call stack for matching you can easily replace the recursion with a non-recursive stack which only stores each candidate matchings. Providing front end access to this allows a lazy substructure matcher which improves memory and has another benefit we'll see later.
So far I've added three matchers VentoFoggia (substructure/identity) and Ullmann (substructure) with scope for others [5].
API
The new structure matching paradigm is similar to Java regular expressions - a Pattern is created for a query and then that pattern can be tested on multiple targets. Whether the pattern is matching SMARTS queries, identity or substructure doesn't matter. Personally I've found my self parsing the query and target molecules in the wrong order with the existing matchers and using the pattern makes it impossible to get the order wrong.
Code 1 - Matching API predicate
IAtomContainer query = ...;
Pattern pattern = Ullmann.findSubstructure(query);
for (IAtomContainer target : targets)
if (pattern.matches(target))
hits++;
With language features due in Java 8 the matching hits can even be done in one line.List<IAtomContainer> hits = targets.filter(Ullmann.findSubstructure(query)::matches);The first atom-mapping can be obtained with match() and all atom mappings with matchAll() (which is lazy).
Code 1 - Matching API, lazy loop over all matches
for (IAtomContainer target : targets)
for (int[] mapping : pattern.matchAll(target))
hits++;
Counting Self Matches
To test the performance I first looked at counting self matches. You'd generally want to do this with a partition refinement procedure (see. Using the CDK's group module) but it serves as a good test for the matchers. To compare the new matchers I've used the CDK's UniversalIsomorphismTester (UIT) and the latest SMSD's Substructure (VF implementation). It should be noted that SMSD does go a lot faster with an API change (see. Benchmarking Substructure Search) but this has not yet been added.
 |
| ferrocene as 200 isomorphisms to it's self |
The average (r=10) time taken to match and count isomorphisms (n) was tested on several symetric structures. All structures were normalised to only have sigma bonds. The SMSD implementation did not get the correct number of matches which I think is a bug propagated from the original ChemKit implementation. If any matcher took longer than 10 minutes I considered it as timed out (shown as n/a).
| Structure | Matched | UIT | SMSD | Ullmann | VentoFoggia |
|---|---|---|---|---|---|
| n | t (s) | t (s) | t (s) | t (s) | |
| ferrocene | 200 | 0.281 | 0.166 (n=180?) | 0.007 | 0.001 |
| fullerene C60 | 120 | 279 | 7.5 (n=1) | 1.2 | 0.017 |
| fullerene C70 | 60 | n/a | 9.8 (n=1) | 2 | 0.058 |
| fullerene C320 | 120 | n/a | n/a | n/a | 0.527 |
VentoFoggia only took half a second to find all matches in fullerene C320. I was curious to see how far it could go. The new matchers allow disconnected structures and matching multiple cyclopropanes we quickly reach a huge number of matches.
 |
| C1CC1.C1CC1.C1CC1.C1CC1.C1CC1.C1CC1 has more than 33 million isomorphisms |
| Structure | Matched | Ullmann | VentoFoggia |
|---|---|---|---|
| (SMILES) | n | t (s) | t (s) |
C1CC1 |
6 | 0 | 0 |
C1CC1.C1CC1 |
72 | 0.002 | 0.002 |
C1CC1.C1CC1.C1CC1 |
1,296 | 0.002 | 0.005 |
C1CC1.C1CC1.C1CC1.C1CC1 |
31,104 | 0.168 | 0.037 |
C1CC1.C1CC1.C1CC1.C1CC1.C1CC1 |
933,120 | 6.5 | 1.1 |
C1CC1.C1CC1.C1CC1.C1CC1.C1CC1.C1CC1 |
33,592,320 | 270.5 | 36.3 |
Dataset Searching
As a more general gauge of performance the NCI Aug 2000 (~250,000 structures) was searched for structures containingO=C1C=CCC=C1. It's a single isolated test so conclusions can't really be drawn but it provides a sense for the timings.  |
| Query structure matched against ~250,000 structures in the NCI Aug 2000 data set. |
| Implementation | t first | t all |
|---|---|---|
| UIT | 43.4 s | 43.6 s |
| SMSD | 28.6 s | 4 m 40 s |
| Ullmann | 6.3 s | 6.5 s |
| VentoFoggia | 5.1 s | 4.9 s |
Stereochemistry Matching
The new matchers also check the tetrahedral and geometric (double-bond) stereo configurations are consistent. Technically it's more efficient to check stereo configurations as the vertices are mapped but I decided to go for a post-matching check. This is more modular and also means we can use the same stereo consistency checking for other matchers (the existing UIT). This is where the lazy matching is really useful as we can get the first match which is stereochemically valid without finding all matchings.The stereo checking also seems to be fast enough and I've left it on by default (including all times in this post).
Using the VentoFoggia implementation, it takes ~38 seconds to find all 4,306,056 matches for C(C)(C)(C)(C) in ChEMBL 17. Matching one with a stereo configuration ([C@](C)(C)(C)(C)) is actually just as quick, taking ~37 seconds to find 1,007,652 matches.
Finally I'll leave you with a sneak peak at the Smiles class I have planned which makes loading and matching molecules from SMILES concise and elegant.
Pattern pattern = VentoFoggia.findSubstructure(Smiles.fromString("[C@](C)(C)(C)(C)"));
for (IAtomContainer target : Smiles.fromPath("/db/chembl/chembl_17.smi.gz")) {
hits += Iterables.size(pattern.matchAll(target)); // Guava utils
}
Code
The matcher code is currently awaiting review on the branch feature/isomorphisms. If the branch is missing then it will be in the package org.openscience.cdk.isomorphism.
References
1. Ullmann JR. An Algorithm for Subgraph Isomorphism. Journal of the Association for Computing Machinery. 20042. Cordella LP et al. A (Sub)Graph Isomorphism Algorithm for Matching Large Graphs. IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE. 2004
3. Rijnbeek M and Steinbeck C. OrChem - An open source chemistry search engine for Oracle. Journal of Cheminformatics. 2009
4. Rahman SA et al. Small Molecule Subgraph Detector (SMSD) toolkit. Journal of Cheminformatics. 2009
5. Gouda K and Hassaan, M. A fast algorithm for subgraph search problem. Informatics and Systems (INFOS). 2012
