GPTomics/bio-pharmacophore-modeling

Version Compatibility

Reference examples tested with: RDKit 2024.09+, Pharmit web service, and PLIP 2.4+ (interaction analysis). Verify the deployed Pharmit/Pharmer interface and query format before automation.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show rdkit then help(rdkit.Chem.Pharm3D) 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.

Pharmacophore Modeling

Build 3D pharmacophore queries that capture the essential interaction features of a ligand-target binding event. A pharmacophore is the spatial arrangement of pharmacophore features (donor, acceptor, hydrophobe, aromatic, charged) sufficient for activity, abstracted from any specific chemotype. Use pharmacophores for scaffold hopping, virtual-screening prefilters, and cross-target SAR transfer. Derive interaction features directly from a co-crystal when available, use apo2ph4 to derive models from an apo pocket (Heider et al. 2023), or align known actives for a ligand-based model. PharmacoForge generates candidate 3D pharmacophores conditioned on a protein pocket; those pharmacophores can then retrieve matching molecules from a library (Flynn et al. 2025).

For 2D scaffold-based searches, see chemoinformatics/scaffold-analysis. For 3D shape similarity, see chemoinformatics/shape-similarity. For protein-ligand interaction analysis, see chemoinformatics/virtual-screening.

Pharmacophore Feature Types

| Feature | Common shorthand | Definition | Geometric tolerance | |---------|------------|------------|----------------------| | H-bond donor | D | -OH, -NH | 1.0-1.5 Å | | H-bond acceptor | A | sp2 O / N (lone pair) | 1.0-1.5 Å | | Hydrophobe | H | sp3 C / aromatic ring centroid | 1.5-2.0 Å | | Aromatic ring | R | Aromatic ring centroid + normal | 1.0-1.5 Å | | Positive ionizable | P | -NH3+, -NR3+ | 1.0-1.5 Å | | Negative ionizable | N | -COO-, -SO3- | 1.0-1.5 Å | | Halogen | X | Cl, Br, I (halogen bond donor) | 1.0-1.5 Å | | Metal coordination | M | sp/sp2 N/O near metal | 0.5-1.0 Å |

Tolerances are pharmacophore-feature distance windows in the search. Tighter tolerances = fewer hits but more specific.

The ranges in this table are repository starting heuristics, not universal feature tolerances. Set final bounds from aligned-feature variability, coordinate uncertainty, and retrospective validation for the selected search engine.

The one-letter labels above are human-readable shorthand, not RDKit API codes. RDKit's shipped BaseFeatures.fdef uses family names such as Donor, Acceptor, Hydrophobe, Aromatic, PosIonizable, and NegIonizable. Its default feature definitions do not provide every halogen-bond or metal-coordination model; add and validate project-specific feature definitions when those interactions matter.

Method Taxonomy

| Method | Origin | Use case | Fails when | |--------|--------|----------|------------| | Ligand-based (LBP) | Catalyst, MOE, RDKit Pharm3D | Multiple actives, no crystal | <3 actives; flexible actives | | Receptor-based (RBP) | apo2ph4, LigandScout, PLIP | Co-crystal or a defined apo pocket | Uncertain pocket conformation | | Common pharmacophore | Validated alignment/feature-consensus workflow; RDKit can represent and query the resulting model | Consensus from active set | Diverse actives or uncertain bioactive conformers confound alignment | | Pocket-conditioned generation (PharmacoForge) | Flynn et al. 2025 | Generate candidate pharmacophores from a protein pocket | Does not directly generate molecules; pretrained model required | | Active learning pharmacophore | Catalyst variant | Iterative refinement | Custom; not standard |

Decision Tree by Scenario

| Scenario | Method | Tools | |----------|--------|-------| | Co-crystal structure available | Interaction-derived receptor model | PLIP or LigandScout + Pharmit | | Apo structure with a defined pocket | Apo receptor model | apo2ph4; export LigandScout PML | | Multiple active compounds, no crystal | Ligand-based common pharmacophore | Alignment plus consensus-feature derivation in validated custom or external tooling; RDKit Pharm3D can apply the resulting model | | Single active compound | Single-conformer pharmacophore | RDKit Pharm3D from bioactive conformer | | Scaffold hopping prospective | Receptor-based + shape filter | apo2ph4 or interaction-derived model + shape search | | Cross-target SAR transfer | Common pharmacophore across targets | Manual + LigandScout | | Generate pocket-conditioned pharmacophores | PharmacoForge | Diffusion model followed by library retrieval | | Library pre-filtering | Pharmacophore screen | Pharmit search |

