Building a 3D Pharmacophore Model from PDB Data: A Free Python Workflow

I used to wonder how computational chemists extract an actionable drug-design blueprint directly from raw biological data. A protein-ligand crystal structure is just a list of atoms and coordinates — how do you get from that to a reusable query you can throw at millions of candidate molecules?

The answer is a 3D pharmacophore: an abstract spatial map of the chemical features a molecule must present to bind a target — hydrogen-bond donors and acceptors, aromatic rings, hydrophobic contacts — and the geometric relationships between them. This post walks through an example Python pipeline that builds a ligand-based pharmacophore from PDB data, end to end, using only free, open-source tools.

The whole idea rests on a simple insight: if many different ligands all bind the same pocket, the chemical features they share in 3D space are probably the ones that matter. So the strategy is to collect lots of bound ligands, physically overlay them, and find the features they agree on.

1. Mining the Protein Data Bank

The pipeline starts by asking the RCSB PDB for every structure of our target protein. Rather than scrape the website, we hit the GraphQL search API and let biotite fetch the structures.

The key here is strict quality control. We don’t want every structure — we want crystallographically reliable ones with a genuine drug-like ligand bound:

• resolved by X-ray crystallography,
• resolution ≤ 3.0 Å (anything coarser and atom positions get unreliable),
• a sensible chain configuration,
• and the presence of a non-polymer ligand > 100 Da (filtering out ions, water, and buffer components).

import requests
import biotite.database.rcsb as rcsb

UNIPROT_ID = "P00533"  # e.g. EGFR

# Build a composite query against the RCSB search API
query = rcsb.CompositeQuery(
    [
        rcsb.FieldQuery(
            "rcsb_polymer_entity_container_identifiers."
            "reference_sequence_identifiers.database_accession",
            exact_match=UNIPROT_ID,
        ),
        rcsb.FieldQuery("exptl.method", exact_match="X-RAY DIFFRACTION"),
        rcsb.FieldQuery(
            "rcsb_entry_info.resolution_combined", less_or_equal=3.0
        ),
        rcsb.FieldQuery(
            "chem_comp.formula_weight", greater=100  # drug-like ligand
        ),
    ],
    operator="and",
)

pdb_ids = rcsb.search(query)
print(f"{len(pdb_ids)} structures passed QC: {pdb_ids[:10]} ...")

The output is a curated list of PDB IDs — the raw material for everything that follows.

2. Extraction and 3D alignment

With the structures in hand, MDAnalysis parses each file and splits it into its components: the protein on one side, the ligand on the other.

The crucial step is 3D superposition. Each crystal structure lives in its own arbitrary coordinate frame, so the ligands are scattered across space and cannot be compared directly. We fix this by aligning every protein’s C-alpha backbone onto a single high-resolution reference structure. Because the binding pockets are rigid relative to the fold, aligning the proteins drags all the bound ligands into a common coordinate space — physically overlaid in the same pocket.

import MDAnalysis as mda
from MDAnalysis.analysis import align

ref = mda.Universe("reference.pdb")          # high-resolution template
aligned_ligands = []

for pdb_id in pdb_ids:
    mobile = mda.Universe(f"{pdb_id}.pdb")

    # superpose Cα backbones onto the reference
    align.alignto(
        mobile, ref,
        select="name CA",
        match_atoms=False,
    )

    # the ligand now sits in the reference frame
    ligand = mobile.select_atoms("not protein and not resname HOH")
    aligned_ligands.append(ligand.positions.copy())

After this stage, every ligand’s atomic coordinates are expressed in the same frame — ready to be compared, clustered, and merged.

3. Spatial clustering with DBSCAN

There’s a catch. Ligands don’t always bind the same place. A protein may have an orthosteric site (the main active site) and one or more allosteric sites, and our QC filter happily collected ligands from all of them. Overlaid, they form distinct clouds of atoms in space.

To isolate the pocket we actually care about, the pipeline runs DBSCAN — a density-based clustering algorithm — on the pooled 3D coordinates of all aligned ligand atoms. DBSCAN is ideal here because it finds clusters of arbitrary shape and labels sparse outliers as noise, without us having to specify the number of sites in advance.

import numpy as np
from sklearn.cluster import DBSCAN

# pool every atom from every aligned ligand
all_coords = np.vstack(aligned_ligands)

