Improved SMILES Substructure Searching

1. Abstract

2. The Merlin substructure search engine

If the query is a complete valid molecule, its fingerprint can be determined and rapidly used to screen the database. One or more bits that are set in the query fingerprint that are not set in the database fingerprint indicates that the query cannot be a substructure of that database entry. Similarly, a SMILES query can be used in a character frequency screen to also prune entries that cannot match. For example, if the query is bromine, and database entry that whose SMILES does not contain atleast one "B" and one "r" is rejected. A database entry that passes both of these screens then requires an explicit pattern match to confirm the existance of the substructure. This requires a dt_smilin and dt_match of each molecule. Although this final stage is relatively quick, thousands of molecules per second, parsing the molecule, SSSR ring detection, aromaticity perception and pattern matching are much slower than fingerprint screening, which can process millions of molecules per second.

On the other hand, generic pattern matching using SMARTS queries is much slower. Unlike SMILES queries, most SMARTS patterns cannot be fingerprinted or character frequency screened, hence every entry in a database must be processed and pattern matched.

To salvage what it can of the situation, Merlin performs a simple check to determine whether a SMARTS query can be transformed into a syntactically valid SMILES, and if so treats it like a SMILES query. This involves deleting all atom and bond expressions that are not SMILES primitives and interpreting the resulting string as a SMILES.

3. Fragment Fingerprints

Examples of fragment fingerprints are shown below:

Fragment  256-bit fingerprint [ASCII representation]
c         .........U...............0..........U....E..2
c:c       ...0.....U...............0..........U....F2.2
c-c       .........U..........E....8..........U....G..2
c1ccccc1  .0.02....U..........6..3.G....2...U.U....F2.2
C=C       .....U......U.............+..........UE...k.2

To demonstrate the improvement, the time to query 126705 structures in NCI95 with the substructure 'c-c' drops from 1:39.47s to 8.73s using the DayCart Oracle cartridge on a 2 CPU 360MHz SUN Ultra60.

4. SMARTS fingerprint enumeration

To describe the steps of the algorithm, consider the case study pattern 'O=[C,N]aa[N,O;!H0]'. This pattern for intra-molecular hydrogen bonds was included in Dave's Daylight cartridge benchmarks presented in his EuroMug98 talk (here). This was consider a pathological case, where only a single 'O' could be fingerprinted and character frequency screened.

Step #1: SMARTS Optimization The first step is to apply SMARTS optimization techniques as first described in my EuroMug97 talk (here). This converts the SMARTS to 'O=[C,N]-a:a-[N!H0,O!H0]', by recognizing that default bonds to atom expressions that must be aliphatic must be single, and that only a single non-aromatic bond is permitted to an aromatic atom.

Step #2: Enumeration Normal Form The next stage is to strip any predicates unsuitable for enumeration or not applicable to fingerprint screening. The internal representation then contains only atomic number and aromaticity atom predicates, and bond type bond predicates. This results in 'O=[C,N]-[c,n,o,s,p,as,se]:[c,n,o,s,p,as,se]-[N,O]'.

Step #3: Enumeration Refinement To prevent a combinatorial explosion, the next stage iteratively removes high degree atoms or bonds until the user-specified limit is reached. For example, this could reduce the above SMARTS to 'O=[C,N].[N,O]'.

Step #4: Pattern Enumeration The next stage is the actual enumeration and fragment fingerprinting (see above). For this case study, this results in four fingerprints 'O=C.N', 'O=C.O', 'O=N.N' and 'O=N.O'.

Step #5: Fingerprint Reduction The final stage reduces the number of fingerprints by eliminating identical and redundant fingerprints. For our example, the fragment fingerprints for 'O=N.O' and 'O=N.N' are identical; the same as the fingerprint for 'O=N'. Additionally, 'O=C.O' is a subset of the fingerprint for 'O=C.N'. Hence, for our case study, only two fingerprints need be screened 'O=C' and 'O=N'.

In the explanation above the enumeration threshold was set artificially low. In practice, the DayCart Oracle cartridge uses a limit of 1024 fingerprints, and the SMARTS query above generates 196 fingerprints to screen. On a 2 CPU 360MHz SUN Ultra/60, the resulting search time drops from 1m09.29s to only 16.96s.

