# Bio Core Protein Structure

> 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.

- Skill: `pavel-kravchenko/bio-core-protein-structure` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-protein-structure`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-protein-structure/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-protein-structure

---


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

```python
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.

```python
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[1]['CA'] - chain[2]['CA']

# CA-CA contact map: pairs within 8 Å, excluding near-sequence neighbors (|i-j| < 4)
ca_list = [res['CA'] for res in chain if res.id[0] == ' ' and 'CA' in res]
contacts = [(i + 1, j + 1, ca_list[i] - ca_list[j])
            for i in range(len(ca_list))
            for j in range(i + 4, len(ca_list))
            if ca_list[i] - ca_list[j] < 8.0]
```

**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.

```python
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.

```python
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`.