bio-core-protein-structure

name	bio-core-protein-structure
description	Parse PDB files with BioPython Bio.PDB, compute distance/angle/dihedral/RMSD, superimpose structures (Kabsch), assign DSSP secondary structure. Use for protein structure analysis, RMSD/alignment, contact maps, Ramachandran plots.
tool_type	python
primary_tool	biopython

Protein Structure Analysis

When to Use

• Parsing a downloaded or local PDB/mmCIF file and navigating its Structure→Model→Chain→Residue→Atom hierarchy.
• Computing bond distances, bond angles, dihedral (phi/psi/omega) angles, or RMSD between two conformers/structures.
• Superimposing two structures (Kabsch algorithm or Bio.PDB.Superimposer) to compare conformations.
• Assigning secondary structure (helix/sheet/coil) and solvent accessibility with DSSP.
• Building a Ramachandran plot or a CA-CA contact map from a chain.

Version Compatibility

• biopython >= 1.81 (Bio.PDB, Bio.PDB.DSSP, Bio.PDB.Polypeptide.PPBuilder)
• mkdssp (DSSP) >= 3.0, from the salilab conda channel
• numpy >= 1.24, matplotlib >= 3.7, Python >= 3.9

Prerequisites

• pip install biopython numpy matplotlib
• conda install -c salilab dssp (provides the mkdssp executable; skip if you don't need secondary structure)
• Familiarity with the PDB fixed-column ATOM format and basic numpy array ops.

Goal: Load a structure file and extract per-chain CA coordinates and sequence. Approach: Download/parse with PDBParser, filter out heteroatoms via residue.id[0] == ' ', walk the SMCRA tree.

from Bio.PDB import PDBParser, PDBList
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# Download a structure from RCSB PDB (crambin, 1CRN -- small, well-resolved, 46 residues)
pdbl = PDBList()
pdb_file = pdbl.retrieve_pdb_file('1CRN', pdir='pdb_files', file_format='pdb')

parser = PDBParser(QUIET=True)
structure = parser.get_structure('crambin', pdb_file)

def extract_ca_trace(chain):
    """Return (res_name, res_num, coord) for every standard-residue CA atom in a chain.

    Skips heteroatoms/waters (residue.id[0] != ' ') and residues missing a CA
    (e.g. some N-terminal or disordered residues).
    """
    trace = []
    for residue in chain:
        if residue.id[0] != ' ':
            continue
        if 'CA' in residue:
            ca = residue['CA']
            trace.append((residue.get_resname(), residue.id[1], ca.get_vector().get_array()))
    return trace

model = structure[0]              # index 0 always exists: sole conformer for X-ray, first NMR model
chain = model['A']
ca_trace = extract_ca_trace(chain)
print(f"{len(ca_trace)} CA atoms in chain A")

# Residue access shorthand (only valid when hetflag=' ' and icode=' '):
residue = chain[10]                # same as chain[(' ', 10, ' ')]

Goal: Compute geometric measures (distance, angle, dihedral, RMSD) and a CA-CA contact map. Approach: Plain numpy vector math; BioPython's Atom.__sub__ gives distance directly for two atoms.

import numpy as np

def distance(coord1, coord2):
    """Euclidean distance between two 3D points."""
    return np.linalg.norm(np.array(coord1) - np.array(coord2))

def angle(coord1, coord2, coord3):
    """Angle in degrees at coord2, formed by coord1-coord2-coord3."""
    v1 = np.array(coord1) - np.array(coord2)
    v2 = np.array(coord3) - np.array(coord2)
    cos_angle = np.clip(np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)), -1.0, 1.0)
    return np.degrees(np.arccos(cos_angle))

def dihedral(p1, p2, p3, p4):
    """Dihedral (torsion) angle in degrees defined by four consecutive points."""
    p1, p2, p3, p4 = [np.array(p) for p in (p1, p2, p3, p4)]
    b1, b2, b3 = p2 - p1, p3 - p2, p4 - p3
    n1, n2 = np.cross(b1, b2), np.cross(b2, b3)
    n1, n2 = n1 / np.linalg.norm(n1), n2 / np.linalg.norm(n2)
    m1 = np.cross(n1, b2 / np.linalg.norm(b2))
    return np.degrees(np.arctan2(np.dot(m1, n2), np.dot(n1, n2)))

def calculate_rmsd(coords1, coords2):
    """RMSD (Angstrom) between two already-aligned Nx3 coordinate arrays.

    RMSD (Å): 0-1 nearly identical | 1-2 same fold | 2-3 similar fold, some
    variation | 3-5 same topology, different details | >5 different structures.
    """
    coords1, coords2 = np.array(coords1), np.array(coords2)
    if coords1.shape != coords2.shape:
        raise ValueError(f"Shape mismatch: {coords1.shape} vs {coords2.shape}")
    return np.sqrt(np.mean(np.sum((coords1 - coords2) ** 2, axis=1)))

# BioPython shorthand: atom1 - atom2 returns distance in Angstroms directly
d = chain[][] - chain[][]


