Lesson 6 of 13 · 11 min

Programmatic Access with Biopython Bio.PDB

Biopython's Bio.PDB module parses macromolecular structures into the SMCRA object model, a strict Structure over Model over Chain over Residue over Atom containment hierarchy. Use MMCIFParser for mmCIF files, the primary format of the modern PDB archive, and PDBParser for legacy fixed-column .pdb files; both build an identical in-memory tree, so downstream code is format-agnostic. Every container behaves like an ordered dictionary keyed by child id, which is why indexing such as structure[0]['A'] resolves a model and then a chain.

Parsing Into SMCRA

python

from Bio.PDB import MMCIFParser, PDBParser

cif_parser = MMCIFParser(QUIET=True)
structure = cif_parser.get_structure("1ake", "1ake.cif")

pdb_parser = PDBParser(QUIET=True)
structure_pdb = pdb_parser.get_structure("1ake", "1ake.pdb")

print(structure.id, "->", "models", [model.id for model in structure])
model = structure[0]
print("chains:", [chain.id for chain in model])
print("same tree from PDB:", [c.id for c in structure_pdb[0]])

1ake -> models [0]
chains: ['A', 'B']
same tree from PDB: ['A', 'B']

Iterating The Hierarchy

python

model = structure[0]
chain = model["A"]

for residue in chain:
    hetflag, resseq, icode = residue.get_id()
    if hetflag != " ":          # " " polymer, "W" water, "H_*" hetero ligand
        continue
    print("residue:", residue.get_resname(), residue.get_full_id())
    for atom in residue:
        print("  ", atom.get_name(), atom.get_coord(), atom.get_coord().dtype)
    break

residue: MET ('1ake', 0, 'A', (' ', 1, ' '))
   N [26.981 53.977 40.085] float32
   CA [26.091 52.849 39.889] float32
   C [26.679 52.163 38.675] float32
   O [27.02  52.865 37.715] float32
   CB [24.677 53.31  39.58 ] float32
   CG [23.624 52.189 39.442] float32
   SD [21.917 52.816 39.301] float32
   CE [21.93  53.926 37.91 ] float32

Flatten And Handle Disorder

python

from Bio.PDB import Selection

atoms = Selection.unfold_entities(structure, "A")     # levels: A R C M S
residues = Selection.unfold_entities(structure, "R")
print("atoms:", len(atoms), "residues:", len(residues))

for atom in atoms:
    if atom.is_disordered():
        res = atom.get_parent()
        for alt in atom.disordered_get_list():
            print(res.get_resname(), res.get_id()[1], atom.get_name(),
                  "altloc", alt.get_altloc(), "occ", alt.get_occupancy(),
                  alt.get_coord())
        break

atoms: 3804 residues: 808
ARG 167 CD altloc A occ 0.5 [24.502 38.811 16.129]
ARG 167 CD altloc B occ 0.5 [24.69  38.671 16.294]

Warning: iterating a Residue yields each DisorderedAtom already resolved to a single selected altloc (highest occupancy, or the first parsed on ties), so a naive get_coord() loop silently discards alternate conformers. Call atom.disordered_get_list() to visit every altloc, or atom.disordered_select('B') to switch the active one before reading coordinates. Also remember that get_coord() returns a numpy float32 array; cast to float64 before accumulating distances or running a least-squares fit to avoid precision loss.

Try it yourself: Load PDB 1AKE with MMCIFParser and flatten chain A to atoms via Selection.unfold_entities(structure[0]['A'], 'A'). Keep only polymer atoms where atom.get_parent().get_id()[0] == ' ', cast each get_coord() to float64, and compute their centroid with numpy.mean. Then loop again, expanding every disordered atom through disordered_get_list(), and report how many extra alternate-conformer coordinates that adds.

PyMOL: Representations, Coloring, and Publication Figures Measuring Distances, Angles, and Dihedrals