bio-substructure-search

What this skill does

## Version Compatibility

Reference examples tested with: RDKit 2024.09+. SMARTS dialect follows Daylight specification with RDKit extensions.

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show rdkit` then `help(rdkit.Chem.MolFromSmarts)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Substructure Search

Search molecular collections for structural patterns using SMARTS. The choice of SMARTS dialect, atom/bond matching mode, and structural-alert catalog determines whether the search is correctly capturing the intended chemistry. PAINS (Baell & Holloway 2010) is the most-cited but most-misunderstood filter -- it identifies patterns of assay interference, not "bad molecules". Knowing when to apply each catalog and how to interpret hits is essential.

For SMARTS-based reactions (transforming matched substructures), see `chemoinformatics/reaction-enumeration`. For 3D pharmacophore matching, see `chemoinformatics/pharmacophore-modeling`.

## SMARTS Grammar Essentials

| Token | Meaning | Example |
|-------|---------|---------|
| `[#6]` | Atom by atomic number | `[#6]` carbon (any hybridization) |
| `c` | Lowercase = aromatic | `c1ccccc1` benzene aromatic |
| `C` | Uppercase = aliphatic only | `C(=O)O` carboxylic acid carbon |
| `[CX4]` | Atom + connection count X | `[CX4]` sp3 carbon (4 connections) |
| `[CX3]=O` | Carbonyl (CX3 = sp2 with 3 bonds) | matches ketone, aldehyde, ester C |
| `[#6;R]` | Atom in ring | `[#6;R]` ring carbon |
| `[#6;!R]` | Atom not in ring | `[#6;!R]` acyclic carbon |
| `[#6;r6]` | Atom in 6-membered ring | `[#6;r6]` six-ring carbon |
| `[a]` | Any aromatic atom | `[a]` |
| `[!#1]` | Anything except H | `[!#1]` heavy atom |
| `[N;H2]` | N with exactly 2 H | `[NH2]` primary amine |
| `[N+]` | Positively charged N | `[N+](=O)[O-]` nitro |
| `[$(...)]` | Recursive SMARTS | `[$(c1ccccc1)]` aromatic 6-ring atom |
| `[c]([F,Cl,Br,I])` | OR within brackets | aryl halide |
| `~` | Any bond type | `c~c` any aromatic-aromatic bond |
| `@` | Aromatic bond | `c@c` |
| `-` | Single bond explicit | `C-C` |
| `=` | Double bond | `C=O` |
| `:` | Aromatic bond explicit | |

## Common SMARTS Patterns

| Pattern | SMARTS | Notes |
|---------|--------|-------|
| Hydroxyl (alcohol + phenol) | `[OX2H]` | OX2H avoids matching O- in OH- |
| Phenol only | `[OX2H][c]` | OH attached to aromatic carbon |
| Aliphatic OH only | `[OX2H][CX4]` | OH attached to sp3 C |
| Carboxylic acid | `[CX3](=O)[OX2H1]` | C(=O)OH |
| Carboxylate | `[CX3](=O)[O-]` | C(=O)O- (deprotonated) |
| Ester | `[CX3](=O)[OX2][!H]` | C(=O)O-R |
| Amide | `[CX3](=[OX1])[NX3]` | C(=O)N-R |
| Primary amine | `[NX3;H2]` | -NH2 |
| Secondary amine | `[NX3;H1]` | -NH-R |
| Tertiary amine | `[NX3;H0;!$(NC=O)]` | -NR2 (not amide N) |
| Quaternary amine | `[NX4+]` | -NR4+ |
| Nitro | `[N+](=O)[O-]` | -NO2 |
| Nitrile | `[CX2]#[NX1]` | -C#N |
| Sulfonamide | `[SX4](=[OX1])(=[OX1])[NX3]` | -S(=O)(=O)N |
| Aryl halide | `[c][F,Cl,Br,I]` | halogen on aromatic |
| Aliphatic halide | `[CX4][F,Cl,Br,I]` | halogen on sp3 C |
| Hydrogen bond donor | `[#7,#8;!H0]` | N or O with at least 1 H |
| Hydrogen bond acceptor | `[#7,#8;!$([NX3]([O-])=O);!$([N+]=O)]` | N/O excluding nitro |
| Michael acceptor | `[CX3]=[CX3][CX3]=O` | enone, acrylamide warhead |
| Aldehyde | `[CX3H1](=O)` | -CHO |
| Ketone | `[CX3](=O)[#6]` | -C(=O)R, both R = C |

## Basic Substructure Match

**Goal:** Test whether a molecule contains a SMARTS pattern and enumerate the matching atom indices.

**Approach:** Parse the molecule with `MolFromSmiles` and the pattern with `MolFromSmarts`, gate with `HasSubstructMatch`, then call `GetSubstructMatches` and map each atom index back to the molecule for inspection.

