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

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

Decision Tree by Scenario

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

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

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

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