Ligand-Based Pharmacophore (RDKit Pharm3D)

Goal: Derive a common pharmacophore from aligned bioactive conformers, then apply that established model to candidate molecules.

Approach: Consensus derivation is a separate modeling step: select or generate plausible bioactive conformers, align them using a documented method, identify conserved feature correspondences, and estimate distance bounds or tolerances. RDKit does not provide a single EmbedPharmacophore call that performs those steps. EmbedPharmacophore instead generates conformations of a molecule that satisfy an already defined pharmacophore.

from rdkit import Chem, Geometry
from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import EmbedLib, Pharmacophore
from rdkit.RDPaths import RDDataDir
import os

fdef_file = os.path.join(RDDataDir, 'BaseFeatures.fdef')
factory = ChemicalFeatures.BuildFeatureFactory(fdef_file)

# This is an already defined model. Coordinates and bounds must come from a
# validated consensus-derivation workflow or another justified source. RDKit
# requires FreeChemicalFeature objects, not feature-family strings.
query_features = [
    ChemicalFeatures.FreeChemicalFeature(
        'Aromatic', Geometry.Point3D(0.0, 0.0, 0.0)),
    ChemicalFeatures.FreeChemicalFeature(
        'Donor', Geometry.Point3D(4.0, 0.0, 0.0)),
]
pharmacophore = Pharmacophore.Pharmacophore(query_features)
pharmacophore.setLowerBound(0, 1, 3.5)
pharmacophore.setUpperBound(0, 1, 5.0)

target = Chem.AddHs(Chem.MolFromSmiles('c1ccc(cc1)CCN'))
can_match, feature_matches = EmbedLib.MatchPharmacophoreToMol(
    target, factory, pharmacophore)
if can_match:
    atom_match = tuple(tuple(matches[0].GetAtomIds())
                       for matches in feature_matches)
    _, embeddings, n_failed = EmbedLib.EmbedPharmacophore(
        target, atom_match, pharmacophore, randomSeed=23, silent=True)

BaseFeatures.fdef (RDKit-shipped) defines feature SMARTS and is a useful starting feature taxonomy. The code above demonstrates applying an existing two-feature model; it does not infer a consensus model from active compounds.

Receptor-Based Pharmacophore (apo2ph4 workflow)

Goal: Derive a pharmacophore from a protein binding-pocket structure without requiring a bound ligand.

Approach: Identify donor, acceptor, and hydrophobic hot spots from apo-pocket geometry, cluster them, and assemble candidate pharmacophores. Heider et al. describe apo2ph4 in J. Chem. Inf. Model. 63:101-110 (2023). Use the source release's documented scripts and environment rather than assuming a packaged apo2ph4 command: the published workflow writes LigandScout PML output, not a generic .ph4 file. Treat conversion to Pharmit, Pharmer, MOE, or Phase as a separate, explicitly validated step because pharmacophore formats are not interchangeable.

When a co-crystal ligand is available, derive pharmacophore directly from the ligand binding pose: each ligand feature in contact with a complementary protein residue is part of the pharmacophore.

from plip.basic import config
from plip.structure.preparation import PDBComplex

mol_complex = PDBComplex()
mol_complex.load_pdb('complex.pdb')
mol_complex.analyze()

for site in mol_complex.interaction_sets.values():
    for interaction in site.all_itypes:
        # Objects are interaction-class-specific. Inspect the documented fields
        # for HydrophobicContact, HydrogenBond, PiStacking, SaltBridge, etc.;
        # there is no universal `.type` or `.ligatom.coords` interface.
        interaction_class = type(interaction).__name__
        print(interaction_class, interaction)

PLIP exposes typed interaction records with class-specific ligand/protein atoms and coordinates. Map those records to pharmacophore features explicitly and retain the interaction class and source atom identifiers.

Pharmacophore Search (Pharmit / Pharmer)

For library screening, configure feature types, centers, radii, and optional shape constraints in Pharmit, or use a Pharmer database and query produced in the format required by the installed release. Do not pass LigandScout PML or a vendor .ph4 file directly unless the selected interface documents that import path. Pharmit reported searching millions of conformers in seconds to minutes; actual runtime depends on query selectivity, database size, and deployment (Sunseri & Koes 2016).

