diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index ec9c45588..328848751 100644 --- a/docs/CHANGELOG.md +++ b/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. diff --git a/rugnux_vs_pdb.py b/rugnux_vs_pdb.py new file mode 100755 index 000000000..3de073204 --- /dev/null +++ b/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()