Lesson 13 of 13 · 12 min

Capstone: An Automated Analysis Pipeline and Where to Go Next

This capstone consolidates the whole course into one reproducible command-line pipeline. Given a single PDB accession it fetches the coordinates, parses them, assigns secondary structure with DSSP, computes a residue contact count, maps the ligand binding site, renders a publication-style PyMOL figure, and emits a machine-readable metrics report. Everything runs from one script so a collaborator can reproduce your numbers and your figure from the accession alone.

Pipeline Overview

The stages are fetch, parse, DSSP, contacts, binding site, figure, and report. Reproducibility depends on pinning the toolchain: mkdssp version 4 reads mmCIF natively and Biopython 1.81 or newer parses its output, while PyMOL open-source gives you headless ray tracing that is deterministic for a fixed version and build. Pin these in an environment file and write all artifacts into one working directory: the numeric JSON report is then bit-for-bit reproducible, and the figure is visually reproducible on the same pinned stack, though ray-traced PNGs are not guaranteed byte-identical across different machines or thread counts.

bash

# Pinned, reproducible environment (conda-forge + bioconda)
conda create -n struct -c conda-forge -c bioconda \
  python=3.11 biopython=1.83 dssp=4 pymol-open-source -y
conda activate struct

# mkdssp and pymol must be discoverable on PATH
mkdssp --version
pymol -cq -d "print('pymol ok')"

python capstone_pipeline.py

The Capstone Script

python

#!/usr/bin/env python3
"""Fetch -> parse -> DSSP -> CB contacts -> binding site -> PyMOL figure -> report."""
import json
import subprocess
import warnings
from pathlib import Path

from Bio.PDB import MMCIFParser, PDBList, DSSP, NeighborSearch, Selection
from Bio.PDB.PDBExceptions import PDBConstructionWarning

warnings.simplefilter("ignore", PDBConstructionWarning)

PDB_ID = "3ptb"          # bovine beta-trypsin
LIGAND = "BEN"           # benzamidine, bound in the S1 pocket
CHAIN = "A"
WORKDIR = Path("work")
WORKDIR.mkdir(exist_ok=True)


def fetch(pdb_id):
    return PDBList().retrieve_pdb_file(
        pdb_id, pdir=str(WORKDIR), file_format="mmCif"
    )


def secondary_structure(model, cif_path):
    dssp = DSSP(model, cif_path, dssp="mkdssp")
    codes = [dssp[k][2] for k in dssp.keys()]   # index 2 = SS code
    n = len(codes)
    helix = sum(c in ("H", "G", "I") for c in codes)
    sheet = sum(c in ("E", "B") for c in codes)
    return {
        "residues_scored": n,
        "helix_fraction": round(helix / n, 3),
        "sheet_fraction": round(sheet / n, 3),
    }


def contact_stats(chain, cutoff=8.0, seq_sep=3):
    reps = []
    for res in chain:
        if res.id[0] != " ":                    # drop HETATM and water
            continue
        atom = res["CB"] if "CB" in res else (res["CA"] if "CA" in res else None)
        if atom is not None:
            reps.append((res.id[1], atom))
    contacts = 0
    for i in range(len(reps)):
        for j in range(i + 1, len(reps)):
            if abs(reps[i][0] - reps[j][0]) <= seq_sep:
                continue
            if reps[i][1] - reps[j][1] < cutoff:   # Biopython atom subtraction = distance
                contacts += 1
    return {"residues": len(reps), "cb_contacts_8A": contacts}


def binding_site(structure, ligand_resname, cutoff=4.5):
    atoms = Selection.unfold_entities(structure, "A")
    ns = NeighborSearch(atoms)
    site = {}
    for a in atoms:
        if a.get_parent().resname != ligand_resname:
            continue
        for res in ns.search(a.coord, cutoff, level="R"):
            if res.id[0] != " ":                # standard residues only
                continue
            chain_id = res.get_parent().id
            site[(chain_id, res.id[1])] = res.resname
    labels = [f"{c}:{site[(c, n)]}{n}" for (c, n) in sorted(site)]
    return {
        "ligand": ligand_resname,
        "cutoff_A": cutoff,
        "n_site_residues": len(labels),
        "site_residues": labels,
    }