```python
from rdkit import Chem

mol = Chem.MolFromSmiles('c1ccc(O)cc1CCO')
pattern = Chem.MolFromSmarts('[OX2H]')

if mol.HasSubstructMatch(pattern):
    matches = mol.GetSubstructMatches(pattern)
    for match in matches:
        atoms = [mol.GetAtomWithIdx(i).GetSymbol() for i in match]
```

`HasSubstructMatch` returns bool, `GetSubstructMatches` returns tuple of tuples of atom indices.

## Recursive SMARTS (key for postdoc-grade patterns)

`[$(pattern)]` matches an atom that *also* matches the entire pattern starting from itself. Critical for context-aware matching.

```python
# Aromatic carbon attached to a carbonyl
pat = Chem.MolFromSmarts('[$(c[C](=O))]')

# Aniline-type N (aromatic carbon-N-H)
pat = Chem.MolFromSmarts('[$([NX3;H2][c])]')

# Hindered amine (N with 2 sp3 neighbors)
pat = Chem.MolFromSmarts('[$([NX3]([CX4])([CX4])[CX4])]')

# H-bond donor (per Lipinski, exclude quaternary)
hbd = Chem.MolFromSmarts('[#7,#8;!H0;!$([NX3+])]')

# H-bond acceptor (per Lipinski, exclude nitro / aniline)
hba = Chem.MolFromSmarts('[$([#7,#8;!H0]);!$([NX3+]=O);!$(N(=O)~O)]')
```

## Structural-Alert Filter Catalogs

| Filter | Origin | Patterns | Use case | Failure mode |
|--------|--------|----------|----------|--------------|
| PAINS_A | Baell & Holloway 2010 (low-quality assay hits) | 480 | Flag known pan-assay interferers | Many false positives in primary screens; legitimate medicines flagged |
| PAINS_B | Baell & Holloway 2010 | 280 | More aggressive PAINS | Similar |
| PAINS_C | Baell & Holloway 2010 | 240 | Most aggressive PAINS | Most permissive |
| BRENK | Brenk 2008 (DDS unsuitable) | 105 | Reactive / toxicity / undesirable | Useful for fragment / virtual library |
| NIH | NIH MLSMR | ~250 | Reactive groups, unstable | Legacy filter |
| ZINC | ZINC clean-leads | ~90 | Drug-like cleanup | Used for library standardization |
| Aldridge | Aldridge medchem rules | ~50 | medchem ugly substructures | Hand-curated |
| Glaxo / Eli Lilly | Vendor lists | varies | Internal "ugly" filters | Often unpublished |
| REOS | Walters & Murcko 2002 | property + structural | Drug-likeness combined filter | Hand-curated thresholds |

## When to Apply Each Filter

| Scenario | Catalog | Reason |
|----------|---------|--------|
| Hit validation from biochemical screen | PAINS_A | Identify assay-interference candidates |
| Library prep for HTS | PAINS_A + Brenk + ZINC | Remove clearly bad |
| Fragment library design | Brenk + ZINC | Remove reactive; PAINS less critical at fragments |
| Lead optimization | None mandatory | Filters can exclude valid leads |
| Natural product analog | None | Filters trained on synthetic chemistry |
| Covalent inhibitor design | Skip warhead filter | Warheads ARE the design |

**Critical:** Capuzzi et al. (2017) showed that 8% of FDA-approved drugs match a PAINS pattern. PAINS is a *flag for assay validation*, not a *killing filter*.

## PAINS Filter

**Goal:** Split a molecule list into PAINS-flagged and PAINS-clean sets using one or more PAINS catalog tiers.

**Approach:** Configure `FilterCatalogParams` with the requested catalog enums, build a `FilterCatalog` once, and for each molecule use `GetFirstMatch` to either bucket it as clean or record the matching pattern description.

```python
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams

def pains_filter(mols, catalogs=('PAINS_A',)):
    params = FilterCatalogParams()
    for cat in catalogs:
        params.AddCatalog(getattr(FilterCatalogParams.FilterCatalogs, cat))
    catalog = FilterCatalog(params)

    flagged = []
    clean = []
    for mol in mols:
        if mol is None:
            continue
        entry = catalog.GetFirstMatch(mol)
        if entry is None:
            clean.append(mol)
        else:
            flagged.append((mol, entry.GetDescription()))
    return clean, flagged
```

Available catalog names: `PAINS_A`, `PAINS_B`, `PAINS_C`, `PAINS` (all), `BRENK`, `NIH`, `ZINC`, `ALL`.

## Reaction-Reactive Group Filter (custom)

For HTS triage, filter electrophilic warheads (acrylamide, chloroacetamide, etc.) unless designing covalent inhib