Bio.PDB

In one line: parse PDB and mmCIF structure files into a Structure → Model → Chain → Residue → Atom hierarchy, and do geometry on it.

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)

from Bio.PDB import PDBParser, MMCIFParser

parser = PDBParser(QUIET=True)                # QUIET suppresses warning spam
structure = parser.get_structure("1TUP", "1tup.pdb")

parser = MMCIFParser(QUIET=True)               # mmCIF is the modern format
structure = parser.get_structure("1tup", "1tup.cif")

model = structure[0]
chain = model["A"]
residue = chain[175]                            # by residue number
atom = residue["CA"]
atom.coord                                       # numpy array [x, y, z]
atom.bfactor; atom.occupancy; atom.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

for model in structure:
    for chain in model:
        for residue in chain:
            for atom in residue:
                ...

structure.get_atoms(); structure.get_residues(); structure.get_chains()
list(structure.get_atoms())

from Bio.PDB import Selection
Selection.unfold_entities(structure, "R")       # all residues
Selection.unfold_entities(chain, "A")            # all atoms in a chain

Residue identifiers

residue.id            # (hetflag, resseq, icode)  e.g. (' ', 175, ' ')
residue.get_resname() # 'ARG'
residue.id[0]         # ' ' = standard amino acid; 'W' = water; 'H_XXX' = hetero

Filter out waters and heteroatoms before geometry:

aa = [r for r in chain if r.id[0] == " "]

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

from Bio.PDB import PDBIO, MMCIFIO, Select

io = PDBIO()
io.set_structure(structure)
io.save("out.pdb")

class ChainSelect(Select):
    def accept_chain(self, chain):
        return chain.id == "A"
    def accept_residue(self, residue):
        return residue.id[0] == " "          # standard residues only

io.save("chainA.pdb", ChainSelect())

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

from Bio.PDB import PDBList

pdbl = PDBList()
pdbl.retrieve_pdb_file("1TUP", pdir="structures", file_format="mmCif")
pdbl.download_pdb_files(["1TUP", "2AC0"], pdir="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+.

Bio.PDB

Bio.PDB

The SMCRA hierarchy

Iterating

Residue identifiers

Geometry

Neighbour search

Superposition and RMSD

Writing

Downloading structures

Bioinformatics examples

Common mistakes

See also