# Structural Bioinformatics

> Parse PDB structures with Bio.PDB; compute RMSD/TM-score via Kabsch superposition; run DSSP/Ramachandran and PWM/PROSITE scans; GO/KEGG enrichment. Use when parsing PDB files, computing RMSD, or enriching genes via GO/KEGG.

- Skill: `pavel-kravchenko/structural-bioinformatics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/structural-bioinformatics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/structural-bioinformatics/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/structural-bioinformatics

---


# Structural Bioinformatics

## When to Use
- Parsing PDB/mmCIF files, navigating the Structure→Model→Chain→Residue→Atom hierarchy
- Computing distances, bond angles, dihedrals, RMSD, or TM-score between structures
- Secondary structure assignment (DSSP) or building a Ramachandran plot
- Building/scanning a PWM (e.g. transcription-factor motif) or a PROSITE regex pattern
- Running GO or KEGG enrichment on a gene list from an experiment

## Version Compatibility
Biopython ≥1.81 (Bio.PDB, Bio.PDB.DSSP, Bio.PDB.Polypeptide.PPBuilder), NumPy ≥1.24, SciPy ≥1.11, Python ≥3.10. External `mkdssp` (DSSP 4.x) for secondary-structure assignment.

## Prerequisites
- `pip install biopython numpy scipy matplotlib`
- `conda install -c salilab dssp` (provides the `mkdssp` binary; required only for the DSSP step)
- Familiarity with PDB fixed-width columns and the SMCRA object model helps but isn't required

**Goal:** Load a PDB structure and pull out CA coordinates for downstream geometry.
**Approach:** `Bio.PDB.PDBParser` builds the SMCRA hierarchy; iterate residues, skip HETATM records, keep CA atoms.

```python
from Bio.PDB import PDBParser, PDBList

def load_ca_trace(pdb_id: str, pdir: str = "pdb_files"):
    """Download a PDB entry and return its CA atoms (chain A) as a list of dicts.

    Residue.id is a tuple (hetflag, resseq, icode); standard amino acids have
    hetflag == ' '. Ligands/water have a non-space hetflag and are skipped.
    """
    pdbl = PDBList()
    pdb_file = pdbl.retrieve_pdb_file(pdb_id, pdir=pdir, file_format="pdb")
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure(pdb_id, pdb_file)

    ca_atoms = []
    for residue in structure[0]["A"]:
        if residue.id[0] != " ":  # skip HETATM (ligands/water)
            continue
        if "CA" in residue:
            ca_atoms.append({
                "res_name": residue.get_resname(),
                "res_num": residue.id[1],
                "coord": residue["CA"].get_vector().get_array(),
            })
    return ca_atoms

# Example: crambin, a small well-resolved 46-residue protein
# ca = load_ca_trace("1CRN")
```

**Goal:** Compare two conformations quantitatively (RMSD, TM-score) after proper superposition.
**Approach:** Superposition (Kabsch or `Bio.PDB.Superimposer`) must precede RMSD — raw RMSD on unaligned coordinates is meaningless. TM-score is length-normalized and comparable across protein sizes.

```python
import numpy as np
from Bio.PDB import Superimposer

def kabsch(mobile: np.ndarray, ref: np.ndarray) -> np.ndarray:
    """Superpose `mobile` (Nx3) onto `ref` (Nx3) and return the rotated/translated coords."""
    mc, rc = mobile - mobile.mean(0), ref - ref.mean(0)
    U, S, Vt = np.linalg.svd(mc.T @ rc)
    d = np.sign(np.linalg.det(Vt.T @ U.T))
    R = Vt.T @ np.diag([1, 1, d]) @ U.T
    return (mc @ R) + ref.mean(0)

def rmsd(c1: np.ndarray, c2: np.ndarray) -> float:
    """Root-mean-square deviation between two equal-length coordinate arrays."""
    diff = np.array(c1) - np.array(c2)
    return np.sqrt(np.mean(np.sum(diff ** 2, axis=1)))

def tm_score(c1: np.ndarray, c2: np.ndarray) -> float:
    """TM-score (length-normalized fold similarity). >0.5 same fold; <0.3 different fold."""
    L = len(c1)
    d0 = max(1.24 * (L - 15) ** (1 / 3) - 1.8, 0.5)
    d = np.sqrt(np.sum((np.array(c1) - np.array(c2)) ** 2, axis=1))
    return np.sum(1 / (1 + (d / d0) ** 2)) / L

