From f1887cf9320544321673543ccb7f3167f6431747 Mon Sep 17 00:00:00 2001 From: "Jason Eastman (windows)" Date: Tue, 8 Sep 2026 23:06:29 -0400 Subject: [PATCH] Build a real protein from the Protein Data Bank A protein is one long chain that folds into one particular shape, and the shape is the entire point -- it decides what the protein does, and the same links in a different arrangement do a different job or none. protein.py reads an entry from RCSB and builds it at a size you can walk round. The default draws the fold and not the atoms. Every atom gives a solid lump, because the inside of a protein is packed and the chain you came to see is buried; one point per amino acid, joined up, is the chain itself. --style spacefill gives the lump, which is worth seeing once for exactly that reason. Colours are the secondary structure, and they are read rather than guessed: a PDB file carries HELIX and SHEET records written by whoever solved the structure, so red helix / yellow sheet / white loops is what the crystallo- graphers said, not what some heuristic here decided. It reproduces the known composition of four very different folds -- myoglobin 118 helix and 0 sheet, haemoglobin 448 and 0, both all-alpha globins; GFP 16 and 107, a beta barrel; crambin 21 and 8. The format is punch-card fixed columns and that is the whole difficulty. An atom NAMED "CA" in columns 13-16 is the alpha carbon every amino acid has in the middle of it; an atom whose ELEMENT is "CA" in columns 77-78 is calcium. Split the line on whitespace and the chain trace grows a spur out to wherever the calcium ion is sitting. There is a test for that, and for altLoc, waters, ligands and NMR models, all against inline text and no network. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Kq1Cjkh3bVT1WLzh5FRVQG --- .gitignore | 4 + examples/protein.py | 435 ++++++++++++++++++++++++++++++++++++++++++ tests/test_protein.py | 196 +++++++++++++++++++ 3 files changed, 635 insertions(+) create mode 100644 examples/protein.py create mode 100644 tests/test_protein.py diff --git a/.gitignore b/.gitignore index 2aeb440..b5f92fb 100644 --- a/.gitignore +++ b/.gitignore @@ -47,3 +47,7 @@ examples/data/exo_build_*.json # purpose, so the examples run with no network at all. examples/data/moon_*.npz examples/data/mars_*.npz + +# PDB entries fetched by protein.py. Small, but there are 200,000 of them and +# they never change, so there is no reason to carry any particular one. +examples/data/pdb_*.pdb diff --git a/examples/protein.py b/examples/protein.py new file mode 100644 index 0000000..ab9ba72 --- /dev/null +++ b/examples/protein.py @@ -0,0 +1,435 @@ +#!/usr/bin/env python3 +"""Build a real protein, from the Protein Data Bank, at a size you can walk round. + + python protein.py # ubiquitin, the standard example + python protein.py --pdb 1MBN # myoglobin, the first ever solved + python protein.py --pdb 4HHB --style spacefill + python protein.py --list # some worth looking at + +The coordinates come from the RCSB Protein Data Bank, live. Every atom in the +file was measured -- by X-ray crystallography, mostly -- and the numbers are in +angstroms, which is a tenth of a nanometre. A carbon atom is about 1.7 A +across. You are about to stand next to one the size of a shed. + +WHY A BACKBONE AND NOT A BALL OF ATOMS + +A protein is a single long chain that folds into one particular shape, and the +shape is the whole point: it is what decides what the protein does, and a chain +with the same links in a different arrangement does a different job or no job +at all. Drawn as every atom, that chain is invisible -- you get a solid lump, +because the inside of a protein is packed. + +So the default draws only the path the chain takes, one point per amino acid, +joined up. That is the fold. --style spacefill gives the lump, which is worth +seeing once, because it is what the protein actually is: the fold you were +admiring has no gaps in it anywhere. + +THE COLOURS ARE THE SECONDARY STRUCTURE, AND THEY ARE MEASURED + +Where the chain coils into a spiral it is an ALPHA HELIX; where it runs +alongside itself in flat ribbons it is a BETA SHEET; the rest is loops holding +those together. Those are not guessed here or worked out from the geometry -- +a PDB file records them, in HELIX and SHEET lines put there by the people who +solved the structure, and this reads them. Red is helix, yellow is sheet, white +is everything else. + +Ubiquitin is a good first look precisely because it has both: one long helix +lying across a sheet of five strands. Myoglobin is nothing but helix. +Green fluorescent protein is a barrel of eleven strands with the bit that +glows threaded up the middle. +""" +import argparse +import io +import math +import os +import sys +import time +import urllib.error +import urllib.request + +try: + import numpy as np +except ImportError: + sys.exit("This needs numpy: pip install numpy") + +import connect + +RCSB = "https://files.rcsb.org/download/%s.pdb" +USER_AGENT = "pyncraft protein (https://github.com/jdeast/pyncraft)" + +# Worth a look, with the number of amino acids, because that is what decides +# how long it takes and how big it comes out. +INTERESTING = [ + ("1UBQ", 76, "Ubiquitin -- one helix on a five-strand sheet. The classic."), + ("1CRN", 46, "Crambin. The smallest thing here; good for a quick try."), + ("1MBN", 153, "Myoglobin -- all helix, and the first protein ever solved."), + ("2LYZ", 129, "Lysozyme, from egg white. It cuts open bacteria."), + ("4INS", 102, "Insulin. Four short chains, two molecules in the file."), + ("1EMA", 221, "Green fluorescent protein -- a barrel of eleven strands."), + ("1TIM", 494, "Triosephosphate isomerase -- the TIM barrel, 8 and 8."), + ("4HHB", 574, "Haemoglobin. Four chains, four hemes, carries your oxygen."), + ("6VXX", 2916, "The SARS-CoV-2 spike. Enormous; use --size and be patient."), +] + +# Standard CPK colours, which every chemistry textbook and every molecular +# viewer has used since Corey, Pauling and Koltun built them out of plastic in +# the sixties. Carbon black, oxygen red, nitrogen blue, sulphur yellow. +ELEMENT_BLOCKS = { + "C": "GRAY_CONCRETE", + "N": "BLUE_CONCRETE", + "O": "RED_CONCRETE", + "S": "YELLOW_CONCRETE", + "P": "ORANGE_CONCRETE", + "H": "WHITE_CONCRETE", + "FE": "BROWN_CONCRETE", + "MG": "LIME_CONCRETE", + "ZN": "LIGHT_GRAY_CONCRETE", + "CA": "CYAN_CONCRETE", + "NA": "PURPLE_CONCRETE", + "CL": "GREEN_CONCRETE", +} +DEFAULT_ELEMENT = "PINK_CONCRETE" + +# Van der Waals radii in angstroms -- how big an atom actually is, which is not +# the same as how far its bonds reach. +VDW = {"C": 1.70, "N": 1.55, "O": 1.52, "S": 1.80, "P": 1.80, "H": 1.20, + "FE": 2.05, "MG": 1.73, "ZN": 1.39, "CA": 2.31, "NA": 2.27, "CL": 1.75} +DEFAULT_VDW = 1.70 + +STRUCTURE_BLOCKS = {"H": "RED_CONCRETE", "E": "YELLOW_CONCRETE", "-": "WHITE_CONCRETE"} + +# A rainbow from the start of the chain to the end, which is the other way +# people colour these: it shows you which way the chain runs. +RAINBOW = ["RED_CONCRETE", "ORANGE_CONCRETE", "YELLOW_CONCRETE", "LIME_CONCRETE", + "GREEN_CONCRETE", "CYAN_CONCRETE", "LIGHT_BLUE_CONCRETE", + "BLUE_CONCRETE", "PURPLE_CONCRETE", "MAGENTA_CONCRETE"] + +CHAIN_BLOCKS = ["RED_CONCRETE", "BLUE_CONCRETE", "YELLOW_CONCRETE", + "LIME_CONCRETE", "MAGENTA_CONCRETE", "CYAN_CONCRETE", + "ORANGE_CONCRETE", "WHITE_CONCRETE"] + + +# ── the Protein Data Bank ────────────────────────────────────────────────── + +def fetch(pdb_id, cache_dir, refresh=False): + """The PDB entry, cached on disk. These files are small and never change.""" + pdb_id = pdb_id.strip().upper() + path = os.path.join(cache_dir, "pdb_%s.pdb" % pdb_id) + if os.path.exists(path) and not refresh: + return io.open(path, encoding="latin-1").read() + req = urllib.request.Request(RCSB % pdb_id, headers={"User-Agent": USER_AGENT}) + try: + text = urllib.request.urlopen(req, timeout=120).read().decode("latin-1") + except urllib.error.HTTPError as e: + if e.code == 404: + raise SystemExit("No PDB entry called %r. They are four characters, " + "like 1UBQ. --list shows some." % pdb_id) + raise SystemExit("The PDB said no (%s)." % e) + os.makedirs(cache_dir, exist_ok=True) + io.open(path, "w", encoding="latin-1", newline="").write(text) + return text + + +def parse(text, waters=False, ligands=True): + """Atoms, secondary structure and a title, out of a PDB file. + + The format is punch-card fixed columns, which is why it is parsed by + character position rather than by splitting: an atom called "CA" in columns + 13-16 is an alpha carbon, and one called "CA" in columns 77-78 is a calcium + ion, and splitting on whitespace loses the difference. + + Only the first MODEL is used. An NMR structure holds twenty or so poses of + the same molecule, and drawing them all at once gives a haystack. + """ + atoms, helix, sheet = [], set(), set() + title, model = [], 0 + for ln in text.splitlines(): + rec = ln[:6] + if rec == "TITLE ": + title.append(ln[10:].strip()) + elif rec == "MODEL ": + model += 1 + if model > 1: + break + elif rec == "HELIX ": + # start and end residue of a helix, inclusive, on one chain + try: + ch, a, b = ln[19], int(ln[21:25]), int(ln[33:37]) + except ValueError: + continue + for r in range(a, b + 1): + helix.add((ch, r)) + elif rec == "SHEET ": + try: + ch, a, b = ln[21], int(ln[22:26]), int(ln[33:37]) + except ValueError: + continue + for r in range(a, b + 1): + sheet.add((ch, r)) + elif rec in ("ATOM ", "HETATM"): + resname = ln[17:20].strip() + if rec == "HETATM": + if resname in ("HOH", "DOD", "WAT") and not waters: + continue + if not ligands: + continue + alt = ln[16] + if alt not in (" ", "A"): + continue # one conformer of a disordered side chain + try: + x, y, z = float(ln[30:38]), float(ln[38:46]), float(ln[46:54]) + resseq = int(ln[22:26]) + except ValueError: + continue + element = (ln[76:78].strip() or ln[12:14].strip()).upper() + atoms.append({ + "name": ln[12:16].strip(), "element": element, + "chain": ln[21], "resseq": resseq, "resname": resname, + "hetatm": rec == "HETATM", "xyz": (x, y, z), + }) + return atoms, helix, sheet, " ".join(title) + + +def backbone(atoms): + """The alpha carbons, in order, one per amino acid, split by chain. + + The alpha carbon is the one every amino acid has in the middle of it, so + joining them up in sequence traces the chain itself and nothing else. + """ + chains = {} + for a in atoms: + if a["name"] == "CA" and not a["hetatm"] and a["element"] in ("C", ""): + chains.setdefault(a["chain"], []).append(a) + for ch in chains: + chains[ch].sort(key=lambda a: a["resseq"]) + return chains + + +# ── drawing ──────────────────────────────────────────────────────────────── + +def connected(a, b): + """Points from a to b, each sharing a face with the one before. + + Rounding along a straight line steps diagonally, and two blocks meeting + only at an edge leave a gap you can see through -- which on a chain reads + as the chain being broken, exactly where it is not. + """ + cur = list(a) + pts = [tuple(cur)] + while tuple(cur) != tuple(b): + far = max(range(3), key=lambda i: abs(b[i] - cur[i])) + cur[far] += 1 if b[far] > cur[far] else -1 + pts.append(tuple(cur)) + return pts + + +def ball(radius): + """Offsets of a solid ball, for stamping along a line or onto an atom.""" + r = max(float(radius), 0.5) + n = int(math.ceil(r)) + ax = np.arange(-n, n + 1) + x, y, z = np.meshgrid(ax, ax, ax, indexing="ij") + keep = (x * x + y * y + z * z) <= (r + 0.25) ** 2 + return np.stack([x[keep], y[keep], z[keep]], axis=1) + + +def to_blocks(atoms, size): + """Angstroms to block coordinates, keeping the shape and centring on zero. + + Returns the positions, the scale in blocks per angstrom, and the shape of + the box they fit in. + """ + xyz = np.array([a["xyz"] for a in atoms], dtype=float) + span = xyz.max(axis=0) - xyz.min(axis=0) + scale = float(size) / max(span.max(), 1e-6) + grid = np.rint((xyz - xyz.min(axis=0)) * scale).astype(int) + return grid, scale, grid.max(axis=0) + 1 + + +def build_backbone(atoms, helix, sheet, size, thickness, colour_by): + """The chain, as a tube through every alpha carbon.""" + chains = backbone(atoms) + if not chains: + raise SystemExit("No alpha carbons in this entry, so there is no chain " + "to trace. Try --style spacefill.") + flat = [a for ch in sorted(chains) for a in chains[ch]] + grid, scale, shape = to_blocks(flat, size) + + palette, index = [], {} + + def label(material): + if material not in index: + palette.append(material) + index[material] = len(palette) + return index[material] + + pad = int(math.ceil(thickness)) + 1 + world = np.zeros(tuple(shape + 2 * pad), dtype=np.int32) + offsets = ball(thickness) + + n = len(flat) + pos = 0 + for ch in sorted(chains): + residues = chains[ch] + for i, atom in enumerate(residues): + if colour_by == "structure": + key = (atom["chain"], atom["resseq"]) + mat = STRUCTURE_BLOCKS["H" if key in helix else + "E" if key in sheet else "-"] + elif colour_by == "rainbow": + mat = RAINBOW[int(len(RAINBOW) * pos / max(n, 1)) % len(RAINBOW)] + else: + mat = CHAIN_BLOCKS[sorted(chains).index(ch) % len(CHAIN_BLOCKS)] + lab = label(mat) + + here = tuple(grid[pos] + pad) + stamp = offsets + here + world[stamp[:, 0], stamp[:, 1], stamp[:, 2]] = lab + + # Join it to the previous residue of the same chain. Alpha carbons + # sit 3.8 A apart, which at any useful scale is several blocks, so + # without this the chain is a string of loose beads. + if i: + prev = tuple(grid[pos - 1] + pad) + for p in connected(prev, here)[1:-1]: + stamp = offsets + np.array(p) + world[stamp[:, 0], stamp[:, 1], stamp[:, 2]] = lab + pos += 1 + return world, palette, scale, len(flat) + + +def build_spacefill(atoms, size, colour_by, chains_present): + """Every atom, at the size it really is.""" + grid, scale, shape = to_blocks(atoms, size) + palette, index = [], {} + + def label(material): + if material not in index: + palette.append(material) + index[material] = len(palette) + return index[material] + + radii = {} + for a in atoms: + radii.setdefault(a["element"], VDW.get(a["element"], DEFAULT_VDW)) + pad = int(math.ceil(max(radii.values()) * scale)) + 1 + world = np.zeros(tuple(shape + 2 * pad), dtype=np.int32) + + stamps = {e: ball(r * scale) for e, r in radii.items()} + order = sorted(range(len(atoms)), + key=lambda i: -radii[atoms[i]["element"]]) + for i in order: + a = atoms[i] + if colour_by == "chain": + mat = CHAIN_BLOCKS[chains_present.index(a["chain"]) % len(CHAIN_BLOCKS)] + else: + mat = ELEMENT_BLOCKS.get(a["element"], DEFAULT_ELEMENT) + lab = label(mat) + stamp = stamps[a["element"]] + grid[i] + pad + world[stamp[:, 0], stamp[:, 1], stamp[:, 2]] = lab + return world, palette, scale, len(atoms) + + +def main(): + p = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + connect.add_arguments(p) + p.add_argument("--pdb", default="1UBQ", help="a four character PDB id (default 1UBQ)") + p.add_argument("--style", choices=("backbone", "spacefill"), default="backbone", + help="the fold, or every atom (default backbone)") + p.add_argument("--colour", choices=("structure", "rainbow", "chain", "element"), + default=None, help="default: structure for a backbone, " + "element for spacefill") + p.add_argument("--size", type=int, default=120, + help="longest side, in blocks (default 120)") + p.add_argument("--thickness", type=float, default=1.6, + help="radius of the backbone tube, in blocks (default 1.6)") + p.add_argument("--at", nargs=3, type=int, metavar=("X", "Y", "Z"), + help="where to put it (default: where you are standing)") + p.add_argument("--waters", action="store_true", + help="include the water molecules, which are usually most of " + "the HETATM records and always in the way") + p.add_argument("--no-ligands", action="store_true", + help="protein only: no hemes, no metals, no drugs") + p.add_argument("--pace", type=int, default=12000, help="blocks per second") + p.add_argument("--budget", type=int, default=900000, + help="refuse to place more blocks than this (default 900000)") + p.add_argument("--force", action="store_true", help="build it anyway") + p.add_argument("--refresh", action="store_true", help="re-download the entry") + p.add_argument("--list", action="store_true", help="list some proteins and stop") + p.add_argument("--dry-run", action="store_true", help="print the plan only") + args = p.parse_args() + + if args.list: + print("%-6s %8s %s" % ("pdb", "residues", "what it is")) + for pid, n, what in INTERESTING: + print("%-6s %8d %s" % (pid, n, what)) + print() + print("Any of the ~200,000 entries at rcsb.org works: --pdb .") + return + + data = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data") + text = fetch(args.pdb, data, args.refresh) + atoms, helix, sheet, title = parse(text, waters=args.waters, + ligands=not args.no_ligands) + if not atoms: + raise SystemExit("No atoms in %s, which should not happen." % args.pdb) + + chains_present = sorted({a["chain"] for a in atoms}) + colour = args.colour or ("structure" if args.style == "backbone" else "element") + + print("%s -- %s" % (args.pdb.upper(), title[:70])) + print(" %d atoms, %d chains (%s)" + % (len(atoms), len(chains_present), ", ".join(chains_present))) + print(" %d residues in a helix, %d in a sheet, as recorded by the people " + "who solved it" % (len(helix), len(sheet))) + + if args.style == "backbone": + world, palette, scale, n = build_backbone( + atoms, helix, sheet, args.size, args.thickness, colour) + print(" tracing %d alpha carbons" % n) + else: + world, palette, scale, n = build_spacefill( + atoms, args.size, colour, chains_present) + print(" drawing all %d atoms at their van der Waals radii" % n) + + blocks = int((world != 0).sum()) + nx, ny, nz = world.shape + print() + print(" 1 block = %.3f angstroms, so a carbon atom is %.1f blocks across" + % (1.0 / scale, 2 * 1.7 * scale)) + print(" %d x %d x %d, %d blocks, about %.0f seconds" + % (nx, ny, nz, blocks, blocks / float(args.pace))) + print(" colours: %s" % ", ".join(palette)) + + if blocks > args.budget and not args.force: + fits = args.size * (float(args.budget) / blocks) ** (1.0 / 3.0) + print() + print(" That is over the budget of %d." % args.budget) + print(" --size %d would fit; --force builds it anyway." % int(fits * 0.95)) + return + + if args.dry_run: + return + + mc = connect.connect(args.host, args.port, args.player) + if args.at: + ox, oy, oz = args.at + else: + pos = mc.player.getTilePos() + ox, oy, oz = pos.x, pos.y + 2, pos.z + print() + print("building at %d, %d, %d" % (ox, oy, oz)) + mc.postToChat("%s: %s" % (args.pdb.upper(), title[:80])) + started = time.time() + sent = mc.buildVoxels(world, palette=palette, origin=(ox, oy, oz), + blocks_per_second=args.pace) + print(" %d commands in %.0fs" % (sent, time.time() - started)) + if args.style == "backbone" and colour == "structure": + mc.postToChat("Red is alpha helix, yellow is beta sheet, white is the " + "loops that hold them together.") + print(" stand at %d %d %d" % (ox + nx // 2, oy + ny + 4, oz + nz // 2)) + + +if __name__ == "__main__": + main() diff --git a/tests/test_protein.py b/tests/test_protein.py new file mode 100644 index 0000000..de3507c --- /dev/null +++ b/tests/test_protein.py @@ -0,0 +1,196 @@ +"""PDB files are punch cards, and that is the whole difficulty. + +The format is fixed columns, seventy years old, and the fields are not +separated by anything. An atom named "CA" in columns 13-16 is the alpha carbon +that every amino acid has in the middle of it; an atom whose ELEMENT is "CA", +in columns 77-78, is a calcium ion. Split the line on whitespace and those two +become the same thing, and the chain trace grows a spur out to wherever the +calcium happens to be sitting. + +So everything here is about reading the right columns, and about the records +that decide what is drawn: HELIX and SHEET, which are not deduced from the +geometry but written down by whoever solved the structure. + +No network. The PDB text is inline. +""" +import os +import sys + +import numpy as np +import pytest + +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "examples")) + +protein = pytest.importorskip("protein") + + +# Columns matter, so this is laid out by column and not by eye. +# 1 2 3 4 5 6 7 +# 1234567890123456789012345678901234567890123456789012345678901234567890123456789 +SAMPLE = """\ +TITLE A SMALL INVENTED PROTEIN +HELIX 1 1 ALA A 2 ALA A 4 1 3 +SHEET 1 A 2 ALA A 6 ALA A 7 0 +ATOM 1 N ALA A 1 10.000 10.000 10.000 1.00 0.00 N +ATOM 2 CA ALA A 1 11.000 10.000 10.000 1.00 0.00 C +ATOM 3 CA ALA A 2 14.800 10.000 10.000 1.00 0.00 C +ATOM 4 CA ALA A 3 18.600 10.000 10.000 1.00 0.00 C +ATOM 5 CA ALA A 4 22.400 10.000 10.000 1.00 0.00 C +ATOM 6 CA ALA A 6 26.200 10.000 10.000 1.00 0.00 C +ATOM 7 CA ALA A 7 30.000 10.000 10.000 1.00 0.00 C +ATOM 8 CA BALA A 8 99.000 99.000 99.000 0.50 0.00 C +ATOM 9 CA ALA B 1 10.000 20.000 10.000 1.00 0.00 C +ATOM 10 CA ALA B 2 13.800 20.000 10.000 1.00 0.00 C +HETATM 11 CA CA A 101 50.000 50.000 50.000 1.00 0.00 CA +HETATM 12 O HOH A 201 60.000 60.000 60.000 1.00 0.00 O +HETATM 13 FE HEM A 301 15.000 15.000 15.000 1.00 0.00 FE +END +""" + +NMR = """\ +MODEL 1 +ATOM 1 CA ALA A 1 0.000 0.000 0.000 1.00 0.00 C +ENDMDL +MODEL 2 +ATOM 2 CA ALA A 1 50.000 50.000 50.000 1.00 0.00 C +ENDMDL +""" + + +def _parse(**kw): + return protein.parse(SAMPLE, **kw) + + +# ── reading the columns ──────────────────────────────────────────────────── + +def test_an_alpha_carbon_is_not_a_calcium_ion(): + """The one that would quietly ruin the picture. + + Atom NAME "CA" in columns 13-16 is the alpha carbon of an amino acid. + ELEMENT "CA" in columns 77-78 is calcium. Both read as "CA" if the line is + split on spaces, and the calcium then gets threaded into the chain trace. + """ + atoms, _, _, _ = _parse() + calciums = [a for a in atoms if a["element"] == "CA"] + assert len(calciums) == 1 + assert calciums[0]["resname"] == "CA" and calciums[0]["hetatm"] + + chains = protein.backbone(atoms) + traced = [(a["chain"], a["resseq"]) for ch in chains for a in chains[ch]] + assert ("A", 101) not in traced, "the calcium ion got into the backbone" + + +def test_the_title_and_the_coordinates_come_out(): + atoms, _, _, title = _parse() + assert title == "A SMALL INVENTED PROTEIN" + first = atoms[0] + assert first["name"] == "N" and first["element"] == "N" + assert first["xyz"] == (10.0, 10.0, 10.0) + assert first["chain"] == "A" and first["resseq"] == 1 + + +def test_helix_and_sheet_ranges_are_inclusive(): + """HELIX 2-4 means 2, 3 and 4 are helix -- not 2 and 3.""" + _, helix, sheet, _ = _parse() + assert helix == {("A", 2), ("A", 3), ("A", 4)} + assert sheet == {("A", 6), ("A", 7)} + + +def test_only_one_conformer_of_a_disordered_side_chain(): + """altLoc B is a second position for the same atom, not another atom.""" + atoms, _, _, _ = _parse() + assert not any(a["resseq"] == 8 for a in atoms) + + +def test_water_is_dropped_and_ligands_are_kept(): + atoms, _, _, _ = _parse() + assert not any(a["resname"] == "HOH" for a in atoms) + assert any(a["element"] == "FE" for a in atoms), "the heme iron should stay" + + with_water, _, _, _ = _parse(waters=True) + assert any(a["resname"] == "HOH" for a in with_water) + + no_lig, _, _, _ = _parse(ligands=False) + assert not any(a["hetatm"] for a in no_lig) + + +def test_only_the_first_model_of_an_nmr_structure(): + """Twenty poses of one molecule drawn at once is a haystack.""" + atoms, _, _, _ = protein.parse(NMR) + assert len(atoms) == 1 + assert atoms[0]["xyz"] == (0.0, 0.0, 0.0) + + +# ── the chain ────────────────────────────────────────────────────────────── + +def test_the_backbone_is_split_by_chain_and_in_sequence(): + atoms, _, _, _ = _parse() + chains = protein.backbone(atoms) + assert sorted(chains) == ["A", "B"] + assert [a["resseq"] for a in chains["A"]] == [1, 2, 3, 4, 6, 7] + assert [a["resseq"] for a in chains["B"]] == [1, 2] + + +def test_the_chain_is_joined_up_with_no_gaps(): + """Blocks meeting only at an edge read as a broken chain.""" + for a, b in (((0, 0, 0), (5, 3, 2)), ((4, 4, 4), (0, 0, 0)), + ((0, 0, 0), (0, 0, 7)), ((2, 2, 2), (2, 2, 2))): + pts = protein.connected(a, b) + assert pts[0] == a and pts[-1] == b + for p, q in zip(pts, pts[1:]): + assert sum(abs(p[i] - q[i]) for i in range(3)) == 1 + + +def test_a_ball_is_round_and_the_right_size(): + for r in (0.5, 1.6, 3.0, 5.0): + off = protein.ball(r) + d = np.sqrt((off ** 2).sum(axis=1)) + assert d.max() <= r + 1.0 + assert (0, 0, 0) in {tuple(p) for p in off} + + +def test_scaling_keeps_the_shape_and_fills_the_size_asked_for(): + atoms = [{"xyz": (0.0, 0.0, 0.0)}, {"xyz": (10.0, 5.0, 2.0)}] + grid, scale, shape = protein.to_blocks(atoms, 100) + assert grid.min() == 0 + assert max(shape) == 101 # 0..100 inclusive + # proportions preserved: 10 : 5 : 2 stays 10 : 5 : 2 + span = grid.max(axis=0) - grid.min(axis=0) + assert span[0] == pytest.approx(2 * span[1], abs=1) + assert span[1] == pytest.approx(2.5 * span[2], abs=1) + + +# ── what gets built ──────────────────────────────────────────────────────── + +def test_a_backbone_is_coloured_by_secondary_structure(): + atoms, helix, sheet, _ = _parse() + world, palette, scale, n = protein.build_backbone( + atoms, helix, sheet, size=60, thickness=1.5, colour_by="structure") + assert n == 8 # six in chain A, two in B + assert set(palette) <= set(protein.STRUCTURE_BLOCKS.values()) + assert "RED_CONCRETE" in palette and "YELLOW_CONCRETE" in palette + assert int((world != 0).sum()) > 0 + + +def test_spacefill_is_coloured_by_element(): + atoms, _, _, _ = _parse() + world, palette, scale, n = protein.build_spacefill( + atoms, size=40, colour_by="element", chains_present=["A", "B"]) + assert n == len(atoms) + assert protein.ELEMENT_BLOCKS["N"] in palette # the nitrogen + assert protein.ELEMENT_BLOCKS["FE"] in palette # the heme iron + + +def test_nothing_built_here_falls_down(): + """These hang in mid-air, so a gravity-affected block would rain down.""" + falling = {"SAND", "RED_SAND", "GRAVEL"} + for name in list(protein.ELEMENT_BLOCKS.values()) + \ + list(protein.STRUCTURE_BLOCKS.values()) + \ + protein.RAINBOW + protein.CHAIN_BLOCKS + [protein.DEFAULT_ELEMENT]: + assert name not in falling + assert not name.endswith("_CONCRETE_POWDER"), name + + +def test_every_element_that_has_a_colour_has_a_radius(): + for element in protein.ELEMENT_BLOCKS: + assert element in protein.VDW, element