5. SMILES vs SMARTS semantics

In version 4.62 of the Daylight smiles toolkit the above molecule contains three aromatic rings, and a central non-aromatic ring. The central ring contains two pyrole-like (not pyridine-like) aromatic nitrogens and hence doesn't contain Huckel's 4n+2 electrons. This is a rare example of a ring containing only aromatic atoms that isn't itself aromatic.

Fortunately, in version 4.71, this problem disappears. The new and improved aromaticity perception algorithm now recognizes that the circuit around the perimeter of the molecule contains 18 (4n+2) electrons, and so now every bond in this particular example is aromatic.

6. On-the-fly Code Generation

int fingertest( unsigned int *qry, unsigned int *db, int len )
{
    while( len-- ) {
        if( *qry != (*qry & *db) )
            return False;
        qry++;
        db++;
    }
    return True;
}

int fingertest( unsigned int *db )
{
    if( (db[ 0]&1074348120) != 1074348120 ) return False;
    if( (db[31]& 168839168) !=  168839168 ) return False;
    if( (db[47]&      1072) !=       1072 ) return False;
    if( (db[19]& 536936452) !=  536936452 ) return False;
    if( (db[ 2]&       514) !=        514 ) return False;
    if( (db[45]&     16384) !=      16384 ) return False;
    if( (db[21]&      8192) !=       8192 ) return False;
    return True;
}

The following assembly language shows the binary machine code generated to test the above fingerprint on the Intel 386 architecture.

8b542404        movl 4(%esp),%edx
8b02            movl (%edx),%eax
2558400940      andl $1074348120,%eax
3d58400940      cmpl $1074348120,%eax
7403            je +3
31c0            xorl %eax,%eax
c3              ret
8b427c          movl 124(%edx),%eax
250048100a      andl $168839168,%eax
3d0048100a      cmpl $168839168,%eax
7403            je +3
31c0            xorl %eax,%eax
c3              ret
...
f642b540        testb $64,181(%edx)
7503            jne +3
31c0            xorl %eax,%eax
c3              ret
f6425520        testb $32,85(%edx)
7503            jne +3
31c0            xorl %eax,%eax
c3              ret
b801000000      movl $0x1,%eax
c3              ret

On an 195MHz R10000 SGI Origin, the time taken to screen the first 2000 molecules of medchem99 against the same database (33046 structures) using 512bit fingerprints drops from 25.7s to 16.7s, and on 2048bit fingerprints from 38.7s to 23.3s.

7. "Triage" finite state machine matching

For example, consider a substructure search of benzene. Any canonical SMILES in the database that contains the substring 'c1ccccc1' must contain a benzene ring as a substructure. And of the 126705 structures in NCI95, UNIX grep command can rapidly identify 14377 such molecules. The advantage is that using a technique called finite state machines (FSMs), or more technically deterministic finite automata (DFAs), searching for substrings can be done incredibly rapidly. DFAs don't perform any backtracking, and can process each character using only five or six machine instructions. In theory, a 360MHz processor can examine 60 million characters per second, or approximately 2 million structures a second, or NCI95 in under a tenth of a second. For the benzene example above, this 0.1s avoids processing 14377 molecules!

The implementation used in DayCart uses a Aho & Corasick's algorithm that converts a trie of words to match into a DFA. This O(n) method inserts the necessary back edges to avoid backtracking much like the Knuth-Morris-Pratt string matching algorithm does for a single string.

In the actual implementation the DFA is represented as a 2-dimensional matrix. The first ordinate being the current state, the second being the next input character and the value in the matrix the sucessor state. The actual matching loop is shown below.

int FSMMatch( FSMType *fsm, char *str )
{
    register char *ptr;
    register int state;
    register int ch;

    ptr = str;
    state = 2;
    while( *ptr ) {
        ch = fsm->map[*ptr++];
        state = fsm->tran[state][ch];
    }
    /* return state == 1; */
    return fsm->tran[state][0] == 1;
}