Pharmacophore Quality Validation

Evaluate a pharmacophore by:

1. Retrospective enrichment: a stated metric on target-relevant actives and inactives/decoys. DUD-E can provide a benchmark with known decoy-construction biases; COCONUT is a natural-products collection, not a target-specific active/decoy benchmark.
2. Geometric tightness: feature distance variance across actives
3. Selectivity: false positives in inactive set should be low
4. Specific consistency: pharmacophore matches each active's bioactive conformer

def pharmacophore_enrichment(query_pharmacophore, actives, inactives,
                             matches_pharmacophore):
    """Return active/inactive match-rate enrichment for a supplied matcher."""
    if not actives or not inactives:
        raise ValueError('actives and inactives must both be non-empty')
    n_active_match = sum(
        bool(matches_pharmacophore(mol, query_pharmacophore))
        for mol in actives)
    n_inactive_match = sum(
        bool(matches_pharmacophore(mol, query_pharmacophore))
        for mol in inactives)
    active_rate = n_active_match / len(actives)
    inactive_rate = n_inactive_match / len(inactives)
    return float('inf') if inactive_rate == 0 else active_rate / inactive_rate

For this repository, enrichment >=5x may be used as a starting triage heuristic only after the active/decoy construction and matching policy are documented. Report the full metric and uncertainty, and calibrate the acceptance threshold on the project dataset.

Pocket-Conditioned Pharmacophore Generation (PharmacoForge)

PharmacoForge (Flynn et al. 2025) applies a diffusion model to a protein pocket and generates candidate 3D pharmacophores. It does not directly generate molecular structures from an input pharmacophore. The validated workflow is:

1. Prepare the protein pocket in the representation required by the published PharmacoForge release.
2. Sample and rank pocket-conditioned pharmacophores.
3. Convert a selected pharmacophore into the query representation used by the search engine.
4. Retrieve matching, purchasable compounds and evaluate them with docking, strain, and physical-validity checks.

The paper compares pharmacophore and downstream retrieval performance with other pocket-based approaches; it does not support a drug-likeness or novelty comparison with REINVENT.

Pharmacophore vs Shape vs 2D Fingerprint

| Method | Captures | Best for | |--------|----------|----------| | ECFP4 Tanimoto | Local atom environments | Lead optimization (same series) | | FCFP4 Tanimoto | Pharmacophore-equivalent atoms | Loose similarity in series | | Shape similarity (ROCS) | 3D shape volume | Scaffold hopping by shape | | Pharmacophore | Discrete features in space | Scaffold hopping with feature specificity | | Combined (Tanimoto + shape) | Multi-objective | Production VS |

Pharmacophore is more interpretable than shape: a hit explains why it matched (donor at position X, hydrophobe at position Y).

Per-Tool Failure Modes

Ligand-based -- diverse actives confound

Trigger: Active set spans multiple scaffolds with different bound conformations.

Mechanism: No common pharmacophore exists; algorithm forces non-consensus features.

Symptom: Pharmacophore matches no actives in retrospective.

Fix: Cluster actives by scaffold first; derive per-cluster pharmacophore.

Receptor-based -- apo structure

Trigger: Protein in apo form (no bound ligand).

Mechanism: Side-chain rotamers differ between apo and holo; "binding site" geometry is wrong.

Symptom: Pharmacophore inferred from apo doesn't match holo experimental data.

Fix: Use AlphaFold3 / Boltz-1 to predict holo conformation; derive pharmacophore from predicted holo.

Pharmacophore -- single conformer bias

Trigger: Active aligned to its first generated conformer, not bioactive conformer.

Mechanism: Crystal structure not available; generated conformer may not be the bound one.

Symptom: Pharmacophore inconsistent across runs (different starting conformer chosen).

Fix: Use conformer ensemble; align all to common scaffold; choose conformer most consistent with other actives.

Tolerance too tight

Trigger: Default geometric tolerance < 0.5 Å.

Mechanism: Real bioactive conformers have flexibility; rigid pharmacophore filters most molecules out.

Symptom: Search returns zero hits.

Fix: Use tolerance 1.0-1.5 Å for drug-like; up to 2 Å for flexible peptide-like.

Pharmacophore search misses bioisostere

