bio-structural-biology-structure-navigation

name	bio-structural-biology-structure-navigation
description	Navigate the Bio.PDB SMCRA hierarchy (Structure-Model-Chain-Residue-Atom) safely, surfacing the heterogeneity it hides by default. Use when deciding how to handle altloc/DisorderedAtom conformers before a distance or RMSD, indexing residues insertion-code-safe with the full (hetflag, resseq, icode) tuple, choosing the ATOM/observed vs SEQRES/canonical vs UniProt sequence, selecting the right Model for an NMR ensemble, filtering waters/hetero/metals correctly, and reconciling auth vs label numbering. Keywords SMCRA, altloc, DisorderedAtom, insertion code, SEQRES, PPBuilder, auth_seq_id.
tool_type	python
primary_tool	Bio.PDB
goal_approach_exempt	true

Version Compatibility

Reference examples tested with: biopython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Structure Navigation

"Walk the chains, residues, and atoms; pull out the sequence" -> Traverse the Structure-Model-Chain-Residue-Atom (SMCRA) tree, but treat every level as a lossy projection that hides heterogeneity unless asked otherwise.

• Python: structure[0]['A'][(' ', 100, ' ')]['CA'].coord for full-tuple direct access

Governing Principle: SMCRA is a convenient tree that HIDES heterogeneity by default

The SMCRA model (Structure > Model > Chain > Residue > Atom) is a readable in-memory tree, but its defaults quietly collapse the very heterogeneity that changes the answer. Five traps recur, and none of them raises an error - the code runs and returns a plausible-looking number computed on the wrong thing.

1. A DisorderedAtom silently forwards every uncaught call to ONE child - the HIGHEST-OCCUPANCY altloc, not literally altloc 'A'. So get_atoms(), .coord, distances, clashes, and RMSDs all use a single conformer, invisibly, even when the active site is 60/40 disordered. This is the #1 correctness trap. Enumerate with is_disordered() and disordered_get_list(); never let both altlocs of one atom enter the same geometric calculation.
2. The residue id is a 3-tuple (hetflag, resseq, icode). Naive residue.id[1] drops both the hetero flag and the INSERTION CODE, so antibody residues 100, 100A, 100B collapse onto one key; chain[100] works only until an insertion code, a hetero residue, or an altloc residue exists at that number, then raises a KeyError that looks like the residue is missing. Key on the full tuple; reduce to id[1] only for display.
3. The sequence PPBuilder extracts is the OBSERVED (ATOM-record) sequence with missing-density gaps silently concatenated away - it is NOT the SEQRES/construct/UniProt canonical sequence. A disordered 12-residue loop becomes 12 vanished characters with no marker, so mapping conservation or alignment columns by string position is off-by-many after the first gap. Map through residue NUMBERS (auth_seq_id) or SIFTS, never by string index.
4. NMR and multi-state files have multiple Model objects. Iterating chains without first selecting a model conflates conformers; structure[0] silently picks one NMR member and over-claims precision. Ask "how many models and why" first.
5. mmCIF carries two numbering schemes: auth (matches the paper, has insertion codes, can be negative/gapped) and label (gapless 1..N, no icodes). MMCIFParser defaults to auth; flip auth_residues=False and residue 100 becomes a different residue. Pick one scheme and stay in it.

Decision: which sequence source

These three "sequences" are routinely conflated; each answers a different question and the wrong one silently misindexes everything downstream.

Source	What it is	Get it via	Use when	Fails when
ATOM / observed	Only residues with modeled coordinates; gaps concatenated away	PPBuilder().build_peptides() then pp.get_sequence()	Per-atom geometry, contacts, extracting exactly what was resolved	Aligning to UniProt/MSA by string position (gaps shift the frame)
SEQRES / declared	Full sequence the depositor says is in the crystal, including unresolved residues	SeqIO.parse(file, 'pdb-seqres') or 'cif-seqres'	Knowing the true construct length, locating missing loops	Assuming every SEQRES residue has coordinates (it does not)
UniProt / canonical	The reference biological sequence (no tags, no engineered mutations)	database-access/uniprot-access + SIFTS residue mapping	Mapping conservation/domains/mutations onto structure positions	Assuming construct == canonical (tags, point mutations, chimeras differ)

Decision: residue selection idiom

Filter on the hetflag (id[0]), not the residue name, and know exactly what each idiom keeps and drops.