Although the details of FSM construction are omitted, the set of substrings to match and modifications to the classic construction are detailed below.

1. The triage algorithm is currently only applied to single component (connected) molecules without stereochemistry. Hence if the canonical SMILES of the query contains a '.' (indicating multiple parts), a '>' (indicating a reaction), or '@', '/' or '\' (indicating stereo chemistry) the triage is skipped. Isotopes are correctly handled. Future implementations may support stereochemistry [with more advanced SMILES generation] and searches of reactant/product sides of the reaction [by simple manipulation of the DFA].

2. Instead of searching for just the canonical SMILES, triage currently attempts to generate all possible SMILES of a query molecule. However, to once again avoid combinatorial explosions, the number of SMILES is currently capped at 100. For example, methylbenzene, currently generates the 12 permutations below:

   Cc1ccccc1       c1(C)ccccc1     c1(ccccc1)C
   c1c(C)cccc1     c1c(cccc1)C     c1ccccc1C
   c1cc(C)ccc1     c1cc(ccc1)C     c1cccc(C)c1
   c1cccc(c1)C     c1ccc(C)cc1     c1ccc(cc1)C
   

3. To reduce the number of SMILES permutations, the algorithm takes advantage of a property of canonical SMILES that higher bond orders are prefered over lower. This reduces the number of permutations of a carboxy query to 3, 'O=CO', 'C(=O)O', 'O=CO'.

4. To avoid semantic differences, single bonds between aromatic atoms, such as in biphenyl, are explicitly written as '-'. This preserves the semantics, but for v471 only matches where the bridge bond is in a ring. Hence, a biphenyl search does not triage itself, but does find 146 superstructures in NCI95.

5. Multiple versions of each SMILES permutation are inserted into the trie by sequentially incrementing ring closure digits, whilst the highest ring closure index is less then %10. The nine benzene cases, increase the number of matches in NCI95 from 14377 to 40029.

6. Additional states are required at the end of SMILES ending in the characters 'C' and 'B' to avoid confusion with 'Cl' and 'Br', i.e. 'CC' is not a substructure of 'CCl'! These transitions are easy enough to add to the DFA, as an extra state where any input character other than the appropriate 'l' or 'r' leads to the sucess state. [And for the very observant explains the commented termination condition in the code above].

7. SMILES that are only a single character long, i.e. 'C', 'N', 'O', 'S', 'P', 'F', 'I' and 'B', require special attention to avoid finding a match inside square brackets, i.e. 'P' is not a substructure of '[Pt]'. This is handled by adding an extra state such that a '[' at the initial state goes to the new state, a ']' at the new state returns to the initial state, and every other character loops in the new state. Note that the query 'Cl' matches '[Cl-]', but the query 'I' does not match '[I-]'.

8. An easy to implement optimization to the DFA is to silently skip open parenthesis at states where they are not otherwise accepted. This doesn't increase the number of states, but greatly increases coverage. So that the query 'CC' will additionally match 'C(C'. With the benzene example, the number of NCI95 matches increases from 40029 to 50893.

9. The DFA is modified to enter the fail state upon encountering either a space, a tab or the string '%99'. The space and tab delimters are only useful in processing files, whilst the '%99' avoids the strange behaviour of wrap-around ring closure indices.

Performance results for a set of queries are shown in the table below:

8. Bibliography

1. A.V. Aho and M.J. Corasick, "Efficient String Matching: An Aid to Bibliographic Search", Communications of the ACM, Vol. 18, No. 6, pp. 333-340, 1975.
2. D. Keppel, S.J. Eggers and R.R. Henry, "Evaluating Runtime-Compiled Value Specific Optimizations", Technical Report UWCSE 93-11-02, Department of Computer Science and Engineering, University of Washington, 1993.
3. D.E. Knuth, "The Art of Computer Programming: Vol. 3, Sorting and Searching", 2nd Edition, Addison-Welsey Publishers, 1998.
4. D.E. Knuth, J.H. Morris and V.R. Pratt, "Fast Pattern Matching in Strings", Siam Journal of Computing, Vol. 6, No. 2, pp. 323-350, 1977.