Lesson 10 of 13 · 12 min

Analyzing Ligand Binding Sites

A ligand binding site is the shell of protein residues that line the cavity a bound small molecule occupies. Characterizing it programmatically turns a static complex into a quantitative description of shape, chemistry, and contacts you can compare across structures. In a PDB file the ligand is stored as HETATM records that share the flag with solvent, ions, and additives, so you must select only the ligand of interest and drop water before defining a pocket. This lesson isolates the HETATM ligand, expands a 5 angstrom residue shell in both PyMOL and Biopython, and catalogs the contacts, using EGFR kinase bound to erlotinib (ligand code AQ4) as the worked example.

Isolating the Ligand

text

# PyMOL command console
fetch 1m17, async=0

# HETATM ligand only: erlotinib carries the residue name AQ4
select ligand, resn AQ4
remove solvent

# byres expands to complete residues; keep only protein atoms
# with any atom within 5 A of the ligand
select pocket, byres (polymer within 5 of ligand)

# one line per pocket residue, plus a count
iterate pocket and name CA, print(f"{chain}/{resn}-{resi}")
print("pocket residues:", cmd.count_atoms("pocket and name CA"))

A/LEU-694
A/GLY-695
A/VAL-702
A/ALA-719
A/LYS-721
A/GLU-738
A/MET-742
A/LEU-764
A/ILE-765
A/THR-766
A/GLN-767
A/LEU-768
A/MET-769
A/PRO-770
A/PHE-771
A/GLY-772
A/CYS-773
A/LEU-820
A/THR-830
A/ASP-831
pocket residues: 20

Pocket Residues With Biopython

python

from Bio.PDB import PDBList, PDBParser, NeighborSearch, Selection

PDBList().retrieve_pdb_file("1m17", pdir=".", file_format="pdb")
structure = PDBParser(QUIET=True).get_structure("1m17", "pdb1m17.ent")
model = structure[0]

# Ligand atoms: HETATM residues whose hetfield is 'H_AQ4'
ligand_atoms = [a for a in model.get_atoms()
                if a.get_parent().id[0] == "H_AQ4"]

# NeighborSearch builds a KD-tree over every atom in the model
atoms = Selection.unfold_entities(model, "A")
ns = NeighborSearch(atoms)

# Any standard residue with an atom <= 5.0 A from any ligand atom
pocket = set()
for lig_atom in ligand_atoms:
    for res in ns.search(lig_atom.coord, 5.0, level="R"):
        if res.id[0] == " ":            # standard amino acid, not HETATM/water
            pocket.add(res)

for res in sorted(pocket, key=lambda r: r.id[1]):
    print(res.get_parent().id, res.resname, res.id[1])
print("pocket residues:", len(pocket))

A LEU 694
A GLY 695
A VAL 702
A ALA 719
A LYS 721
A GLU 738
A MET 742
A LEU 764
A ILE 765
A THR 766
A GLN 767
A LEU 768
A MET 769
A PRO 770
A PHE 771
A GLY 772
A CYS 773
A LEU 820
A THR 830
A ASP 831
pocket residues: 20

Characterizing Pocket Chemistry

python

import numpy as np
from collections import Counter

# Side-chain chemical classes
CLASS = {
    "ARG": "positive", "LYS": "positive", "HIS": "positive",
    "ASP": "negative", "GLU": "negative",
    "SER": "polar", "THR": "polar", "ASN": "polar", "GLN": "polar",
    "CYS": "polar", "TYR": "polar",
    "GLY": "nonpolar", "ALA": "nonpolar", "VAL": "nonpolar",
    "LEU": "nonpolar", "ILE": "nonpolar", "MET": "nonpolar",
    "PRO": "nonpolar", "PHE": "nonpolar", "TRP": "nonpolar",
}

lig_coords = np.array([a.coord for a in ligand_atoms])

rows = []
for res in sorted(pocket, key=lambda r: r.id[1]):
    res_coords = np.array([a.coord for a in res])
    # minimum ligand-to-residue heavy-atom distance
    dmin = np.linalg.norm(
        res_coords[:, None, :] - lig_coords[None, :, :], axis=-1
    ).min()
    rows.append((res.resname, res.id[1], CLASS[res.resname], dmin))

classes = Counter(c for _, _, c, _ in rows)
polar = classes["polar"] + classes["positive"] + classes["negative"]
frac_polar = polar / len(rows)

for resname, resi, klass, dmin in rows:
    print(f"{resname}{resi:<4} {klass:<9} {dmin:5.2f} A")
print(f"\ncomposition: {dict(classes)}")
print(f"polar fraction: {frac_polar:.2f}")

LEU694  nonpolar   3.27 A
GLY695  nonpolar   4.50 A
VAL702  nonpolar   4.19 A
ALA719  nonpolar   3.31 A
LYS721  positive   3.70 A
GLU738  negative   3.81 A
MET742  nonpolar   4.63 A
LEU764  nonpolar   3.43 A
ILE765  nonpolar   4.45 A
THR766  polar      3.41 A
GLN767  polar      3.15 A
LEU768  nonpolar   3.50 A
MET769  nonpolar   2.70 A
PRO770  nonpolar   3.11 A
PHE771  nonpolar   3.45 A
GLY772  nonpolar   3.24 A
CYS773  polar      4.73 A
LEU820  nonpolar   3.45 A
THR830  polar      3.65 A
ASP831  negative   3.11 A

composition: {'nonpolar': 13, 'positive': 1, 'negative': 2, 'polar': 4}
polar fraction: 0.35

The profile reads as a classic ATP pocket: thirteen nonpolar residues dominate, holding the polar fraction to 0.35, while the closest contact at 2.70 angstrom is the hinge residue Met769, whose backbone amide donates a hydrogen bond to the quinazoline N1 of erlotinib (labeled N2 in the AQ4 atom nomenclature). Thr766 is the gatekeeper, Asp831 is the DFG aspartate, and Glu738 is the alpha-C helix glutamate that salt-bridges the catalytic Lys721, so the handful of polar and charged residues sit exactly where selectivity and catalysis are decided. A contact near 2.7 to 3.2 angstrom to a nitrogen or oxygen flags a likely hydrogen bond, whereas 3.5 to 4.7 angstrom contacts are van der Waals packing.

Try it yourself: Re-run the Biopython analysis on a more polar site, such as 1STP (streptavidin bound to biotin, ligand BTN) or 3ERT (estrogen receptor with 4-hydroxytamoxifen, ligand OHT). Tighten the cutoff to 4.0 angstrom to keep only first-shell contacts, then compare the polar fraction against the hydrophobic ATP pocket above. Which residue makes the closest contact, and is the nearest atom on the backbone or the side chain?

Contacts, Interfaces, and Hydrogen Bonds Structural Superposition and Alignment