rugnux_vs_pdb.py: measure a corpus against the models its depositors published

A companion to rugnux_vs_xds.py, asking a different question. That one compares merging
statistics against another program's reduction of the same images; this one asks what an
independent party already decided the images were worth - the deposited model, its published
resolution and its published R-free - and compares against that.

The runs are deliberately bare: --model and nothing else. No space group, no cell, no resolution
limit. Anything told to the program is a thumb on the scale, and the question is what a bare
invocation reaches on its own.

The R-free it reports is a validation number rather than a refinement result: the deposited model
is placed rigid-body against our merge and the free residual of that placement is what is
compared. The depositor's published R-free came from full refinement against their own reduction,
so it is the harder number, and the comparison is deliberately unfair in their favour.

Entries are identified by directory name, which must be the PDB entry ID; metadata and
coordinates come from the RCSB API. The sweep is resumable, since a hundred datasets is hours.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N

docs/CHANGELOG.md

@@ -11,6 +11,7 @@
 * Where the cell metric carries more rotational symmetry than the Bravais class the indexer named, the extra rotations are put to the intensities and the space-group search is asked again on the metric's own cell - adopted only where the intensities confirm the higher symmetry - so a lattice that is nearly but not exactly hexagonal, or whose reduction landed in a sub-cell, still reaches its true point group.
 * Systematic-absence calls rest on the evidence rather than on counts: a screw axis whose absent class the data show extinct is no longer refused because a handful of reflections in it read as present, and `SPACE_GROUP_ALTERNATIVES=` no longer drops a candidate that differs only on a zone the sweep never measured.
 * A reference correlation measured on too few reflections is refused instead of scored zero, so a run given a reference MTZ is no longer reindexed on an operator that mapped almost everything outside the reference's coverage.
+* `rugnux_vs_pdb.py` processes a set of crystals whose data are in the PDB, validates each merge against the depositor's own model with `--model`, and tabulates and plots the measured resolution and R-free against the deposited ones.
 * A frame counts as indexed from 6 spots on its lattice rather than 9, so a weakly diffracting crystal whose frames cannot carry 9 is no longer refused the lattice it fits; `--min-indexed-spots` overrides it.
 * Where the two-pass rotation indexer chooses between the file's geometry and its own refined one, a refined pass that measures far fewer of the low-order axial reflections is rejected unless it is also better overall, so a crystal no longer loses a screw axis to a pass that merged marginally better while missing the reflections a space group is read from.
 * `-C` accepts a known cell in any equivalent description - conventional or primitive, centred or not - instead of only the reduced primitive form, so a centred cell given the way it is published no longer makes the run report that it found no lattice.

rugnux_vs_pdb.py