# Or let BioPython do the alignment + RMSD in one step:
# sup = Superimposer(); sup.set_atoms(fixed_atoms, moving_atoms)
# sup.apply(moving_structure.get_atoms()); print(sup.rms)
```

**Goal:** Assign secondary structure and phi/psi angles, then draw a Ramachandran plot.
**Approach:** Run DSSP for SS + solvent accessibility; use `PPBuilder` for phi/psi (works even without `mkdssp` installed).

```python
from Bio.PDB.DSSP import DSSP
from Bio.PDB.Polypeptide import PPBuilder
import numpy as np

def dssp_summary(structure, pdb_file: str):
    """Return (ss_string, helix_frac, sheet_frac) from DSSP. Requires mkdssp on PATH."""
    dssp = DSSP(structure[0], pdb_file, dssp="mkdssp")
    ss = "".join(dssp[key][2] for key in dssp.keys())  # H/E/G/I/B/T/S/-
    helix = sum(ss.count(c) for c in "HGI")
    sheet = sum(ss.count(c) for c in "EB")
    return ss, helix / len(ss), sheet / len(ss)

def phi_psi_angles(chain):
    """Compute (phi, psi) in degrees for every residue in `chain` via PPBuilder."""
    ppb = PPBuilder()
    angles = []
    for pp in ppb.build_peptides(chain):
        for i, (phi, psi) in enumerate(pp.get_phi_psi_list()):
            if phi is not None and psi is not None:
                angles.append((pp[i].get_resname(), np.degrees(phi), np.degrees(psi)))
    return angles

# phi_psi = phi_psi_angles(structure[0]['A'])
# plt.scatter([p for _, p, _ in phi_psi], [s for _, _, s in phi_psi]); plt.xlim(-180,180); plt.ylim(-180,180)
```

**Goal:** Build/scan a PWM and translate a PROSITE pattern into a regex for motif search.
**Approach:** Log-odds PWM with pseudocounts; PROSITE `x`→`.`, `{P}`→`[^P]`, `(n)`/`(n,m)`→ regex repeats.

```python
import re
import numpy as np

BASES = ["A", "C", "G", "T"]

def build_pwm(seqs: list[str], pseudocount: float = 0.1) -> np.ndarray:
    """Build a log2-odds PWM from a list of equal-length aligned sequences."""
    pfm = np.zeros((4, len(seqs[0])))
    for seq in seqs:
        for i, b in enumerate(seq.upper()):
            if b in BASES:
                pfm[BASES.index(b), i] += 1
    ppm = (pfm + pseudocount) / (len(seqs) + 4 * pseudocount)
    return np.log2(ppm / 0.25)

def scan_pwm(pwm: np.ndarray, sequence: str, threshold: float = None) -> list[tuple]:
    """Slide `pwm` across `sequence`, return (pos, subseq, score) hits above threshold."""
    L = pwm.shape[1]
    thresh = threshold or 0.6 * np.sum(np.max(pwm, axis=0))
    hits = []
    for i in range(len(sequence) - L + 1):
        s = sum(pwm[BASES.index(b), j]
                for j, b in enumerate(sequence[i:i + L].upper()) if b in BASES)
        if s >= thresh:
            hits.append((i, sequence[i:i + L], s))
    return sorted(hits, key=lambda x: -x[2])

def prosite_to_regex(pattern: str) -> str:
    """Convert a PROSITE pattern (e.g. 'N-{P}-[ST]-{P}') to a Python regex."""
    parts = []
    for elem in pattern.strip(".").split("-"):
        m = re.match(r"^(.+?)\((\d+)(?:,(\d+))?\)$", elem)
        core, low, high = (m.group(1), m.group(2), m.group(3)) if m else (elem, None, None)
        r = ("." if core == "x" else core if core.startswith("[")
             else f"[^{core[1:-1]}]" if core.startswith("{") else core)
        if low:
            r += f"{{{low},{high}}}" if high else f"{{{low}}}"
        parts.append(r)
    return "".join(parts)

