Bio.PDB

Core Principle: Bio.PDB parses macromolecular structures into a hierarchical Structure-Model-Chain-Residue-Atom (SMCRA) object model for geometric analysis and spatial measurement.

The SMCRA hierarchy

Structure
 └─ Model          (NMR ensembles have many; X-ray usually one)
     └─ Chain      ('A', 'B', ...)
         └─ Residue  (het_flag, resseq, icode)
             └─ Atom  (name, coord, bfactor, occupancy, element)

Use mmCIF for anything modern. The legacy PDB format cannot represent structures with more than 62 chains or 99,999 atoms, which excludes large complexes and most cryo-EM structures.

Iterating

Residue identifiers

Filter out waters and heteroatoms before geometry:

The insertion code (icode) exists because antibody numbering schemes insert residues like 100A, 100B. Ignoring it silently merges distinct residues.

Geometry

import numpy as np
from Bio.PDB import calc_angle, calc_dihedral

d = atom1 - atom2                    # Atom subtraction gives DISTANCE in Å
calc_angle(a1.get_vector(), a2.get_vector(), a3.get_vector())
calc_dihedral(a1.get_vector(), a2.get_vector(), a3.get_vector(), a4.get_vector())

from Bio.PDB.Polypeptide import PPBuilder
ppb = PPBuilder()
for pp in ppb.build_peptides(chain):
    print(pp.get_sequence())                 # → a Seq object
    phi_psi = pp.get_phi_psi_list()          # backbone dihedrals, in radians

atom1 - atom2 returning a distance is a nice piece of operator overloading, and easy to misread as vector subtraction. For the actual displacement vector use atom1.coord - atom2.coord.

Neighbour search

from Bio.PDB import NeighborSearch

atoms = [a for a in structure.get_atoms()]
ns = NeighborSearch(atoms)

ns.search(center_coord, radius=5.0, level="R")     # residues within 5 Å
ns.search_all(4.0, level="R")                       # all residue pairs within 4 Å

NeighborSearch builds a KD-tree, so it is O(log n) per query rather than O(n). For contact maps and binding-site analysis, use it: a nested loop over 5,000 atoms is 25 million distance calculations.

Superposition and RMSD

from Bio.PDB import Superimposer

sup = Superimposer()
sup.set_atoms(fixed_atoms, moving_atoms)     # must be the same length, in order
print(sup.rms)                                # RMSD in Å
sup.apply(moving_structure.get_atoms())       # transform in place

The two atom lists must correspond one-to-one. For structures with different residue numbering, align the sequences first and use the alignment to pair the residues.

Writing

Biopython 1.86 made PDBIO enforce the wwPDB B-factor field width (max 6 characters). If you have written custom values into bfactor: a common trick for colouring a structure by conservation or by pLDDT: values outside the legal range now raise instead of writing a malformed file.

Downloading structures

For AlphaFold models, fetch directly from the AlphaFold DB: https://alphafold.ebi.ac.uk/files/AF-P04637-F1-model_v4.cif

In AlphaFold structures the B-factor column holds pLDDT (0–100 confidence), not a temperature factor. Regions below 50 are usually disordered and should not be interpreted structurally.

Bioinformatics examples

# where does a missense variant sit in the structure?
res = structure[0]["A"][175]
print(res.get_resname())                       # confirm it matches the expected ref AA

ns = NeighborSearch(list(structure.get_atoms()))
contacts = ns.search(res["CA"].coord, 8.0, level="R")
print(f"{len(contacts)} residues within 8 Å")

# is it buried or exposed? (needs DSSP installed)
from Bio.PDB.DSSP import DSSP
dssp = DSSP(structure[0], "1tup.pdb")
key = ("A", (" ", 175, " "))
ss, rsa = dssp[key][2], dssp[key][3]
print(f"secondary structure {ss}, relative accessibility {rsa:.2f}")

# distance to a bound ligand or DNA
zn = [a for a in structure.get_atoms() if a.element == "ZN"]
print(min(res["CA"] - a for a in zn))

# per-residue B-factor / pLDDT profile
bf = pd.DataFrame([
    {"resnum": r.id[1], "resname": r.get_resname(),
     "bfactor": np.mean([a.bfactor for a in r])}
    for r in structure[0]["A"] if r.id[0] == " "
])
sns.lineplot(data=bf, x="resnum", y="bfactor")

# contact map
ca = [r["CA"] for r in structure[0]["A"] if r.id[0] == " " and "CA" in r]
coords = np.array([a.coord for a in ca])
dist = np.linalg.norm(coords[:, None, :] - coords[None, :, :], axis=-1)
sns.heatmap(dist < 8, cmap="Greys", square=True, cbar=False)

# map a variant onto conservation, then colour the structure
for r in structure[0]["A"]:
    for a in r:
        a.bfactor = conservation.get(r.id[1], 0.0)     # keep within field width!
io.save("coloured.pdb")

The contact-map line uses Broadcasting to compute all pairwise distances in one expression: (n,1,3) - (1,n,3) → (n,n,3). Fine for a few hundred residues; use scipy.spatial.distance for anything larger.

Common mistakes

• Not filtering waters and heteroatoms.
• Ignoring insertion codes, merging distinct residues.
• Assuming one model. NMR ensembles have 20+.
• atom1 - atom2 misread as vector subtraction. It is a distance.
• Missing residues. Crystal structures have disordered gaps; residue numbering is not contiguous and does not necessarily match UniProt numbering.
• PDB numbering ≠ UniProt numbering. Map through SIFTS before placing a variant.
• Using the legacy PDB format for large complexes.
• Treating AlphaFold B-factors as temperature factors. They are pLDDT.
• Nested loops for contacts instead of NeighborSearch.
• Out-of-range B-factor values now raising in 1.86+.

See also

Seq · SeqRecord · Bio.Entrez · Linear Algebra with NumPy · Broadcasting · Applied - Heatmaps and Clustermaps

Exercise

Biopython is not installed here, so instead of PDBParser you will read the coordinate lines with fixed-column slicing: which is what the legacy PDB format actually is, and why it cannot hold more than 99,999 atoms. The columns you need, 0-based and half-open, are: record name [0:6], atom name [12:16], altLoc [16], residue name [17:20], chain [21], residue number [22:26], insertion code [26], and the coordinates at [30:38], [38:46] and [46:54].

Write trace(chain_id, text) returning the report lines for one chain. Consider only lines whose record name is ATOM, so waters and other heteroatoms in HETATM records are excluded, and only atoms whose altLoc is either blank or A, so an alternate conformer is not counted twice. A residue is identified by its number together with its insertion code; treating the number alone as the identifier merges residues such as 100 and 100A. Print residues=<n> counting distinct residues in first-seen order, then ca_atoms=<m> counting the first CA atom of each residue, then trace=<x> for the summed distance between consecutive CA atoms and span=<x> for the straight-line distance from the first CA to the last, both rounded to two decimals. With fewer than two CA atoms both distances are 0.00. The chain id arrives on the first line of input and the file follows.