@@ -0,0 +1,271 @@
+#!/usr/bin/env python3
+"""
+rugnux_vs_pdb.py -- process a set of crystals whose data are deposited in the PDB, validate
+each merge against the depositor's own model, and compare the result with what they published.
+
+For every dataset given it:
+  * looks the entry up in the RCSB data API and records the deposited high-resolution limit,
+    R-work and R-free;
+  * downloads the deposited coordinates (gzipped mmCIF -- rugnux recognises the format by
+    content, so no conversion);
+  * runs rugnux de novo with --model and NOTHING else: no space group, no cell, no resolution
+    limit. That is the point. The question is what a bare invocation gets to, so anything told
+    to the program is a thumb on the scale;
+  * reads R_FREE / R_WORK and the resolution range back out of the results report, which is a
+    stable interface, and not out of the console log, which is not;
+  * tabulates deposited against measured, and plots the two scatters that answer the question:
+    resolution, and R-free.
+
+R-free here is not a refinement result. rugnux places the deposited model rigid-body and reports
+the free residual of THAT placement against its own merge, so it is a validation number: how well
+the depositor's model explains data this pipeline produced from their images. It is compared with
+the depositor's published R-free, which came from full refinement against their own reduction --
+a harder target, and deliberately so.
+
+Datasets are identified by directory name, which must be the four-character PDB entry ID.
+Small-molecule sets and anything without a deposited model are skipped with a reason.
+
+  ./rugnux_vs_pdb.py --rugnux build/rugnux/rugnux --data /data/raw --work /tmp/modelrun
+  ./rugnux_vs_pdb.py --work /tmp/modelrun --score-only        # re-score without re-running
+"""
+
+import argparse
+import json
+import os
+import re
+import subprocess
+import sys
+
+RCSB_ENTRY = "https://data.rcsb.org/rest/v1/core/entry/{}"
+RCSB_MODEL = "https://files.rcsb.org/download/{}.cif.gz"
+PDB_ID = re.compile(r"^[0-9][A-Za-z0-9]{3}$")
+
+
+def curl(url, dest=None, timeout=180):
+    cmd = ["curl", "-sfL", "--max-time", str(timeout), url]
+    if dest:
+        cmd += ["-o", dest]
+    r = subprocess.run(cmd, capture_output=True)
+    return r.returncode == 0, r.stdout
+
+
+def deposited(pdb, models_dir):
+    """Published resolution and R-factors, plus the coordinate file, for one entry."""
+    rec = {"id": pdb}
+    ok, body = curl(RCSB_ENTRY.format(pdb))
+    if ok:
+        try:
+            d = json.loads(body)
+            # The refine category is a list; a multi-stage refinement has more than one block and
+            # the last is the one the depositor stopped at.
+            refine = d.get("refine") or []
+            r = refine[-1] if refine else {}
+            rec["dep_dmin"] = r.get("ls_d_res_high")
+            rec["dep_rfree"] = r.get("ls_R_factor_R_free")
+            rec["dep_rwork"] = r.get("ls_R_factor_R_work")
+            rec["method"] = (d.get("exptl") or [{}])[0].get("method")
+        except Exception as e:
+            rec["meta_error"] = str(e)
+    else:
+        rec["meta_error"] = "entry not found"
+    dest = os.path.join(models_dir, f"{pdb}.cif.gz")
+    if not os.path.exists(dest) or os.path.getsize(dest) == 0:
+        ok, _ = curl(RCSB_MODEL.format(pdb), dest)
+        if not ok and os.path.exists(dest):
+            os.remove(dest)
+    if os.path.exists(dest) and os.path.getsize(dest) > 0:
+        rec["model"] = dest
+    else:
+        rec["model_error"] = "no deposited coordinates"
+    return rec
+
+
+def find_input(directory):
+    """The image file to hand rugnux: a master file wins, otherwise the first frame of the
+    lowest-numbered template, which is how the readers pick up a numbered series."""
+    best = None
+    for root, _, files in os.walk(directory):
+        for fn in sorted(files):
+            p = os.path.join(root, fn)
+            if fn.endswith(".h5") and "master" in fn:
+                return p
+            if fn.endswith(".nxs"):
+                best = best or p
+            elif best is None and (fn.endswith(".cbf") or fn.endswith(".h5")):
+                best = p
+    return best
+
+
+def parse_report(path):
+    rec = {}
+    with open(path, errors="replace") as fh:
+        for line in fh:
+            if "=" in line:
+                k, _, v = line.partition("=")
+                rec[k.strip()] = v.strip()
+    return rec
+
+
+def num(rec, key):
+    try:
+        return float(rec[key])
+    except (KeyError, TypeError, ValueError):
+        return None
+
+
+def run_one(rugnux, model, image, outdir, threads, gpulock):
+    os.makedirs(outdir, exist_ok=True)
+    cmd = []
+    if gpulock:
+        cmd += ["flock", gpulock]
+    cmd += [rugnux, "-o", "p", "-N", str(threads), "--no-export-unmerged",
+            "--spot-low-resolution", "50", "--scaling-low-resolution", "50",
+            "--model", model, image]
+    with open(os.path.join(outdir, "run.log"), "w") as log:
+        return subprocess.run(cmd, cwd=outdir, stdout=log, stderr=subprocess.STDOUT).returncode
+
+
+def collect(work, meta):
+    rows = []
+    for pdb in sorted(meta):
+        report = os.path.join(work, "runs", pdb, "p_report.txt")
+        if not os.path.exists(report):
+            continue
+        rep = parse_report(report)
+        rng = rep.get("INCLUDE_RESOLUTION_RANGE", "").split()
+        rows.append(dict(
+            id=pdb,
+            dep_dmin=meta[pdb].get("dep_dmin"),
+            dep_rfree=meta[pdb].get("dep_rfree"),
+            dep_rwork=meta[pdb].get("dep_rwork"),
+            our_dmin=float(rng[1]) if len(rng) == 2 else None,
+            our_rfree=num(rep, "R_FREE"),
+            our_rwork=num(rep, "R_WORK"),
+            model_fit=rep.get("MODEL_FIT"),
+            space_group=rep.get("SPACE_GROUP_NAME"),
+            cc_half=num(rep, "CC_HALF"),
+        ))
+    return rows
+
+
+def scatter(pts, xlabel, ylabel, title, path, floor_at_zero=False):
+    try:
+        import matplotlib
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+    except ImportError:
+        print(f"  (matplotlib absent, skipping {os.path.basename(path)})")
+        return
+    if not pts:
+        return
+    xs, ys = [p[0] for p in pts], [p[1] for p in pts]
+    lo = 0.0 if floor_at_zero else min(min(xs), min(ys)) * 0.95
+    hi = max(max(xs), max(ys)) * 1.05
+    fig, ax = plt.subplots(figsize=(6, 6))
+    ax.plot([lo, hi], [lo, hi], color="0.6", lw=1, zorder=1)
+    ax.scatter(xs, ys, s=26, alpha=0.8, edgecolor="none", zorder=2)
+    ax.set_xlim(lo, hi)
+    ax.set_ylim(lo, hi)
+    ax.set_aspect("equal")
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel(ylabel)
+    ax.set_title(title)
+    ax.grid(alpha=0.25, lw=0.5)
+    fig.tight_layout()
+    fig.savefig(path, dpi=140)
+    plt.close(fig)
+    print(f"  wrote {path}")
+
+
+def report(work, rows):
+    both = [r for r in rows if r["our_rfree"] is not None and r["dep_rfree"] is not None]
+    res = [r for r in rows if r["our_dmin"] is not None and r["dep_dmin"] is not None]
+    print(f"\n{len(rows)} processed; R-free on both sides for {len(both)}; resolution for {len(res)}\n")
+    if both:
+        n = len(both)
+        below = sum(1 for r in both if r["our_rfree"] <= r["dep_rfree"])
+        w02 = sum(1 for r in both if r["our_rfree"] <= r["dep_rfree"] + 0.02)
+        w05 = sum(1 for r in both if r["our_rfree"] <= r["dep_rfree"] + 0.05)
+        d = sorted(r["our_rfree"] - r["dep_rfree"] for r in both)
+        print(f"  R-free at or below deposited   {below:3d}/{n}  ({100 * below / n:.0f}%)")
+        print(f"  within +0.02                   {w02:3d}/{n}  ({100 * w02 / n:.0f}%)")
+        print(f"  within +0.05                   {w05:3d}/{n}  ({100 * w05 / n:.0f}%)")
+        print(f"  median delta R-free            {d[len(d) // 2]:+.4f}   ({d[0]:+.4f} .. {d[-1]:+.4f})")
+    if res:
+        deeper = sum(1 for r in res if r["our_dmin"] < r["dep_dmin"])
+        d = sorted(r["our_dmin"] - r["dep_dmin"] for r in res)
+        print(f"  d_min deeper than deposited    {deeper:3d}/{len(res)}  ({100 * deeper / len(res):.0f}%)")
+        print(f"  median delta d_min             {d[len(d) // 2]:+.3f} A  ({d[0]:+.3f} .. {d[-1]:+.3f})")
+    print()
+    scatter([(r["dep_dmin"], r["our_dmin"]) for r in res],
+            "deposited $d_{min}$ (A)", "rugnux $d_{min}$ (A)",
+            f"Resolution, n={len(res)} (below the line: rugnux cuts deeper)",
+            os.path.join(work, "dmin.png"))
+    scatter([(r["dep_rfree"], r["our_rfree"]) for r in both],
+            "deposited $R_{free}$", "rugnux $R_{free}$",
+            f"R-free against the deposited model, n={len(both)} (below the line: better)",
+            os.path.join(work, "rfree.png"), floor_at_zero=True)
+    no_rfree = [r for r in rows if r["our_rfree"] is None]
+    if no_rfree:
+        print(f"\nran but reported no R-free ({len(no_rfree)}):")
+        for r in no_rfree:
+            print(f"  {r['id']:<8}MODEL_FIT={r['model_fit']}")
+
+
+def main():
+    ap = argparse.ArgumentParser(
+        description="Process deposited datasets with --model and compare with the deposition.")
+    ap.add_argument("--work", required=True, help="working directory (models, runs, results)")
+    ap.add_argument("--rugnux", help="rugnux binary (required unless --score-only)")
+    ap.add_argument("--data", help="directory whose subdirectories are named by PDB entry ID")
+    ap.add_argument("--datasets", help="JSON file mapping entry ID to image file, instead of --data")
+    ap.add_argument("--threads", type=int, default=6)
+    ap.add_argument("--gpulock", help="lock file to serialise GPU use with other runs")
+    ap.add_argument("--score-only", action="store_true", help="re-score existing runs, run nothing")
+    args = ap.parse_args()
+
+    work = args.work
+    os.makedirs(os.path.join(work, "models"), exist_ok=True)
+    os.makedirs(os.path.join(work, "runs"), exist_ok=True)
+    meta_path = os.path.join(work, "deposited.json")
+    meta = json.load(open(meta_path)) if os.path.exists(meta_path) else {}
+
+    if not args.score_only:
+        if not args.rugnux:
+            ap.error("--rugnux is required unless --score-only")
+        if args.datasets:
+            inputs = json.load(open(args.datasets))
+        elif args.data:
+            inputs = {}
+            for name in sorted(os.listdir(args.data)):
+                d = os.path.join(args.data, name)
+                if os.path.isdir(d) and PDB_ID.match(name):
+                    image = find_input(d)
+                    if image:
+                        inputs[name.lower()] = image
+        else:
+            ap.error("one of --data or --datasets is required")
+
+        print(f"{len(inputs)} datasets with a PDB-shaped name and an image file")
+        for i, (pdb, image) in enumerate(sorted(inputs.items()), 1):
+            if pdb not in meta or "model" not in meta[pdb]:
+                meta[pdb] = deposited(pdb, os.path.join(work, "models"))
+                json.dump(meta, open(meta_path, "w"), indent=1)
+            if "model" not in meta[pdb]:
+                print(f"[{i}/{len(inputs)}] {pdb}: {meta[pdb].get('model_error', 'no model')}")
+                continue
+            outdir = os.path.join(work, "runs", pdb)
+            if os.path.exists(os.path.join(outdir, "p_report.txt")):
+                continue                      # resumable: an interrupted sweep picks up here
+            print(f"[{i}/{len(inputs)}] {pdb}", flush=True)
+            rc = run_one(args.rugnux, meta[pdb]["model"], image, outdir, args.threads, args.gpulock)
+            if rc != 0:
+                print(f"         exited {rc}")
+
+    rows = collect(work, meta)
+    json.dump(rows, open(os.path.join(work, "results.json"), "w"), indent=1)
+    report(work, rows)
+
+
+if __name__ == "__main__":
+    main()