# prosite_to_regex("N-{P}-[ST]-{P}") -> 'N[^P][ST][^P]'  (N-glycosylation site)
```

**Goal:** Test whether a gene list is enriched for GO terms or KEGG pathways.
**Approach:** Hypergeometric test per term, then Benjamini-Hochberg FDR (terms are correlated, so Bonferroni is too conservative).

```python
from scipy import stats
import urllib.request

def go_enrichment(gene_list: list[str], term_to_genes: dict[str, set], N: int = 20000) -> list[dict]:
    """Hypergeometric enrichment of gene_list against a term->gene-set mapping.

    N is the background gene universe size (e.g. ~20000 for human protein-coding genes).
    """
    query = set(gene_list)
    n = len(query)
    results = []
    for term, tgenes in term_to_genes.items():
        K, k = len(tgenes), len(query & tgenes)
        if k == 0:
            continue
        p = stats.hypergeom.sf(k - 1, N, K, n)
        results.append({"term": term, "k": k, "K": K, "p": p})
    results.sort(key=lambda r: r["p"])
    for i, r in enumerate(results):
        r["fdr"] = min(r["p"] * len(results) / (i + 1), 1.0)
    return results

def kegg_get(op: str, *args: str) -> str:
    """Call the KEGG REST API, e.g. kegg_get('get', 'hsa04210') or kegg_get('find', 'pathway', 'apoptosis')."""
    url = "https://rest.kegg.jp/" + "/".join([op] + list(args))
    with urllib.request.urlopen(url, timeout=15) as r:
        return r.read().decode()
# Organism codes: hsa=human, mmu=mouse, sce=yeast, eco=E.coli
```

## Quick Reference

### Protein Structure Levels
| Level | Stabilized by |
|-------|--------------|
| Primary | Peptide bonds (N→C sequence) |
| Secondary | Backbone H-bonds (α-helix i→i+4, β-sheet) |
| Tertiary | Hydrophobic core, H-bonds, disulfides (single chain 3D) |
| Quaternary | Same forces, multiple subunits |

### Secondary Structure Geometry
| Element | Phi/Psi | Rise/res |
|---------|---------|----------|
| α-helix | -57/-47° | 1.5 Å |
| 3₁₀-helix | -49/-26° | 2.0 Å |
| β-strand (antiparallel) | -139/+135° | 3.4 Å |

DSSP codes: `H`=α-helix `G`=3₁₀ `I`=π `E`=β-strand `B`=β-bridge `T`=turn `S`=bend `-`=coil

### GO Evidence Hierarchy
- **Experimental** (EXP, IDA, IMP, IGI): highest quality
- **Computational** (ISS, ISO, IBA): medium
- **Automatic** (IEA): lowest — exclude from stringent analyses

## Pitfalls
- **PDB is fixed-width, not whitespace-delimited.** Use `line[30:38]` for X coordinate, not `line.split()`.
- **Residue ID is a tuple, not an int.** `structure[0]['A'][10]` is shorthand for `(' ', 10, ' ')`. Ligands have a non-space hetflag.
- **NMR structures have multiple models.** Use `structure[0]` for the first model; iterate over `structure` for ensemble analysis.
- **DSSP requires the external `mkdssp` binary.** Install via `conda install -c salilab dssp`; `PPBuilder` phi/psi works without it.
- **RMSD without superposition is meaningless.** Always Kabsch-align (or `Superimposer`) first.
- **PWM zero probabilities → -inf log-odds.** Always add a pseudocount before taking the log.
- **GO true-path rule must be applied before enrichment.** Propagate each annotation to all ancestor terms first.
- **Multiple testing in GO/pathway analysis.** Use BH FDR, not Bonferroni — terms are correlated.
- **IEA annotations are auto-assigned and lower quality.** Filter them out for experimental conclusions.
- **PROSITE `{P}` is a negative class, not a quantifier.** It translates to `[^P]` in regex, not `{1}`.

## See Also
- `bio-pathway-analysis-go-enrichment` — dedicated GO enrichment workflows
- `bio-pathway-analysis-kegg-pathways` — KEGG pathway enrichment and mapping
- `bio-structural-biology-structure-io` — broader structure file I/O (mmCIF, multi-format)
- `bio-structural-biology-geometric-analysis` — additional geometric/structural analyses