db = DBSCAN(eps=4.0, min_samples=10).fit(all_coords)
labels = db.labels_

# keep the densest cluster — our target pocket
valid = labels[labels != -1]
target_label = np.bincount(valid).argmax()
pocket_mask = labels == target_label

print(f"Target pocket holds {pocket_mask.sum()} atoms "
      f"across {len(set(valid))} detected site(s)")

We keep only the ligands sitting in the densest cluster — the orthosteric pocket — and move on with those.

4. Feature extraction and bond correction

Now a notorious gotcha. PDB files do not reliably store bond orders. An aromatic ring, a carbonyl, a carboxylate — they all come back as ambiguous single bonds. Feed that straight into a cheminformatics toolkit and your chemistry is silently wrong, which poisons every feature you extract afterward.

The fix is to treat the trustworthy 2D chemistry (a SMILES string from the PDB’s chemical component dictionary) as ground truth, and transfer its bond orders onto the 3D coordinates. RDKit does this with AssignBondOrdersFromTemplate, with rdFMCS available to reconcile partial or mismatched atoms via a maximum-common-substructure match.

Once the molecule is sanitized, RDKit’s feature factory sweeps the 3D structure and tags every pharmacophoric point — donors, acceptors, aromatics, hydrophobes — with its type and coordinates.

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import ChemicalFeatures
from rdkit import RDConfig
import os

# 1. fix bond orders using the 2D SMILES as a template
template = Chem.MolFromSmiles(ligand_smiles)          # correct chemistry
mol_3d   = Chem.MolFromPDBFile("ligand.pdb", removeHs=False)
mol      = AllChem.AssignBondOrdersFromTemplate(template, mol_3d)
Chem.SanitizeMol(mol)

# 2. map pharmacophoric features in 3D
fdef = os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")
factory = ChemicalFeatures.BuildFeatureFactory(fdef)

features = []
conf = mol.GetConformer()
for f in factory.GetFeaturesForMol(mol):
    pos = f.GetPos()
    features.append((f.GetFamily(), (pos.x, pos.y, pos.z)))

# e.g. ('Donor', (12.3, 4.1, -2.0)), ('Aromatic', (...)), ...

Repeat this for every ligand in the target cluster and you accumulate a big cloud of typed feature points — many noisy, many redundant, but collectively encoding what binding looks like for this target.

5. Generating the consensus model with k-means

A cloud of thousands of individual features isn’t a model — it’s the raw data for one. The final step distills it into a clean, minimal ensemble that represents the binding mode shared by the majority of active ligands.

The pipeline groups the pooled features by type (all acceptors together, all donors together, and so on) and runs k-means on the 3D coordinates within each group. Each resulting centroid is a candidate consensus feature. Centroids supported by only a handful of ligands are discarded as low-frequency noise; the ones the population agrees on survive.

import numpy as np
from collections import defaultdict
from sklearn.cluster import KMeans

# group every feature point by family
by_type = defaultdict(list)
for family, coord in all_features:
    by_type[family].append(coord)

n_ligands = len(target_cluster_ligands)
consensus = []

for family, coords in by_type.items():
    coords = np.array(coords)
    k = max(1, round(len(coords) / n_ligands))   # ~1 site per ligand
    km = KMeans(n_clusters=k, n_init="auto").fit(coords)

    for c in range(k):
        members = (km.labels_ == c).sum()
        # keep only features shared by the majority
        if members >= 0.5 * n_ligands:
            consensus.append((family, km.cluster_centers_[c]))

for family, center in consensus:
    print(f"{family:10s} at {np.round(center, 2)}")

What comes out is a short, interpretable list: a handful of typed features with precise 3D coordinates and well-defined distances between them.

The payoff

The final output is a compact, accurate 3D spatial map of the essential binding features for the target — derived entirely from public structural data, with no licensed software anywhere in the chain.

From here it’s ready to deploy in a virtual screening campaign: encode the consensus features and their inter-feature distances as a geometric query, then scan large compound libraries (ZINC, Enamine REAL, your in-house collection) for novel scaffolds that match the same spatial pattern — fresh chemistry that satisfies the same binding logic as the known actives. 🚀

The beauty of the approach is that every stage is a well-supported open-source library doing what it does best:

No proprietary suite required — just Python and a curated stack of free libraries, turning raw PDB coordinates into a drug-design blueprint.