Goal	Idiom	Keeps / drops correctly?
Standard amino acids only	r.id[0] == ' '	Correct; drops water, ligands, and modified residues
Water	r.id[0] == 'W'	Correct; water hetflag is 'W', NOT 'H_' - a startswith('H_') water strip MISSES water
Ligands and hetero groups	r.id[0].startswith('H_')	Also catches modified residues (MSE, SEP, PTR) mid-chain - not just free ligands
Strip hetero blindly	r.id[0] != ' '	DANGEROUS - deletes catalytic metals, cofactors, AND selenomethionine (MSE) out of the chain
A specific residue type	r.resname == 'ARG'	Fine for type queries; never use resname to classify water vs ligand vs standard

Required Imports

from Bio.PDB import PDBParser, MMCIFParser, PPBuilder, CaPPBuilder, Selection
from Bio.Data.PDBData import protein_letters_3to1, protein_letters_3to1_extended

Accessing Hierarchy Levels

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

model = structure[0]                       # first model (see NMR caveat below)
chain = model['A']
residue = chain[(' ', 100, ' ')]           # full id tuple - insertion-code and hetero safe
residue_bare = chain[100]                   # convenience path; breaks on icode/hetero/altloc at 100
atom = residue['CA']

Iterating Over Structure

for model in structure:
    for chain in model:
        for residue in chain:
            hetflag, resseq, icode = residue.id   # keep the whole tuple, not resseq alone
            for atom in residue:
                print(f'{chain.id}:{resseq}{icode}:{atom.name}')

for chain in structure.get_chains():
    print(f'Chain: {chain.id}')

Enumerating Disordered Atoms (the highest-value pattern)

# A DisorderedAtom forwards uncaught calls to its highest-occupancy child by default,
# so plain iteration measures ONE conformer. Enumerate every altloc explicitly.
for atom in residue:
    if atom.is_disordered():
        for alt in atom.disordered_get_list():   # each alt is a real Atom with its own coord/occupancy
            print(f'{atom.name} altloc {alt.altloc}: occ={alt.occupancy} coord={alt.coord}')
    else:
        print(f'{atom.name}: coord={atom.coord}')

# disordered_select(altloc) mutates the active child GLOBALLY - reset it when done
if atom.is_disordered():
    atom.disordered_select(atom.disordered_get_id_list()[0])

Disordered Residues (microheterogeneity / point mutation at one site)

# Residue.is_disordered() returns 1 for a NORMAL residue that merely holds altloc atoms and 2 for a
# true DisorderedResidue (two resnames at one position); disordered_get_id_list/disordered_select
# exist ONLY on the latter, so gate on == 2 (or isinstance DisorderedResidue) or this AttributeErrors.
if residue.is_disordered() == 2:
    names = residue.disordered_get_id_list()      # alternative resnames at this position
    residue.disordered_select('ALA')              # pick one by resname before any geometry

Extracting the Observed Sequence (know it has silent gaps)

# PPBuilder builds peptides from OBSERVED atoms via a C-N distance criterion and breaks at gaps;
# the returned Seq has missing-density loops concatenated away with no gap marker.
ppb = PPBuilder()
for pp in ppb.build_peptides(structure):
    print(f'observed segment len={len(pp.get_sequence())}: {pp.get_sequence()}')

# CaPPBuilder connects residues whose CA atoms are within ~4.3 Angstroms, so it bridges
# small backbone breaks - useful for CA-only/broken chains, but it can mis-join true gaps.
ca_ppb = CaPPBuilder()
segments = ca_ppb.build_peptides(structure)

Reading the Declared (SEQRES) Sequence and Locating Gaps

from Bio import SeqIO

# SEQRES = the full declared sequence, including residues with no coordinates.
for record in SeqIO.parse('protein.pdb', 'pdb-seqres'):
    print(f'{record.id} declared length {len(record.seq)}')

# header['missing_residues'] (populated when get_header=True) lists unmodeled residues -
# the difference between SEQRES and observed. Reconcile before mapping to UniProt.
parser = PDBParser(QUIET=True, get_header=True)
structure = parser.get_structure('protein', 'protein.pdb')
missing = structure.header.get('missing_residues', [])

Converting Residue Names (modified residues map to X or drop)

# protein_letters_3to1 is the strict 20-aa map (unknown resname -> KeyError, so use .get).
# protein_letters_3to1_extended additionally maps modified residues (MSE->M, SEP->S, PTR->Y).
seq = ''
for residue in chain:
    if residue.id[0] == ' ' or residue.resname in protein_letters_3to1_extended:
        seq += protein_letters_3to1_extended.get(residue.resname, 'X')

