stage_structures.py

#!/usr/bin/env python3
"""Download and pre-orient every structure used in the tRNA life-cycle figure.

Panels a-d must show the tRNA in ONE orientation, so each complex is superposed
onto the free tRNA (6UGG chain A) by least squares on the P atoms of the four
conserved helical stems.  Panel e (the three-site ribosome) gets its own frame,
built from the plane through its three tRNA centroids, plus a front cut that
opens the intersubunit space.

Writes staged/*.cif; the ChimeraX scripts in scripts/ open those files.

    python stage_structures.py
"""
import os
import urllib.request

import gemmi
import numpy as np

ENTRIES = {"6ugg": "6UGG", "1c0a": "1C0A", "1b23": "1B23", "1vy5": "1VY5"}

# hero tRNA chain, and the offset between its residue numbering and 1..76
TRNA_CHAIN = {"6ugg": "A", "1c0a": "B", "1b23": "R"}
OFFSET = {"6ugg": 0, "1c0a": 600, "1b23": 0}

# conserved helical stems, in canonical tRNA numbering
CORE = (list(range(1, 8)) + list(range(66, 73))        # acceptor stem
        + list(range(10, 14)) + list(range(22, 26))    # D stem
        + list(range(27, 32)) + list(range(39, 44))    # anticodon stem
        + list(range(49, 54)) + list(range(61, 66)))   # T stem

# landmarks defining the tRNA's own frame: elbow, acceptor end, anticodon loop
LM_ELBOW = [8, 9, 10, 11, 12, 25, 26, 54, 55, 56, 57, 58]
LM_ACCEPTOR = [1, 2, 3, 4, 71, 72, 73, 74, 75, 76]
LM_ANTICODON = [33, 34, 35, 36, 37]

Z_CUT = 8.0          # panel e: drop ribosome residues in front of this plane


def fetch(pdb_id, path):
    if not os.path.exists(path):
        urllib.request.urlretrieve(f"https://files.rcsb.org/download/{pdb_id}.cif", path)
    return path


def read_clean(path):
    s = gemmi.read_structure(path)
    s.setup_entities()
    s.remove_alternative_conformations()
    s.remove_hydrogens()
    s.remove_waters()
    return s


def atom_xyz(s, chain, seqid, name="P"):
    for ch in s[0]:
        if ch.name != chain:
            continue
        for r in ch:
            if r.seqid.num == seqid:
                a = r.find_atom(name, "*")
                if a:
                    return np.array([a.pos.x, a.pos.y, a.pos.z])
    return None


def centroid(s, chain, nums, name="P"):
    pts = [atom_xyz(s, chain, n, name) for n in nums]
    pts = [p for p in pts if p is not None]
    return np.mean(pts, axis=0) if pts else None


def chain_xyz(s, names):
    return np.array([[a.pos.x, a.pos.y, a.pos.z]
                     for ch in s[0] if ch.name in names for r in ch for a in r])


def apply_tr(s, R, t):
    for model in s:
        for ch in model:
            for r in ch:
                for a in r:
                    q = R @ np.array([a.pos.x, a.pos.y, a.pos.z]) + t
                    a.pos = gemmi.Position(*q)


def kabsch(P, Q):
    """Rotation+translation minimising |R@P + t - Q|."""
    Pc, Qc = P.mean(0), Q.mean(0)
    H = (P - Pc).T @ (Q - Qc)
    U, _, Vt = np.linalg.svd(H)
    D = np.diag([1, 1, np.sign(np.linalg.det(Vt.T @ U.T))])
    R = Vt.T @ D @ U.T
    return R, Qc - R @ Pc


def trna_frame(s, chain, off=0):
    """Orthonormal frame from elbow -> acceptor (u) and elbow -> anticodon (v)."""
    E = centroid(s, chain, [n + off for n in LM_ELBOW])
    A = centroid(s, chain, [n + off for n in LM_ACCEPTOR])
    B = centroid(s, chain, [n + off for n in LM_ANTICODON])
    u = A - E
    u /= np.linalg.norm(u)
    v = B - E
    v -= u * np.dot(v, u)
    v /= np.linalg.norm(v)
    return E, np.vstack([u, v, np.cross(u, v)])


def target_frame():
    """Textbook L: acceptor arm up-right, anticodon arm down-left, in the xy plane."""
    tu = np.array([0.94, 0.30, 0.0])
    tu /= np.linalg.norm(tu)
    tv = np.array([0.30, -0.94, 0.0])
    tv -= tu * np.dot(tv, tu)
    tv /= np.linalg.norm(tv)
    return np.vstack([tu, tv, np.cross(tu, tv)])


