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.

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.

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

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.

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

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

Secondary Structure Geometry

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