Selecting Entities and Full Identifiers

residues = Selection.unfold_entities(structure, 'R')   # S/M/C/R/A level codes
atoms = Selection.unfold_entities(chain, 'A')

atom = structure[0]['A'][(' ', 100, ' ')]['CA']
print(atom.get_full_id())   # ('protein', 0, 'A', (' ', 100, ' '), ('CA', ' '))

Working with NMR Ensembles (do not conflate models)

n_models = len(structure)                  # NMR = conformer ensemble, X-ray/cryo-EM usually 1
if n_models > 1:
    # Compute per-model and report the distribution; never average coordinates across models.
    for model in structure:
        ca = [r['CA'].coord for r in model.get_residues() if r.has_id('CA')]
        print(f'model {model.id}: {len(ca)} CA atoms')

Reading mmCIF with an Explicit Numbering Scheme

# MMCIFParser defaults to auth numbering (matches the paper, has insertion codes).
# label numbering is gapless 1..N with no icodes - a different residue at the same number.
cif_parser = MMCIFParser(QUIET=True, auth_residues=True)   # keep auth for literature/UniProt cross-ref
structure = cif_parser.get_structure('protein', 'protein.cif')

Common Errors

Symptom	Cause	Fix
Distance/RMSD subtly off on a disordered site	get_atoms() returned only the highest-occupancy altloc of a DisorderedAtom	Enumerate disordered_get_list(); pick one altloc consistently before geometry
Impossibly close contacts / inflated atom count	Both altlocs of one atom entered the same calculation	Select a single altloc per site; never mix conformers
KeyError on chain[100] for a residue that is clearly present	Residue has an insertion code, hetero flag, or altloc, so the bare-int path misses it	Index with the full tuple chain[(' ', 100, ' ')] or iterate and match id[1]/id[2]
Antibody CDR residues 100/100A/100B collapse to one	Keyed on residue.id[1] (resseq) and dropped id[2] (icode)	Key on the full (hetflag, resseq, icode) tuple
Structure sequence one residue shorter than expected after each loop	PPBuilder returns the observed sequence with missing-density gaps concatenated away	Use SEQRES (pdb-seqres) for length; map to UniProt by residue number or SIFTS, not string index
Selenomethionine protein reads full of X or has holes in the chain	MSE is a hetero (H_MSE) residue; strict protein_letters_3to1 returns X or a blanket hetero strip deleted it	Use protein_letters_3to1_extended (MSE->M); filter hetero by explicit allow/deny list
Water strip leaves waters behind	Filtered with startswith('H_'); water hetflag is 'W'	Strip water with r.id[0] == 'W'
Stripping hetero removed a catalytic metal or cofactor	Blanket r.id[0] != ' ' deletes metals, heme, FAD, ions, and mid-chain MSE	Remove only water/buffer by explicit list; keep ligands and metals
NMR metric multiplied ~20x or averaged to nonsense	Iterated all models, or averaged coordinates across the ensemble	Select one representative model, or compute per-model and report the spread
mmCIF residue numbers do not match the paper

Related Skills

• structure-io - Parse and write PDB/mmCIF; auth vs label numbering at the I/O layer
• geometric-analysis - Distances, angles, RMSD, SASA once heterogeneity is resolved
• structure-modification - Strip waters/hetero, edit coordinates and B-factors safely
• interface-analysis - Requires the biological assembly, not the deposited asymmetric unit
• sequence-manipulation/seq-objects - Work with the extracted Seq objects
• alignment/msa-parsing - Map SEQRES/ATOM sequences onto alignment columns
• database-access/uniprot-access - Fetch the canonical sequence for SIFTS-based mapping

References

• Hamelryck T, Manderick B (2003). PDB file parser and structure class implemented in Python. Bioinformatics 19(17):2308-2310. DOI 10.1093/bioinformatics/btg332.
• Cock PJA, Antao T, Chang JT, et al. (2009). Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics 25(11):1422-1423. DOI 10.1093/bioinformatics/btp163.
• Berman HM, Westbrook J, Feng Z, et al. (2000). The Protein Data Bank. Nucleic Acids Research 28(1):235-242. DOI 10.1093/nar/28.1.235.
• Velankar S, Dana JM, Jacobsen J, et al. (2013). SIFTS: Structure Integration with Function, Taxonomy and Sequences resource. Nucleic Acids Research 41(D1):D483-D489. DOI 10.1093/nar/gks1258.