def render_figure(cif_path, ligand):
    png = WORKDIR / f"{PDB_ID}_site.png"
    pml = f"""
load {cif_path}, cplx
hide everything
bg_color white
show cartoon, polymer
color grey80, polymer
select site, byres (polymer within 4.5 of resn {ligand})
show sticks, site and not (name C+N+O)
show sticks, resn {ligand}
color marine, site
color yellow, resn {ligand}
set ray_shadows, 0
orient resn {ligand}
zoom resn {ligand}, 8
ray 1200, 900
png {png}, dpi=300
"""
    script = WORKDIR / "figure.pml"
    script.write_text(pml)
    subprocess.run(["pymol", "-cq", str(script)], check=True)
    return png


def main():
    cif_path = fetch(PDB_ID)
    structure = MMCIFParser(QUIET=True).get_structure(PDB_ID, cif_path)
    model = structure[0]

    report = {
        "pdb_id": PDB_ID,
        "secondary_structure": secondary_structure(model, cif_path),
        "contacts": contact_stats(model[CHAIN]),
        "binding_site": binding_site(structure, LIGAND),
    }
    report["figure"] = str(render_figure(cif_path, LIGAND))

    (WORKDIR / f"{PDB_ID}_report.json").write_text(json.dumps(report, indent=2))
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()

{
  "pdb_id": "3ptb",
  "secondary_structure": {
    "residues_scored": 223,
    "helix_fraction": 0.099,
    "sheet_fraction": 0.323
  },
  "contacts": {
    "residues": 223,
    "cb_contacts_8A": 1387
  },
  "binding_site": {
    "ligand": "BEN",
    "cutoff_A": 4.5,
    "n_site_residues": 9,
    "site_residues": [
      "A:ASP189",
      "A:SER190",
      "A:GLN192",
      "A:SER195",
      "A:TRP215",
      "A:GLY216",
      "A:GLY219",
      "A:GLY226",
      "A:VAL227"
    ]
  },
  "figure": "work/3ptb_site.png"
}

Try it yourself: Point the script at the EGFR kinase-inhibitor complex 3POZ (set PDB_ID to 3poz and LIGAND to 03P). Compare its helix and sheet fractions against trypsin, list the residues within 4.5 A of the inhibitor, and extend binding_site() to also record each site residue's relative solvent accessibility, which is index 3 of every dssp[key] tuple (thread the dssp object into the function so it has access). You should see a far higher helix fraction than the trypsin beta-barrel and hinge residues lining the ATP pocket.

Where To Go Next

A crystal structure is one static, low-energy snapshot. Molecular dynamics with OpenMM or GROMACS propagates the system under a force field to sample conformational fluctuations, side-chain rotamer exchange, and induced-fit motion that a static contact map cannot reveal. Molecular docking with AutoDock Vina predicts ligand poses in the pocket you just mapped; validate the protocol by redocking the crystallographic ligand and checking pose RMSD against the deposited coordinates before you trust any novel compound. For larger or lower-resolution assemblies, cryo-EM density maps from the EMDB let you fit and refine models in Coot, Phenix, or ChimeraX, judging quality by local resolution and map-to-model FSC rather than crystallographic R factors.

When no experimental structure exists, predicted models from AlphaFold or a local ColabFold run fill the gap, but read their confidence metrics first. pLDDT is a per-residue self-confidence on a 0 to 100 scale: above 90 is high accuracy, 70 to 90 is a reliable backbone with less certain side chains, and below 50 usually flags intrinsically disordered or unconstrained regions rather than a wrong fold. The Predicted Aligned Error matrix gives expected positional error between residue pairs, so high-error off-diagonal blocks mean two individually confident domains whose relative orientation is uncertain. Keep the caveats in view: pLDDT is the model grading itself, not experimental validation; predictions are typically apo, single-conformation, and omit ligands, cofactors, ions, and post-translational modifications; and a high-confidence fold can still be biologically wrong for allosteric, membrane, or complex-dependent states.

Structure Quality and Validation