ca_list = [res[]  res  chain  res.[] ==     res]
contacts = [(i + , j + , ca_list[i] - ca_list[j])
             i  ((ca_list))
             j  (i + , (ca_list))
             ca_list[i] - ca_list[j] < ]

Goal: Superimpose two conformers and report RMSD after optimal alignment (not just raw RMSD, which is meaningless before alignment). Approach: Kabsch algorithm (SVD-based rotation) for full control; Bio.PDB.Superimposer for a one-line production path. Always align on Cα atoms only and report the number of residues aligned alongside the RMSD.

import numpy as np
from Bio.PDB import Superimposer

def kabsch_superpose(mobile, reference):
    """Superimpose `mobile` (Nx3) onto `reference` (Nx3) via the Kabsch algorithm.

    Returns (rotated_coords, rmsd, rotation_matrix). Both inputs must already
    be in 1:1 correspondence (same residue order/count).
    """
    mobile, reference = np.array(mobile, dtype=float), np.array(reference, dtype=float)
    mobile_c = mobile - mobile.mean(axis=0)
    ref_c = reference - reference.mean(axis=0)

    H = mobile_c.T @ ref_c
    U, S, Vt = np.linalg.svd(H)
    d = np.sign(np.linalg.det(Vt.T @ U.T))          # correct for reflection
    R = Vt.T @ np.diag([1.0, 1.0, d]) @ U.T

    rotated = (mobile_c @ R) + reference.mean(axis=0)
    rmsd = np.sqrt(np.mean(np.sum((rotated - reference) ** 2, axis=1)))
    return rotated, rmsd, R

# Production shortcut: BioPython's built-in Superimposer (needs matching Atom objects)
sup = Superimposer()
sup.set_atoms(fixed_atoms, moving_atoms)             # lists of Bio.PDB.Atom, same length/order
sup.apply(moving_atoms)                              # rotates moving_atoms in place onto fixed_atoms
print(f"RMSD: {sup.rms:.3f} A over {len(fixed_atoms)} atoms")

Goal: Assign secondary structure (helix/sheet/coil) and per-residue solvent accessibility. Approach: Run DSSP via BioPython's wrapper (requires the external mkdssp binary); fall back gracefully if it isn't installed.

from Bio.PDB.DSSP import DSSP

def summarize_secondary_structure(model, pdb_file):
    """Run DSSP on `model` and return (ss_sequence, helix_frac, sheet_frac, coil_frac).

    DSSP codes: H=alpha-helix, G=3-10 helix, I=pi-helix (helix group);
    E=beta-strand, B=isolated bridge (sheet group); T/S/'-' = turn/bend/coil.
    Raises if `mkdssp` is not on PATH -- install with `conda install -c salilab dssp`.
    """
    dssp = DSSP(model, pdb_file, dssp='mkdssp')
    ss_sequence = ''.join(dssp[key][2] for key in dssp.keys())
    total = len(ss_sequence)
    helix = sum(ss_sequence.count(c) for c in 'HGI') / total
    sheet = sum(ss_sequence.count(c) for c in 'EB') / total
    return ss_sequence, helix, sheet, 1 - helix - sheet

try:
    ss_sequence, helix_frac, sheet_frac, coil_frac = summarize_secondary_structure(structure[0], pdb_file)
    print(f"Helix {helix_frac:.1%}  Sheet {sheet_frac:.1%}  Coil {coil_frac:.1%}")
except Exception as e:
    print(f"DSSP unavailable ({e}); install with: conda install -c salilab dssp")

Pitfalls

• PDB uses fixed-width columns, not whitespace: ATOM records are column-delimited (name cols 13-16, chain col 22, resseq 23-26, x/y/z 31-54). Never .split() an ATOM line by hand — use column slicing or BioPython.
• SMCRA hierarchy: Structure → Model → Chain → Residue → Atom. X-ray = one model (index 0); NMR = multiple models (each a conformer). Residue IDs are tuples (hetflag, resseq, icode): (' ', 10, ' ') for standard residues, ('H_NAG', 1, ' ') for heteroatoms.
• Resolution: lower is better. <2.0 Å resolves atoms and waters; 3.0-4.0 Å only backbone is reliable, side chains are approximate.
• B-factor vs pLDDT: in experimental structures, B-factor >60 Ų means flexible/disordered. In AlphaFold2 output, the same column stores pLDDT (0-100); pLDDT >90 is very confident, 50-70 is low confidence — don't compare the two scales directly.
• Raw RMSD is meaningless without alignment: two identical structures translated apart have huge RMSD. Always superimpose first (Kabsch/Superimposer), align on Cα only, and report the residue count aligned alongside the value.
• DSSP needs the external mkdssp binary, not just the BioPython wrapper — pip install biopython alone will raise FileNotFoundError at DSSP(...) call time.

See Also

• bio-applied-structural-methods — downstream structural methods (docking, SASA, structure-based ML).
• bio-core-nucleic-acid-structure — same SMCRA/geometry patterns applied to DNA/RNA structures.
• alphafold-structure-prediction — generating a predicted structure when no PDB entry exists.
• bio-core-biopython-essentials — general BioPython patterns beyond Bio.PDB.