def stage_panels_ad(outdir="staged", structdir="structures"):
    """Orient 6UGG, then superpose 1C0A and 1B23 onto it. Returns fit RMSDs."""
    s = {pid: read_clean(fetch(ENTRIES[pid], f"{structdir}/{pid}.cif"))
         for pid in ("6ugg", "1c0a", "1b23")}

    T = target_frame()
    E, M = trna_frame(s["6ugg"], "A", 0)
    R0 = T.T @ M
    apply_tr(s["6ugg"], R0, -R0 @ E)

    rmsds = {}
    for pid in ("1c0a", "1b23"):
        P, Q = [], []
        for n in CORE:
            a = atom_xyz(s["6ugg"], "A", n)
            b = atom_xyz(s[pid], TRNA_CHAIN[pid], n + OFFSET[pid])
            if a is not None and b is not None:
                Q.append(a)
                P.append(b)
        P, Q = np.array(P), np.array(Q)
        R, t = kabsch(P, Q)
        apply_tr(s[pid], R, t)
        d = np.linalg.norm((R @ P.T).T + t - Q, axis=1)
        rmsds[pid] = (len(d), float(np.sqrt(np.mean(d ** 2))))

    for pid, st in s.items():
        st.setup_entities()
        st.make_mmcif_document().write_file(f"{outdir}/{pid}_aln.cif")
    return rmsds


def stage_panel_e(outdir="staged", structdir="structures"):
    """1VY5: keep one 70S copy, orient on the three-tRNA plane, cut the front off."""
    s = read_clean(fetch(ENTRIES["1vy5"], f"{structdir}/1vy5.cif"))

    # two 70S copies per asymmetric unit: A*/B* is the first, C*/D* the second
    for name in [ch.name for ch in s[0] if ch.name[0] not in ("A", "B")]:
        s[0].remove_chain(name)
    s.setup_entities()

    s30 = [ch.name for ch in s[0] if ch.name.startswith("A")
           and ch.name not in ("AV", "AW", "AX", "AY")]
    s50 = [ch.name for ch in s[0] if ch.name.startswith("B")]

    cA = chain_xyz(s, {"AW"}).mean(0)      # A-site tRNA
    cP = chain_xyz(s, {"AX"}).mean(0)      # P-site tRNA
    cE = chain_xyz(s, {"AY"}).mean(0)      # E-site tRNA
    n = np.cross(cP - cE, cA - cP)
    n /= np.linalg.norm(n)
    xax = cA - cE                          # E on the left, A on the right
    xax -= n * np.dot(xax, n)
    xax /= np.linalg.norm(xax)
    yax = np.cross(n, xax)
    mid = (cA + cP + cE) / 3
    if np.dot(chain_xyz(s, set(s30)).mean(0) - mid, yax) > 0:
        yax, n = -yax, -n                  # force the 30S to the bottom of the frame
    R = np.vstack([xax, yax, n])
    apply_tr(s, R, -R @ mid)
    s.setup_entities()
    s.make_mmcif_document().write_file(f"{outdir}/1vy5_epa.cif")

    # cutaway: remove ribosome residues in front of the tRNAs (never the tRNAs/mRNA)
    cut = read_clean(f"{outdir}/1vy5_epa.cif")
    for ch in cut[0]:
        if ch.name in ("AV", "AW", "AX", "AY"):
            continue
        doomed = [r.seqid.num for r in ch
                  if np.mean([a.pos.z for a in r]) > Z_CUT]
        for num in doomed:
            for i, r in enumerate(ch):
                if r.seqid.num == num:
                    del ch[i]
                    break
    cut.setup_entities()
    cut.make_mmcif_document().write_file(f"{outdir}/1vy5_epa_cut.cif")
    return {"30S_chains": len(s30), "50S_chains": len(s50)}


if __name__ == "__main__":
    os.makedirs("structures", exist_ok=True)
    os.makedirs("staged", exist_ok=True)
    fits = stage_panels_ad()
    for pid, (n, rms) in fits.items():
        print(f"{pid}: core-stem fit on {n} P atoms, rmsd {rms:.2f} A")
    print(stage_panel_e())
    # 6UGG side view (panel b) is the same file, rotated by ChimeraX at render time