Bio Structural Biology Structure Navigation

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

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']

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}')

# 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])

# 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

# 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)

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', [])

# 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')

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', ' '))

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')

# 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')