"""Generate a readable, explicitly provisional-or-complete scientific report.""" import collections import datetime import hashlib import json from pathlib import Path from score_ablation import uncertainty B = Path(__file__).resolve().parent E = B.parent SOURCES = {} def load(path): raw = path.read_bytes() SOURCES[str(path.relative_to(E))] = hashlib.sha256(raw).hexdigest() return json.loads(raw) def percent(n, d): return f"{n}/{d} ({100*n/d:.1f}%)" if d else "not yet available" def supported_baseline_summary(summary): """HLA*LA G classes cannot be truncated into numeric two-field calls.""" return {method: ({'G': resolutions['G']} if method == 'hlala' else resolutions) for method, resolutions in summary.items()} def interval_verdict(effect, interval, pvalue=None): if effect is None or interval is None: return "pending" if interval[0] > 0 and pvalue is not None and pvalue < .05: return "positive aggregate difference in this tested configuration" if interval[1] < 0 and pvalue is not None and pvalue < .05: return "negative aggregate difference in this tested configuration" return "no clear aggregate evidence under both family interval and sign-flip checks" def copy_contrast(rows, a, b, group_of): """Family-resampled recovery-rate difference; denominator is truth copies.""" index = {(r["donor"], r["gene"], r["arm"]): r for r in rows} pairs, recovered_a, recovered_b = [], 0, 0 for row in rows: if row["arm"] != a: continue other = index.get((row["donor"], row["gene"], b)) if other is None: continue left, right = row["unmatched"], other["unmatched"] n = left["eligible_copies"] assert n == right["eligible_copies"] delta = left["exact_copies_recovered"] - right["exact_copies_recovered"] assert abs(delta) <= n recovered_a += left["exact_copies_recovered"] recovered_b += right["exact_copies_recovered"] # Only family sums/denominators enter uncertainty(). This expansion # preserves both; it does not assert a phased identity for a gain. differences = [1 if delta > 0 else -1] * abs(delta) + [0] * (n - abs(delta)) pairs.extend({"donor": row["donor"], "difference": d} for d in differences) return {"a": a, "b": b, "eligible_copies": len(pairs), "recovered_a": recovered_a, "recovered_b": recovered_b, "donors_with_unnamed_truth": len({r["donor"] for r in pairs}), **uncertainty(pairs, group_of)} def main(): status_path = E / "execution-20261020/FINAL_RESULTS.json" status = load(status_path) if status_path.exists() else {"status": "in_progress"} track = load(E / "track-a-scoring-20261018/TRACK_A_RESULT.json") catalogue = load(B / "CATALOGUE_AUDIT.json") folds = load(E / "pangenome-reference-cv-20261014/FOLDS.json") group_of = {d: tuple(sorted(g)) for g in folds["groups"] for d in g} graph_planned = 0 for name in ("genomic-graph-rebase-20261020", "genomic-graph-rebase-remaining26-20261020"): plan_path, staging_path = E / name / "PLAN.json", E / name / "STAGING_PLAN.json" if plan_path.exists(): graph_planned += len(load(plan_path)["arms"]) elif staging_path.exists(): graph_planned += load(staging_path)["arms"] graph_planned = graph_planned or 120 lines = [ "DogoHLA: catalogue content and graph contribution", "Updated UTC: " + datetime.datetime.now(datetime.timezone.utc).isoformat(), "Execution status: " + status["status"], "", f"Scope: 86 Asian donors for catalogue tests; {graph_planned//2} for genomic graph-retrieval;", "60 for original-reference retrieval and 10%/25% read-fraction linkage.", "This is exposed developmental CV.", "An execution-complete result is not an untouched confirmatory benchmark.", "Counts below cover completed available outputs; compare methods through", "their matched contrasts, not through unequal availability denominators.", "Known family groups are resampled together. Donors absent from the", "relationship metadata are treated as singletons; independence is not certified.", "", "1. Catalogue typing", ] for arm in ["ipd", "genomic", "genomic_hprc", "genomic_hprc_asian", "genomic_nonasian", "genomic_core_hprc", "genomic_core_hprc_asian", "genomic_core_nonasian"]: r = track["summary"]["all"].get(arm) if not r: continue lines.append( f"{arm}: donors={r['donors']}; named={percent(r['correct'], r['eligible'])}; " f"unnamed recovered={r['unmatched_recovered']}/{r['unmatched_copies']}; " f"represented in training={r['in_panel_recovered']}/{r['in_panel']}.") primary = [] for r in track["contrasts"]: if r["cohort"] == "all" and (r["a"], r["b"]) in [ ("genomic", "ipd"), ("genomic_hprc_asian", "genomic"), ("genomic_core_hprc_asian", "genomic"), ("genomic_core_hprc_asian", "genomic_core_nonasian"), ("genomic_core_hprc_asian", "genomic_hprc_asian"), ] and r["pairs"]: primary.append(r) lines.append( f"Matched {r['a']} minus {r['b']}: {r['diff_pp']:+.2f} percentage points " f"(n={r['pairs']}; family bootstrap95={r['boot95']}; " f"family sign-flip p={r.get('family_sign_flip_p')}).") lines += [ "", "The core-trimmed diagnostic trims flanks and deduplicates sequences; the", "builder also skips exact IPD-identical sequences. It is a combined", "reference-normalization intervention, not an isolated flank experiment.", "Some retained cores have known genomic labels but different sequence", "boundaries; this is not a strictly unnamed-sequence-only augmentation.", "The pre-existing flanked contrast remains reported; core-trimmed results", "are a later diagnostic, not a retroactively chosen primary endpoint.", "", "2. Population content and its limits", "hprc = East/South Asian HPRC training donors; hprc_asian additionally", "includes Arab/Saudi donors. nonasian is matched by donor count, not by", "unique sequence count. These labels must not be conflated.", ] coverage = [r for r in catalogue["coverage_totals"] if r["gene"] == "all"] for r in coverage: lines.append( f"{r['region']}: {r['donors']} donors; {r['unnamed']} gene copies lack a genomic IPD match; " f"training-core coverage={r['unnamed_covered']}; " f"extra from Arab/Saudi={r['unnamed_extra_from_arab_saudi']}.") composition = collections.defaultdict(collections.Counter) for r in catalogue["composition"]: for key in ("unique_flanked_sequences", "unique_cores", "unique_unnamed_cores"): composition[r["panel"]][key] += r[key] for panel, counts in composition.items(): nfold = len({r["fold"] for r in catalogue["composition"] if r["panel"] == panel}) lines.append(f"{panel} mean per-fold catalogue size across eight loci: " + ", ".join(f"{k}={v/nfold:.1f}" for k, v in counts.items()) + ".") lines.append("This is an oracle coverage ceiling, not typing accuracy or certification of novel alleles.") copy_contrasts = [] for a, b in (("genomic_hprc", "genomic_nonasian"), ("genomic_hprc_asian", "genomic_nonasian"), ("genomic_hprc_asian", "genomic_hprc"), ("genomic_hprc_asian", "genomic"), ("genomic_core_hprc_asian", "genomic_core_nonasian"), ("genomic_core_hprc_asian", "genomic")): c = copy_contrast(track["rows"], a, b, group_of) copy_contrasts.append(c) if c["eligible_copies"]: lines.append( f"Matched unnamed-copy recovery, {a} minus {b}: " f"{c['recovered_a']} vs {c['recovered_b']} / {c['eligible_copies']}; " f"change={c['difference_pp']:+.2f}pp; CI95={c['bootstrap95_pp']}; " f"family p={c['sign_flip_p']:.5g}.") provenance_path = B / "CORE_PROVENANCE_AUDIT.json" if provenance_path.exists(): provenance = load(provenance_path) lines.append("Core-reference source provenance was reconstructed across merged contexts; " "no held-out donors were found in the reconstructed training sources.") sensitivity = {} for arm in ("genomic_hprc_asian", "genomic_core_hprc_asian", "genomic_core_nonasian"): selected = [r for r in track["rows"] if r["arm"] == arm and r["gene"] != "DRB1"] sensitivity[arm] = { "donors": len({r["donor"] for r in selected}), "unnamed": sum(r["unmatched"]["eligible_copies"] for r in selected), "recovered": sum(r["unmatched"]["exact_copies_recovered"] for r in selected), } s = sensitivity[arm] lines.append(f"Excluding DRB1, {arm}: unnamed recovery {s['recovered']}/{s['unnamed']} ({s['donors']} donors).") lines += ["", "3. Graph-proposed alignment contribution"] graph_reports = {} for filename in ("GRAPH_RETRIEVAL_ABLATION.json", "GRAPH_GENOMIC_RETRIEVAL_ABLATION.json"): path = B / filename if not path.exists(): lines.append(filename + ": pending") continue r = load(path) graph_reports[filename] = {"completed_donor_panel_arms": r["completed_donor_panel_arms"], "contrasts": r["contrasts"]} expected = graph_planned if filename.startswith("GRAPH_GENOMIC") else 120 lines.append(f"{filename}: {r['completed_donor_panel_arms']}/{expected} donor-panel arms.") for c in r["contrasts"]: if not c["eligible"]: continue lines.append( f"{c['panel']}/{c['mode']}: union {c['union_correct']}/{c['eligible']} vs native-only " f"{c['native_only_correct']}/{c['eligible']}; change={c['difference_pp']:+.2f}pp " f"CI95={c['bootstrap95_pp']}; p={c['sign_flip_p']:.5g}; " f"unnamed recovery {c['union_unnamed_recovered']} vs {c['native_only_unnamed_recovered']}.") graph_competitive = [] graph_path = B / "GRAPH_RESULTS.json" if graph_path.exists(): graph = load(graph_path) assert graph["track_result_sha256"] == SOURCES["track-a-scoring-20261018/TRACK_A_RESULT.json"] index = {(r["donor"], r["gene"], r["method"]): r for r in graph["rows"]} lines.append("Competitive comparison to full-genomic IPD (matched donors):") for contrast in graph["contrasts"]: if contrast["baseline"] != "T1K/genomic": continue method = contrast["method"] paired = [] for row in graph["rows"]: if row["method"] != method or not row["four_field"]["eligible"]: continue other = index.get((row["donor"], row["gene"], "T1K/genomic")) if other is None: continue assert other["four_field"]["eligible"] paired.append({"donor": row["donor"], "difference": int(row["four_field"]["correct"]) - int(other["four_field"]["correct"])}) if not paired: continue stats = {**contrast, **uncertainty(paired, group_of)} assert stats["pairs"] == len(paired) graph_competitive.append(stats) lines.append( f"{method} minus genomic: n={len(paired)}; gains={stats['gains']}, " f"losses={stats['losses']}; change={stats['difference_pp']:+.2f}pp; " f"CI95={stats['bootstrap95_pp']}; p={stats['sign_flip_p']:.5g}.") lines += [ "Both arms retain the same catalogue, fallback, rescoring model and gates.", "Only graph-proposed placements are removed in the native-only arm.", "A positive result identifies a contribution of this graph mapper within", "this pipeline; it does not show that every graph method is superior.", "The original-reference replay does not establish superiority to genomic IPD.", "", "4. Linkage at reduced read fractions", ] linkage = [] linkage_diagnostics = [] for tag in ("dsf100", "dsf250"): path = E / f"linkage-inference-{tag}-20261020/SCORES.json" if not path.exists(): lines.append(tag + ": pending") continue score = load(path) idx = {(r["donor"], r["panel"], r["mode"], r["gene"]): r for r in score["rows"]} collection_path = path.parent / "COLLECTION.json" proof = load(collection_path) diagnostics = collections.defaultdict(collections.Counter) for entry in proof["verified"]: arm = entry["arm"] folder = path.parent / "collected" / arm["donor"] / arm["panel"] receipt = load(folder / "RESULT.json") assert SOURCES[str((folder / "RESULT.json").relative_to(E))] == entry["result_sha256"] modes = {} for mode in ("independent", "linked", "shuffled"): modes[mode] = load(folder / (mode + ".json")) assert SOURCES[str((folder / (mode + ".json")).relative_to(E))] == receipt["output_sha256"][mode + ".json"] for baseline in ("independent", "shuffled"): stats = diagnostics[arm["panel"], baseline] stats["donors"] += 1 for gene, call in modes["linked"]["calls"].items(): other = modes[baseline]["calls"][gene] assert call["assigned_fragments"] == other["assigned_fragments"] stats["loci"] += 1 stats["selected_pair_changes"] += call["best_candidate_pairs"] != other["best_candidate_pairs"] stats["slot_changes"] += call["slots"] != other["slots"] stats["score_changes"] += abs(call["unphased_pair_score"] - other["unphased_pair_score"]) > 1e-12 left = idx[arm["donor"], arm["panel"], "linked", gene]["numeric_four_field"] right = idx[arm["donor"], arm["panel"], baseline, gene]["numeric_four_field"] if left["eligible"]: stats["eligible_loci"] += 1 stats["corrected"] += left["correct"] and not right["correct"] stats["harmed"] += right["correct"] and not left["correct"] for (panel, baseline), counts in diagnostics.items(): linkage_diagnostics.append({"tag": tag, "panel": panel, "baseline": baseline, **counts}) lines.append( f"{tag}/{panel} linked vs {baseline}: {counts['selected_pair_changes']}/{counts['loci']} " f"selected candidate pairs changed; eligible named calls corrected={counts['corrected']}, " f"harmed={counts['harmed']}.") for panel in ("hprc", "hprc_asian"): for baseline in ("independent", "shuffled"): paired = [] for r in score["rows"]: if r["panel"] != panel or r["mode"] != "linked": continue other = idx[r["donor"], panel, baseline, r["gene"]] assert r["numeric_four_field"]["eligible"] == other["numeric_four_field"]["eligible"] if r["numeric_four_field"]["eligible"]: paired.append({"donor": r["donor"], "difference": int(r["numeric_four_field"]["correct"]) - int(other["numeric_four_field"]["correct"])}) c = {"tag": tag, "panel": panel, "baseline": baseline, "status": score["status"], "donors": len({r["donor"] for r in paired}), "pairs": len(paired), **uncertainty(paired, group_of)} linkage.append(c) if paired: lines.append( f"{tag}/{panel} linked minus {baseline}: donors={c['donors']}; " f"n={c['pairs']}; change={c['difference_pp']:+.2f}pp; " f"CI95={c['bootstrap95_pp']}; p={c['sign_flip_p']:.5g}; status={score['status']}.") verdicts = [] for c in primary: verdicts.append({ "contrast": c["a"] + " minus " + c["b"], "verdict": interval_verdict(c["diff_pp"], c["boot95"], c.get("family_sign_flip_p")), "pairs": c["pairs"], "difference_pp": c["diff_pp"], "interval": c["boot95"], }) for filename, report in graph_reports.items(): for c in report["contrasts"]: if c["mode"] == "E3": verdicts.append({ "contrast": filename + " " + c["panel"] + "/E3 union minus native-only", "verdict": interval_verdict(c["difference_pp"], c["bootstrap95_pp"], c["sign_flip_p"]), "pairs": c["eligible"], "difference_pp": c["difference_pp"], "interval": c["bootstrap95_pp"], }) for c in linkage: verdicts.append({ "contrast": c["tag"] + " " + c["panel"] + " linked minus " + c["baseline"], "verdict": interval_verdict(c["difference_pp"], c["bootstrap95_pp"], c["sign_flip_p"]), "pairs": c["pairs"], "difference_pp": c["difference_pp"], "interval": c["bootstrap95_pp"], }) for c in copy_contrasts: verdicts.append({ "contrast": "unnamed recovery: " + c["a"] + " minus " + c["b"], "verdict": interval_verdict(c["difference_pp"], c["bootstrap95_pp"], c["sign_flip_p"]), "pairs": c["eligible_copies"], "difference_pp": c["difference_pp"], "interval": c["bootstrap95_pp"], }) for c in graph_competitive: if c["method"].startswith("DogoHLA_genomic/"): verdicts.append({ "contrast": c["method"] + " minus full-genomic IPD", "verdict": interval_verdict(c["difference_pp"], c["bootstrap95_pp"], c["sign_flip_p"]), "pairs": c["pairs"], "difference_pp": c["difference_pp"], "interval": c["bootstrap95_pp"], }) lines += ["", "Current statistical interpretation (provisional until execution is complete):"] for v in verdicts: lines.append(f"- {v['contrast']}: {v['verdict']} (n={v['pairs']}).") lines += [ "Read fractions are 10% and 25% of the original paired reads, not measured", "absolute coverage levels. Linked/independent/shuffled arms share the", "same read likelihoods and marginal candidate priors.", "", "5. Basic explanation and limits of any negative result", "- Full-genomic IPD repairs a major reference-completeness confound.", " Improvements over shortened default IPD must not be assigned to graphs.", "- Population haplotypes add relevant gene sequences missing from IPD.", " This can improve unnamed-sequence recovery even when named accuracy", " does not improve; the named endpoint excludes unnamed truth.", "- The current graph caller selects observed candidate paths. It cannot", " synthesize an allele absent from every candidate, even if some of its", " component variation is in the graph. Missing-catalogue failures are", " therefore distinct from mapping or inference failures.", "- Graph retrieval can help particular loci; aggregate gains must survive", " the native-only placement control and the strong genomic baseline.", "- Regional recruitment starts from linear-reference MHC, alternative", " and HLA alignments, with mate rescue. Graph inference cannot recover", " pairs discarded at recruitment. These tests therefore do not measure", " a possible benefit of graph-based recruitment from the full WGS.", "- A null LD result applies to the tested prior, read evidence and depths.", " It does not establish absence of biological linkage disequilibrium.", "- DRB1 dominates many unmatched counts. Exact annotation-core matches", " need independent raw-read/boundary validation before novelty claims.", "- All comparisons remain developmental. Multiple modes and diagnostic", " endpoints are reported rather than selecting the best result.", ] a3_path = E / "baseline-completion-20261020/A3_SCORES.json" a3_summary = None if a3_path.exists(): a3 = load(a3_path) supported = supported_baseline_summary(a3["summary"]) a3_summary = {"status": a3["status"], "summary": supported, "limitations": a3["limitations"]} lines += ["", "6. Supported-resolution comparator benchmarks", "OptiType/SpecHLA use two-field comparisons; HLA*LA uses common G classes.", "Unresolved numeric projections of HLA*LA G classes are not wrong two-field calls.", "SpecHLA uses 3.65 annotation with 3.38 mapping/binning; OptiType", "uses its bundled class-I database; HLA*LA uses its historical graph.", "They are not fully matched 3.65-database benchmarks.", "Detailed locus, region and stricter-resolution results: baseline-completion-20261020/A3_SCORES.txt."] for method, resolutions in supported.items(): endpoint = "G" if method == "hlala" else "2" label = "common-G" if endpoint == "G" else "two-field" c = resolutions[endpoint]["all"] lines.append( f"{method}: completed donors={c['completed_donors']}/86; " f"{label} correct={c['correct']}/{c['truth_eligible_pairs']}; " f"called={c['called']}; " f"matched with genomic T1K={c['paired_with_t1k']} eligible pairs.") if c["paired_with_t1k"]: u = c["contrast_comparator_minus_t1k"] lines.append( f" {method} minus genomic T1K: {u['difference_pp']:+.2f}pp; " f"CI95={u['bootstrap95_pp']}; family p={u['sign_flip_p']:.5g}.") lines.append("HLA-HD remains blocked by unavailable distribution; no access-request message was sent.") lines += ["", "Broader research plan disposition:"] broad = status.get("broader_limitations", {}) for name, disposition in broad.items(): lines.append(f"- {name}: {disposition}") report = { "status": status["status"], "catalogue_contrasts": primary, "coverage": coverage, "unnamed_copy_contrasts": copy_contrasts, "non_DRB1_sensitivity": sensitivity, "graph_retrieval": graph_reports, "graph_competitive": graph_competitive, "linkage": linkage, "linkage_diagnostics": linkage_diagnostics, "current_verdicts": verdicts, "supported_resolution_baselines": a3_summary, "broader_plan_complete": False, "source_sha256": SOURCES, "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "interpretation": lines[lines.index("5. Basic explanation and limits of any negative result") + 1:], } (B / "SCIENTIFIC_VERDICT.json").write_text(json.dumps(report, indent=2) + "\n") (B / "RESULTS.txt").write_text("\n".join(lines) + "\n") print(status["status"], "report updated") if __name__ == "__main__": main()