--- 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 ` 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 ```python from Bio.PDB import PDBParser, MMCIFParser, PPBuilder, CaPPBuilder, Selection from Bio.Data.PDBData import protein_letters_3to1, protein_letters_3to1_extended ``` ## Accessing Hierarchy Levels ```python 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 ```python 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) ```python # 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) ```python # 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) ```python # 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 ```python 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) ```python # 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 ```python 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) ```python 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 ```python # 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 | Read label numbering instead of auth (or mixed a label index into an auth structure) | Set `auth_residues=True` (default) and stay in one scheme | | Missing loop treated as a covalent chain break | Gap in coordinates is unmodeled disorder, not a real break | Reconcile against `header['missing_residues']`/SEQRES; flag gaps as disorder | | Silent atom drops / merged chains on a messy file | `PDBParser(QUIET=True)` suppressed the discontinuous-chain warnings | For unfamiliar files parse without QUIET (or capture warnings) first | ## 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.