Trigger: Bioisostere replacement (e.g., -COOH replaced by tetrazole).

Mechanism: Tetrazole functions as acid bioisostere but RDKit features may not classify identically.

Symptom: Known bioisosteric active not found.

Fix: Use ChemAxon-style bioisosteric feature equivalence; or pharmacophore feature class expansion (acid generic vs -COOH specific).

PLIP -- water bridge absent from output

Trigger: Bridging water between ligand donor and protein acceptor.

Mechanism: PLIP can report water bridges, but the required crystallographic water must be present in the input and satisfy its geometric criteria.

Symptom: Pharmacophore missing critical H-bond feature.

Fix: Retain relevant crystallographic waters, inspect PLIP water-bridge output, and review borderline geometry manually.

Reconciliation: Ligand-Based vs Receptor-Based

| Aspect | Ligand-based | Receptor-based | |--------|--------------|----------------| | Data needed | Multiple actives with defensible conformers/alignment | A defined pocket, optionally with a co-crystal ligand | | Main bias | Known active chemotypes, conformer choice, and alignment | Pocket structure, protonation, retained waters, and interaction-detection/modeling rules | | Hit-set behavior | Depends on feature abstraction and tolerances | Depends on selected pocket interactions, excluded volumes, and tolerances | | Confidence evidence | Retrospective recovery across held-out actives/inactives | Recovery of known interaction geometry and retrospective or prospective validation |

Choose between ligand- and receptor-based models using the available structural/activity evidence and target-relevant validation. Neither approach is universally more reliable, diverse, or suitable for scaffold hopping.

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | Pharm3D.EmbedPharmacophore fails | Bounds matrix infeasible | Review/loosen justified bounds and, when more attempts are warranted, increase the documented count argument; inspect n_failed | | Pharmacophore matches everything | Too few features | Add features; tighten tolerances | | Pharmacophore matches nothing | Too many features or tight bounds | Reduce feature count; loosen tolerances | | BaseFeatures.fdef not found | RDKit installation issue | Check from rdkit.RDPaths import RDDataDir | | Pharmacophore-conformer mismatch | Wrong conformer used | Use bioactive conformer from crystal | | Pharmit search timeout | Library too large | Pre-filter by 2D fingerprint Tanimoto | | apo2ph4 PML has no useful model | No robust pocket hot spots at selected settings | Recheck pocket definition and documented thresholds; inspect alternative models |

References

• Wolber & Langer, J. Chem. Inf. Model. 45:160-169 (2005) -- LigandScout pharmacophores. https://doi.org/10.1021/ci049885e
• Heider et al., J. Chem. Inf. Model. 63:101-110 (2023; published online 2022) -- apo2ph4. https://doi.org/10.1021/acs.jcim.2c00814
• Flynn EL, Shah R, Dunn I, Aggarwal R, Koes DR, Front. Bioinform. 5:1628800 (2025) -- PharmacoForge. https://doi.org/10.3389/fbinf.2025.1628800
• RDKit, Chem.Pharm3D API documentation. https://www.rdkit.org/docs/source/rdkit.Chem.Pharm3D.html
• RDKit, EmbedPharmacophore API documentation -- embedding molecules against an existing pharmacophore. https://www.rdkit.org/docs/source/rdkit.Chem.Pharm3D.EmbedLib.html#rdkit.Chem.Pharm3D.EmbedLib.EmbedPharmacophore
• COCONUT, official resource -- open natural-products collection. https://coconut.naturalproducts.net/
• Adasme et al., Nucleic Acids Res. 49:W530-W534 (2021) -- PLIP interaction profiler. https://doi.org/10.1093/nar/gkab294
• Sunseri & Koes, Nucleic Acids Res. 44:W442-W448 (2016) -- Pharmit interactive search. https://doi.org/10.1093/nar/gkw287

Related Skills

• chemoinformatics/molecular-io - Parse molecules
• chemoinformatics/conformer-generation - Generate 3D for pharmacophore
• chemoinformatics/shape-similarity - 3D shape adjacent to pharmacophore
• chemoinformatics/virtual-screening - Pharmacophore as docking pre-filter
• chemoinformatics/scaffold-analysis - 2D scaffold-hopping context
• chemoinformatics/generative-design - Generate or optimize molecules after pharmacophore-based retrieval
• structural-biology/structure-io - PDB handling