This package accelerates the matching of a query molecule to a large batch of patterns, given as SMARTS strings. This is a common subroutine in template-based retrosynthesis, where one is interested in which of 10-100k reaction templates are applicable for any given molecule. The key insight is that we can first use cheap heuristics to rule out patterns that can never match to avoid solving an expensive subgraph isomorphism problem.

Installation

Install batchmatch using pip:

pip install batchmatch 

Usage

The batchmatch library has one core object PatternLibrary, which takes as input a list of SMARTS strings and allows users to quickly query molecules against them using its match() method:

import batchmatch
from rdkit import Chem

with open("smarts.txt", "r") as f:
    smarts = f.read().splitlines()
query = Chem.MolFromSmiles("Cn1c(=O)c2c(ncn2C)n(C)c1=O")

# Instantiate the pattern library 
# Importantly, this only has to be done once, as it computes some indices  
patterns = batchmatch.PatternLibrary(
    smarts=smarts,
    consolidate=True,  # default; whether to deduplicate SMARTS 
    filters=batchmatch.DEFAULT_FILTERS,  # default; the filter heuristics used 
    chiral=True,  # default; whether to use chirality when matching
)

# Match the query molecule
for i in patterns.match(query):
    print(f"Query matches {smarts[i]}")

How does it work?

Under the hood, PatternLibrary first consolidates (deduplicates) the input SMARTS. Then, it applies a series of heuristics (filters) to rule out blatantly unmatching patterns, which can be run in isolation with the could_match() method. Currently, batchmatch supports 4 filters:

1. RingFilter: rules out patterns with rings not present in the query.
2. AtomCountFilter: rules out patterns based on their atom type counts. For example, a pattern with three oxygens can never match a molecule with only one.
3. RDKitFilter: uses a modified version of the RDKit fingerprint. Specifically, we rule out a pattern if their fingerprint's on bits are not a subset of the query's.
4. PatternFilter: the same as above but using RDKit's PatternFingerprint, which was adapted from MHN-React.

Empirically, we find the Ring > Atom > RDKit filter combination (in that order) works best and is what we use by default. Note that the bulk of the computation used in the filters (namely, computing statistics about the patterns) is done and cached on instantiation of PatternLibrary. Finally, we run an exhaustive substructure match with the remaining patterns to yield the final results.

Limitations

While batchmatch is well-tested on templates from AIZynthFinder and queries from ChEMBL, it is difficult to anticipate every use case due to the expansiveness of SMARTS. The strategies used herein are highly tailored towards SMARTS extracted from reaction templates, which only use a narrow subset of the specification language. I would always recommend running some spot checks and benchmarks on your data to see if batchmatch works well. For convenience,

batchmatch.PatternLibrary(smarts=smarts, consolidate=False, filters=[])

performs the brute-force strategy of matching against every SMARTS with no extra tricks!

Benchmarks

Below, I have run an informal benchmarks (test/benchmark.py) on my laptop (Intel Core i7-12800H processor), using the AIZynthFinder template set and 1k random molecules from ZINC. For each combination of 0-4 filters, I report the best permutation below in terms of match time (s) per molecule. I also consider the precision, the number of true matches divided by number of exhaustive substructure matches. An ideal filter (Ideal) would filter out exactly the true matches and therefore have a precision of 1, but not a match time of 0 because one still has to perform the matching routine.

Overall, we find that our default setting is almost 20x faster than a naive implementation and over 80x as precise! Compared to the Pattern strategy in MHN-React, we are 2x faster. However, there is still a lot of potential for improvement with a 15x slowdown and 5x decreased precision compared to the theoretical optimum. Likely, porting our filters to C++ will yield easy speed improvements, though I lack the knowledge to do this currently.