From e9403943cbd53913e57227703a61e144762249fa Mon Sep 17 00:00:00 2001 From: Dave Lawrence Date: Tue, 18 Aug 2026 05:08:09 +0000 Subject: [PATCH] Paper: VariantValidator/Mutalyzer head-to-head and LOVD checker embedding note - #112 Extend the injection-corpus comparison to the two sequence-aware validation services over their public REST APIs (measured 2026-08-18, VariantValidator 4.0.1, Mutalyzer; script caches and resumes, gives up if a service blocks sustained automated use and reports coverage). Weighted by the production error mix VariantValidator recovers 56.1% and Mutalyzer 36.2% vs 100% for clean_hgvs(); neither falsely corrects a valid input. Table S5 gains the two service columns and rejection notes. Discussion and Methods note that VariantValidator 4.0.0 embeds the LOVD syntax checker's suggestions in its API responses without applying them, so the local LOVD measurement also characterises that layer. --- paper/Snakefile | 28 +- paper/discussion.md | 18 + .../vv_mutalyzer_comparison.csv | 2 + paper/methods.md | 18 +- paper/references.bib | 8 + paper/scripts/vv_mutalyzer_head_to_head.py | 550 ++++++++++++++++++ paper/supplementary.md | 85 ++- 7 files changed, 676 insertions(+), 33 deletions(-) create mode 100644 paper/empirical_results/vv_mutalyzer_comparison.csv create mode 100644 paper/scripts/vv_mutalyzer_head_to_head.py diff --git a/paper/Snakefile b/paper/Snakefile index 790b465..50f961a 100644 --- a/paper/Snakefile +++ b/paper/Snakefile @@ -75,7 +75,8 @@ GEN_FACTS = "output/facts" # inputs it needs are absent the rule copies the committed EMPIRICAL file instead, so # a build always renders; `full` writes these back to EMPIRICAL after a good render. COMPUTED_FACTS = ["coverage.csv", "benchmark.csv", "cleaning.csv", - "lovd_comparison.csv", "sources.csv", + "lovd_comparison.csv", "vv_mutalyzer_comparison.csv", + "sources.csv", "version_stability.csv", "positional_drift.csv"] # FROZEN - measured once (or transcribed from published papers) and committed under @@ -276,6 +277,31 @@ rule lovd_comparison: shutil.copy(f"{EMPIRICAL}/lovd_comparison.csv", output[0]) +rule vv_mutalyzer_comparison: + """Head-to-head: clean_hgvs() vs the VariantValidator and Mutalyzer public + REST services on the injection corpus (issue #112). + + Same cases and scoring as the lovd_comparison rule, driven by + paper/scripts/vv_mutalyzer_head_to_head.py over the public APIs + (rest.variantvalidator.org, mutalyzer.nl). Remote services cannot be + version-pinned, so the facts record the service date and the version + metadata the APIs report (measured 2026-08-18, VariantValidator 4.0.1). + The live run needs network access and about an hour (requests are cached + and resumable in output/rest_cache/); opt in with --config + vv_mutalyzer=true. Otherwise the frozen measured constants from + paper/empirical_results/ are copied. + """ + output: "output/facts/vv_mutalyzer_comparison.csv" + run: + import os, shutil + os.makedirs("output/facts", exist_ok=True) + if config.get("vv_mutalyzer", False): + shell(f"{PYTHON} paper/scripts/vv_mutalyzer_head_to_head.py") + else: + print("No vv_mutalyzer=true - copying the frozen measured service comparison facts") + shutil.copy(f"{EMPIRICAL}/vv_mutalyzer_comparison.csv", output[0]) + + rule sources: """Count annotation releases ingested per consortium/build (the sources.*_releases facts in Methods) from generate_transcript_data/cdot_transcripts.yaml - the diff --git a/paper/discussion.md b/paper/discussion.md index 380eac5..99dd774 100644 --- a/paper/discussion.md +++ b/paper/discussion.md @@ -71,6 +71,24 @@ descriptions on a transcript reference ({{ lovd_comparison.lovd_flagged_invalid_pct | dp(0) }}% of the valid originals) as requiring a genomic reference, a deliberate design position rather than a defect. +The same corpus puts the sequence-aware services in context. Run over their public REST +APIs ({{ vv_mutalyzer_comparison.service_date }}; Methods), VariantValidator's validated +description recovered {{ vv_mutalyzer_comparison.vv_weighted_pct | dp(0) }}% of cases +weighted by the production error mix, and Mutalyzer +{{ vv_mutalyzer_comparison.mut_weighted_pct | dp(0) }}% (Supplementary Table S5). Each +repairs the damage it can interpret and rejects the rest as invalid: whitespace and +letter case for both, with VariantValidator also handling quotes, doubled colons, +accession re-casing and swapped gene and transcript forms at 94 to 99% per category, +capped by transcripts absent from its database rather than by the injected error. +Neither falsely corrected a valid input. These services answer a different question from +`clean_hgvs()`: they validate a description against the reference sequence, so a string +that cannot be parsed is an error to report, not text to repair. The roles are +converging, though. Since release 4.0.0 (June 2026) VariantValidator embeds the LOVD +syntax checker [@VariantValidator400], and its API responses for input it cannot parse +carry the checker's ranked suggestions alongside the validation error, although the +service does not apply them. The LOVD measurement above therefore also characterises +the syntax-repair layer surfaced by a current VariantValidator pipeline. + Beyond HGVS resolution, the JSON representation is useful in its own right. It parses far faster than the GTF/GFF files it is built from and loads trivially over HTTP, so cdot doubles as a lightweight, queryable gene/transcript reference. We publish the per-release diff --git a/paper/empirical_results/vv_mutalyzer_comparison.csv b/paper/empirical_results/vv_mutalyzer_comparison.csv new file mode 100644 index 0000000..1c806fb --- /dev/null +++ b/paper/empirical_results/vv_mutalyzer_comparison.csv @@ -0,0 +1,2 @@ +service_date,vv_version,n_cases,vv_scored_cases,vv_coverage_pct,cdot_pct,vv_pct,mut_pct,cdot_weighted_pct,vv_weighted_pct,mut_weighted_pct,vv_refseq_pct,vv_ensembl_pct,mut_refseq_pct,mut_ensembl_pct,n_nonintronic,cdot_nonintronic_pct,vv_nonintronic_pct,mut_nonintronic_pct,originals_n,vv_originals_scored,vv_false_corrections,vv_rejected_valid,mut_false_corrections,mut_rejected_valid,mut_representation_changes,vv_lovd_suggestion_cases +2026-08-18,4.0.1.dev7+gbdab9c72f,3419,3419,100.0,100.0,37.7,12.5,100.0,56.1,36.2,32.9,44.2,11.8,13.6,2478,100.0,36.9,15.5,998,998,0,37,0,112,2,1528 diff --git a/paper/methods.md b/paper/methods.md index 61b8f19..3b32d71 100644 --- a/paper/methods.md +++ b/paper/methods.md @@ -236,7 +236,23 @@ list gave identical results, so top-1 ranking cost LOVD nothing. Every injected is string-repairable by construction (`inject_and_clean.py` does not inject errors whose repair needs transcript data, such as a missing accession prefix), so both tools are eligible on every case; the same false-correction check (does the tool alter a valid -input) is applied to both over the uncorrupted originals. Version-fallback safety is measured by `compute_version_stability.py` +input) is applied to both over the uncorrupted originals. + +`vv_mutalyzer_head_to_head.py` extends the same cases and scoring rule to the two +sequence-aware validation services, VariantValidator [@Freeman2018] and Mutalyzer +[@Lefter2021], over their public REST APIs +({{ vv_mutalyzer_comparison.service_date }}; remote services cannot be version-pinned, +so the facts record the service date and the version metadata each API reports). +VariantValidator requests are routed by accession family (its Ensembl transcripts live +on a separate endpoint) and a case counts as recovered when the validated transcript +description matches the target; for Mutalyzer a match by either its +`corrected_description` or its `normalized_description` counts, since the normalizer +may legitimately re-shift a representation. On the uncorrupted originals, a valid input +the service *alters* is a false correction, while one it *rejects* (for example +Mutalyzer's `EINTRONIC` for intronic positions on a transcript reference, or a +transcript absent from VariantValidator's database) is counted separately as a +validity or coverage position, matching how LOVD's flagged-invalid originals are +reported. Version-fallback safety is measured by `compute_version_stability.py` on GRCh38, using a seeded {{ version_stability.sample_n | commas }}-accession sample drawn from accessions cdot holds at two or more versions (the only accessions where a version bump can be assessed). The same run bins preserved coding bases by relative CDS position diff --git a/paper/references.bib b/paper/references.bib index 2f9a0a5..13fb9be 100644 --- a/paper/references.bib +++ b/paper/references.bib @@ -43,6 +43,14 @@ @misc{LovdHgvsChecker note = {Accessed 2026}, } +@misc{VariantValidator400, + author = {{VariantValidator}}, + title = {{VariantValidator} release 4.0.0: integration of the {LOVD} {HGVS} syntax checker}, + year = {2026}, + howpublished = {\url{https://github.com/openvar/variantValidator/releases/tag/v4.0.0}}, + note = {Released 5 June 2026}, +} + @article{Cingolani2012, author = {Cingolani, Pablo and Platts, Adrian and Wang, Le Lily and Coon, Melissa and Nguyen, Tung and Wang, Luan and Land, Susan J. and Lu, Xiangyi and Ruden, Douglas M.}, title = {A program for annotating and predicting the effects of single nucleotide polymorphisms, {SnpEff}}, diff --git a/paper/scripts/vv_mutalyzer_head_to_head.py b/paper/scripts/vv_mutalyzer_head_to_head.py new file mode 100644 index 0000000..d5f08c1 --- /dev/null +++ b/paper/scripts/vv_mutalyzer_head_to_head.py @@ -0,0 +1,550 @@ +#!/usr/bin/env python3 +""" +Head-to-head: cdot ``clean_hgvs()`` vs VariantValidator and Mutalyzer (issue #112). + +Extends ``lovd_head_to_head.py`` to the two sequence-aware validation services, +run over their public REST APIs on the exact injection corpus from +``inject_and_clean.py`` (same seed, caps and cases; see that script). Because +these are remote services there is no version pin; the service date and the +version metadata each API reports are recorded in the facts instead +(VariantValidator returns ``variantvalidator_version``; measured 2026-08-18). + +Services +-------- +- VariantValidator (rest.variantvalidator.org): batch pipe-delimited GET, + ``/VariantValidator/variantvalidator/GRCh38/{v1|v2|...}/all``. Ensembl + transcripts need the ``variantvalidator_ensembl`` endpoint, so requests are + routed by the accession family of the case's canonical target. Responses + for unparseable input embed the LOVD syntax checker's suggestions + (``lovd_corrections``, via api.lovd.nl), which is recorded per case. +- Mutalyzer (mutalyzer.nl): ``/api/normalize/{description}``, one call per + string, a few polite concurrent workers. + +Scoring (identical rule to lovd_head_to_head.py) +------------------------------------------------ +A case is recovered iff the service returns a description equal to the known +canonical target, compared gene-annotation-insensitively (see +``lovd_head_to_head.norm``). For VariantValidator the compared value is the +validated transcript description; for Mutalyzer either +``corrected_description`` or ``normalized_description`` may match (the +normalizer may legitimately re-shift a representation; accepting either is the +service's best answer). False-correction check on the uncorrupted originals: +the service returns a description differing from the valid input. + +Requests are cached to JSONL (``output/rest_cache/``) keyed by input string, +so an interrupted run resumes without re-querying; re-running only re-scores. +VariantValidator's public service blocks sustained automated use (it tarpits +the connection after a few hundred requests, which is how any scripted user +experiences it); after repeated dead batches the fetcher gives up, VV is +scored over the cases it answered, and the facts record that coverage +(``vv_scored_cases`` / ``vv_coverage_pct``). A later resume refetches the gap. + +Usage:: + + python paper/scripts/vv_mutalyzer_head_to_head.py # full run + python paper/scripts/vv_mutalyzer_head_to_head.py --limit 3 # smoke test +""" + +import argparse +import concurrent.futures +import csv +import json +import sys +import time +import urllib.error +import urllib.parse +import urllib.request +from pathlib import Path + +import hgvs.parser + +from cdot.hgvs.clean import clean_hgvs + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +import inject_and_clean as iac # noqa: E402 +from lovd_head_to_head import build_cases, norm # noqa: E402 + +FACTS_DIR = iac.FACTS_DIR +CACHE_DIR = iac.REPO / "output/rest_cache" +SERVICE_DATE = "2026-08-18" +USER_AGENT = "cdot-paper-benchmark (github.com/SACGF/cdot; davmlaw@gmail.com)" + +VV_BASE = "https://rest.variantvalidator.org/VariantValidator" +MUT_BASE = "https://mutalyzer.nl/api/normalize/" +VV_BATCH = 5 + + +def _get(url, timeout=300): + req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) + with urllib.request.urlopen(req, timeout=timeout) as r: + return json.loads(r.read()) + + +def _get_deadline(url, deadline): + """_get with a hard wall-clock deadline. Socket timeouts only bound + *inactivity*, so a tarpitting server that drip-feeds bytes can hold a + request open indefinitely; this abandons the request outright instead.""" + ex = concurrent.futures.ThreadPoolExecutor(1) + try: + return ex.submit(_get, url, deadline).result(timeout=deadline) + finally: + ex.shutdown(wait=False, cancel_futures=True) + + +# --------------------------------------------------------------------------- +# Cache +# --------------------------------------------------------------------------- + +class Cache: + def __init__(self, path): + self.path = path + self.data = {} + if path.exists(): + for line in path.read_text().splitlines(): + if line.strip(): + rec = json.loads(line) + if rec.get("error"): + continue # transient failure: refetch on resume + self.data[rec["input"]] = rec + path.parent.mkdir(parents=True, exist_ok=True) + self.fh = open(path, "a") + + def put(self, rec): + self.data[rec["input"]] = rec + self.fh.write(json.dumps(rec) + "\n") + self.fh.flush() + + +# --------------------------------------------------------------------------- +# VariantValidator driver +# --------------------------------------------------------------------------- + +def vv_endpoint(target): + fam = "variantvalidator_ensembl" if target.startswith("ENST") else "variantvalidator" + return f"{VV_BASE}/{fam}/GRCh38/" + + +def vv_parse(response, meta_out): + """Map submitted_variant -> record from a (possibly batch) VV response.""" + out = {} + meta_out.update(response.get("metadata") or {}) + for key, obj in response.items(): + if key in ("flag", "metadata") or not isinstance(obj, dict): + continue + sub = obj.get("submitted_variant") + if sub is None: + continue + rec = out.setdefault(sub, {"descriptions": [], "warnings": [], "lovd_suggestions": False}) + desc = obj.get("hgvs_transcript_variant") or "" + if desc: + rec["descriptions"].append(desc) + rec["warnings"].extend(obj.get("validation_warnings") or []) + if obj.get("lovd_corrections"): + rec["lovd_suggestions"] = True + return out + + +def _collapse(s): + return "".join(s.split()) + + +def vv_fetch_batch(cache, strings, endpoint, state): + """Fetch a batch. Timeouts and connection failures are throttling signals: + retry the same batch with growing backoff (never bisect, which multiplies + request load). Only an HTTP error response (a content problem) bisects, + down to singles, so one bad case cannot poison the batch. Three batches in + a row exhausting their retries abandons the VV fetch for this run (the + service is blocking sustained automated use; whatever is cached gets + scored, with coverage reported, and a later resume refetches the rest).""" + todo = [s for s in strings if s not in cache.data] + if not todo or state["abandoned"]: + return + url = (endpoint + urllib.parse.quote("|".join(todo), safe="") + + "/all?content-type=application%2Fjson") + err = None + for attempt in range(4): + try: + meta = {} + parsed = vv_parse(_get_deadline(url, deadline=120), meta) + # VV normalises whitespace before echoing submitted_variant, so + # fall back to a whitespace-collapsed lookup (batches never contain + # two strings sharing a collapsed form; vv_fetch_all guards that). + collapsed = {_collapse(k): v for k, v in parsed.items()} + for s in todo: + rec = parsed.get(s) or collapsed.get(_collapse(s)) or \ + {"descriptions": [], "warnings": ["NO_RESULT"], + "lovd_suggestions": False} + rec = dict(rec) + rec["input"] = s + rec["vv_version"] = meta.get("variantvalidator_version", "") + cache.put(rec) + state["consec_fail"] = 0 + time.sleep(2) # pacing between batch requests + return + except urllib.error.HTTPError as e: + err = e + break + except Exception as e: # timeout / tarpit / connection reset: back off + err = e + wait = min(600, 60 * 2 ** attempt) + print(f" VV throttled ({str(e)[:60] or type(e).__name__}); " + f"backing off {wait}s", flush=True) + time.sleep(wait) + else: + # Deadline exhaustion is a service-level block, not a content problem: + # bisecting would just multiply load. Leave the batch uncached and + # count the failure toward giving up. + state["consec_fail"] += 1 + if state["consec_fail"] >= 3: + state["abandoned"] = True + print(" VV: giving up for this run (service is blocking sustained " + "automated use); scoring the cases fetched so far", flush=True) + return + if len(todo) == 1: + cache.put({"input": todo[0], "descriptions": [], "warnings": [], + "lovd_suggestions": False, "error": str(err)[:200]}) + else: + mid = len(todo) // 2 + vv_fetch_batch(cache, todo[:mid], endpoint, state) + vv_fetch_batch(cache, todo[mid:], endpoint, state) + + +def vv_fetch_all(cache, items, log_every=20): + """items: list of (input_string, target). Routed and batched; a batch is + flushed early if it would contain a duplicate input or predicted target + (batch responses are keyed by corrected description, which must stay + unambiguous within one request).""" + groups = {} + for s, target in items: + groups.setdefault(vv_endpoint(target), []).append((s, target)) + n_req = 0 + state = {"consec_fail": 0, "abandoned": False} + for endpoint, pairs in groups.items(): + batch, seen = [], set() + pending = [p for p in pairs if p[0] not in cache.data] + for s, target in pending: + if state["abandoned"]: + return + key = norm(target) + if len(batch) >= VV_BATCH or _collapse(s) in seen or key in seen: + vv_fetch_batch(cache, batch, endpoint, state) + n_req += 1 + if n_req % log_every == 0: + done = sum(1 for p in pairs if p[0] in cache.data) + print(f" VV {endpoint.split('/')[-3]}: {done}/{len(pairs)}", + flush=True) + batch, seen = [], set() + batch.append(s) + seen.add(_collapse(s)) + seen.add(key) + if batch and not state["abandoned"]: + vv_fetch_batch(cache, batch, endpoint, state) + + +# --------------------------------------------------------------------------- +# Mutalyzer driver +# --------------------------------------------------------------------------- + +def mut_fetch_one(s): + url = MUT_BASE + urllib.parse.quote(s, safe="") + for attempt in range(3): + try: + d = _get(url, timeout=120) + return {"input": s, + "normalized": d.get("normalized_description", ""), + "corrected": d.get("corrected_description", ""), + "errors": []} + except urllib.error.HTTPError as e: + if e.code < 500: + try: + body = json.loads(e.read()) + body = body.get("custom", body) + codes = [x.get("code") for x in (body.get("errors") or [])] + except Exception: + codes = [f"HTTP{e.code}"] + return {"input": s, "normalized": "", "corrected": "", "errors": codes} + time.sleep(2 ** (attempt + 1)) + except Exception as e: + if attempt == 2: + return {"input": s, "normalized": "", "corrected": "", + "errors": [f"FAIL:{str(e)[:100]}"]} + time.sleep(2 ** (attempt + 1)) + return {"input": s, "normalized": "", "corrected": "", "errors": ["RETRIES"]} + + +def mut_fetch_all(cache, strings, workers=5, log_every=200): + todo = [s for s in dict.fromkeys(strings) if s not in cache.data] + with concurrent.futures.ThreadPoolExecutor(workers) as pool: + for i, rec in enumerate(pool.map(mut_fetch_one, todo)): + cache.put(rec) + if (i + 1) % log_every == 0: + print(f" Mutalyzer: {i + 1}/{len(todo)}", flush=True) + + +# --------------------------------------------------------------------------- +# Scoring +# --------------------------------------------------------------------------- + +def vv_hit(rec, target): + return any(norm(d) == norm(target) for d in rec.get("descriptions", [])) + + +def mut_hit(rec, target): + t = norm(target) + return norm(rec.get("corrected") or "") == t or norm(rec.get("normalized") or "") == t + + +def source_family(target): + return "ensembl" if target.startswith("ENST") else "refseq" + + +_INTRONIC = __import__("re").compile(r"\d+[+-]\d+") + + +def is_intronic(target): + """Intronic (or UTR-offset) position: Mutalyzer and LOVD both reject these + on a transcript reference by design (EINTRONIC / EWRONGREFERENCE).""" + return bool(_INTRONIC.search(target.split(":", 1)[-1])) + + +def score(cases, vv_cache, mut_cache): + per_class = {} + per_source = {"refseq": {}, "ensembl": {}} + for name, _ in iac.INJECTORS: + per_class[name] = {"weight": iac.REAL_RESCUE_OP_COUNTS[name], + "n_attempted": 0, "n_vv": 0, "cdot": 0, "vv": 0, + "mut": 0, "vv_lovd_suggestions": 0} + for src in per_source.values(): + src.update(n=0, n_vv=0, cdot=0, vv=0, mut=0) + nonintronic = {"n": 0, "n_vv": 0, "mut": 0, "vv": 0, "cdot": 0} + + for name, orig, perturbed, expected in cases: + c = per_class[name] + s = per_source[source_family(expected)] + c["n_attempted"] += 1 + s["n"] += 1 + if not is_intronic(expected): + nonintronic["n"] += 1 + cleaned, _ = clean_hgvs(perturbed) + if norm(cleaned) == norm(expected): + c["cdot"] += 1 + s["cdot"] += 1 + if not is_intronic(expected): + nonintronic["cdot"] += 1 + # VV is scored only over the cases it answered before blocking + # sustained automated use; coverage is reported alongside. + vrec = vv_cache.data.get(perturbed) + if vrec is not None: + c["n_vv"] += 1 + s["n_vv"] += 1 + if not is_intronic(expected): + nonintronic["n_vv"] += 1 + if vv_hit(vrec, expected): + c["vv"] += 1 + s["vv"] += 1 + if not is_intronic(expected): + nonintronic["vv"] += 1 + if vrec.get("lovd_suggestions"): + c["vv_lovd_suggestions"] += 1 + if mut_hit(mut_cache.data[perturbed], expected): + c["mut"] += 1 + s["mut"] += 1 + if not is_intronic(expected): + nonintronic["mut"] += 1 + + def pct(num, den): + return round(100.0 * num / den, 1) if den else None + + for c in per_class.values(): + c["cdot_pct"] = pct(c["cdot"], c["n_attempted"]) + c["mut_pct"] = pct(c["mut"], c["n_attempted"]) + c["vv_pct"] = pct(c["vv"], c["n_vv"]) + for s in per_source.values(): + s["cdot_pct"] = pct(s["cdot"], s["n"]) + s["mut_pct"] = pct(s["mut"], s["n"]) + s["vv_pct"] = pct(s["vv"], s["n_vv"]) + nonintronic["cdot_pct"] = pct(nonintronic["cdot"], nonintronic["n"]) + nonintronic["mut_pct"] = pct(nonintronic["mut"], nonintronic["n"]) + nonintronic["vv_pct"] = pct(nonintronic["vv"], nonintronic["n_vv"]) + return per_class, per_source, nonintronic + + +def weighted(per_class, key): + wsum = wnum = 0.0 + for c in per_class.values(): + if c[key + "_pct"] is not None: + wsum += c["weight"] + wnum += c["weight"] * c[key + "_pct"] + return round(wnum / wsum, 1) if wsum else None + + +def false_corrections(pool, vv_cache, mut_cache): + """On valid inputs, distinguish *altered* (service returns a different + description: a true false correction) from *rejected* (service returns + nothing: a validity design position or coverage gap, eg Mutalyzer's + EINTRONIC for intronic positions on a transcript reference, or a + transcript missing from VV's database).""" + counts = {"vv_scored": 0, "vv_altered": 0, "vv_rejected": 0, + "mut_altered": 0, "mut_rejected": 0, "mut_repr_changes": 0} + for orig in pool: + t = norm(orig) + vrec = vv_cache.data.get(orig) + if vrec is not None: + counts["vv_scored"] += 1 + if not vv_hit(vrec, orig): + if vrec.get("descriptions"): + counts["vv_altered"] += 1 + else: + counts["vv_rejected"] += 1 + mrec = mut_cache.data[orig] + if not mut_hit(mrec, orig): + if mrec.get("normalized") or mrec.get("corrected"): + counts["mut_altered"] += 1 + else: + counts["mut_rejected"] += 1 + elif norm(mrec.get("normalized") or "") != t: + counts["mut_repr_changes"] += 1 # accepted, representation re-shifted + return counts + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--limit", type=int, default=None, + help="cap cases per category (smoke test only; not for facts)") + ap.add_argument("--cache-dir", type=Path, default=CACHE_DIR) + args = ap.parse_args() + + parser = hgvs.parser.Parser() + pool = iac.load_sample(parser, regenerate=False) + cases = build_cases(parser, pool) + if args.limit: + by_cat = {} + cases = [c for c in cases + if by_cat.setdefault(c[0], []).append(c) or len(by_cat[c[0]]) <= args.limit] + pool = pool[: args.limit * 10] + + vv_cache = Cache(args.cache_dir / "vv.jsonl") + mut_cache = Cache(args.cache_dir / "mutalyzer.jsonl") + + items = [(perturbed, expected) for _, _, perturbed, expected in cases] + items += [(orig, orig) for orig in pool] + + t0 = time.time() + print(f"Mutalyzer: {sum(1 for s, _ in items if s not in mut_cache.data)} of " + f"{len(items)} strings uncached") + print(f"VariantValidator: {sum(1 for s, _ in items if s not in vv_cache.data)} of " + f"{len(items)} strings uncached") + # Mutalyzer first: it tolerates steady traffic, and the sequencing gives + # VariantValidator's anti-abuse tarpit a quiet period before VV resumes. + mut_fetch_all(mut_cache, [s for s, _ in items]) + vv_fetch_all(vv_cache, items) + fetch_minutes = round((time.time() - t0) / 60, 1) + + per_class, per_source, nonintronic = score(cases, vv_cache, mut_cache) + n = sum(c["n_attempted"] for c in per_class.values()) + n_vv = sum(c["n_vv"] for c in per_class.values()) + totals = {k: sum(c[k] for c in per_class.values()) for k in ("cdot", "vv", "mut")} + fc = false_corrections(pool, vv_cache, mut_cache) + vv_version = next((r.get("vv_version") for r in vv_cache.data.values() + if r.get("vv_version")), "") + + facts = { + "issue": "SACGF/cdot#112", + "tier": 1, + "description": ( + "Head-to-head on the reproducible injection corpus: cdot clean_hgvs() " + "vs VariantValidator and Mutalyzer public REST services, measured " + f"{SERVICE_DATE}; same cases and scoring as the LOVD comparison " + "(see paper/scripts/vv_mutalyzer_head_to_head.py)." + ), + "service_date": SERVICE_DATE, + "vv_version": vv_version, + "n_cases": n, + "sample_size": len(pool), + "vv_scored_cases": n_vv, + "vv_coverage_pct": round(100.0 * n_vv / n, 1) if n else None, + "cdot_pct": round(100.0 * totals["cdot"] / n, 1), + "vv_pct": round(100.0 * totals["vv"] / n_vv, 1) if n_vv else None, + "mut_pct": round(100.0 * totals["mut"] / n, 1), + "cdot_weighted_pct": weighted(per_class, "cdot"), + "vv_weighted_pct": weighted(per_class, "vv"), + "mut_weighted_pct": weighted(per_class, "mut"), + "vv_refseq_pct": per_source["refseq"]["vv_pct"], + "vv_ensembl_pct": per_source["ensembl"]["vv_pct"], + "mut_refseq_pct": per_source["refseq"]["mut_pct"], + "mut_ensembl_pct": per_source["ensembl"]["mut_pct"], + "n_nonintronic": nonintronic["n"], + "cdot_nonintronic_pct": nonintronic["cdot_pct"], + "vv_nonintronic_pct": nonintronic["vv_pct"], + "mut_nonintronic_pct": nonintronic["mut_pct"], + "originals_n": len(pool), + "vv_originals_scored": fc["vv_scored"], + "vv_false_corrections": fc["vv_altered"], + "vv_rejected_valid": fc["vv_rejected"], + "mut_false_corrections": fc["mut_altered"], + "mut_rejected_valid": fc["mut_rejected"], + "mut_representation_changes": fc["mut_repr_changes"], + "vv_lovd_suggestion_cases": sum(c["vv_lovd_suggestions"] for c in per_class.values()), + "fetch_minutes": fetch_minutes, + "per_class": per_class, + "per_source": per_source, + } + + FACTS_DIR.mkdir(parents=True, exist_ok=True) + json_path = FACTS_DIR / "vv_mutalyzer_comparison.json" + json_path.write_text(json.dumps(facts, indent=2) + "\n") + csv_path = FACTS_DIR / "vv_mutalyzer_comparison.csv" + row = {k: facts[k] for k in ( + "service_date", "vv_version", "n_cases", "vv_scored_cases", + "vv_coverage_pct", "cdot_pct", "vv_pct", "mut_pct", + "cdot_weighted_pct", "vv_weighted_pct", "mut_weighted_pct", + "vv_refseq_pct", "vv_ensembl_pct", "mut_refseq_pct", "mut_ensembl_pct", + "n_nonintronic", "cdot_nonintronic_pct", "vv_nonintronic_pct", + "mut_nonintronic_pct", "originals_n", "vv_originals_scored", + "vv_false_corrections", "vv_rejected_valid", "mut_false_corrections", + "mut_rejected_valid", "mut_representation_changes", + "vv_lovd_suggestion_cases")} + with open(csv_path, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=list(row)) + w.writeheader() + w.writerow(row) + + def s(x): + return "-" if x is None else str(x) + + print(f"\nVV/Mutalyzer head-to-head ({SERVICE_DATE}, {n} injected cases, " + f"VV answered {n_vv} ({facts['vv_coverage_pct']}%), " + f"fetch {fetch_minutes} min, VV {vv_version})") + print(f" {'category':<34} {'n':>4} {'cdot%':>7} {'n_vv':>5} {'vv%':>7} {'mut%':>7}") + for name, c in per_class.items(): + if not c["n_attempted"]: + continue + print(f" {name:<34} {c['n_attempted']:>4} {s(c['cdot_pct']):>7} " + f"{c['n_vv']:>5} {s(c['vv_pct']):>7} {s(c['mut_pct']):>7}") + print(f" {'-'*68}") + print(f" overall : cdot {s(facts['cdot_pct'])}% vv {s(facts['vv_pct'])}% " + f"(of answered) mut {s(facts['mut_pct'])}%") + print(f" weighted: cdot {s(facts['cdot_weighted_pct'])}% " + f"vv {s(facts['vv_weighted_pct'])}% mut {s(facts['mut_weighted_pct'])}%") + print(f" by source: vv refseq {s(facts['vv_refseq_pct'])}% ensembl " + f"{s(facts['vv_ensembl_pct'])}% | mut refseq {s(facts['mut_refseq_pct'])}% " + f"ensembl {s(facts['mut_ensembl_pct'])}%") + print(f" non-intronic ({nonintronic['n']}): cdot " + f"{s(facts['cdot_nonintronic_pct'])}% vv {s(facts['vv_nonintronic_pct'])}% " + f"mut {s(facts['mut_nonintronic_pct'])}%") + print(f" originals ({len(pool)}, vv answered {fc['vv_scored']}): " + f"vv altered {fc['vv_altered']} / rejected {fc['vv_rejected']}; " + f"mut altered {fc['mut_altered']} / rejected {fc['mut_rejected']} " + f"(+{fc['mut_repr_changes']} accepted with re-shifted representation)") + print(f" vv responses embedding LOVD suggestions: " + f"{facts['vv_lovd_suggestion_cases']}/{n_vv}") + print(f"\nWrote: {json_path}") + print(f"Wrote: {csv_path}") + + +if __name__ == "__main__": + main() diff --git a/paper/supplementary.md b/paper/supplementary.md index e095b8b..d7e8bb1 100644 --- a/paper/supplementary.md +++ b/paper/supplementary.md @@ -100,41 +100,50 @@ at most 200 cases per category), and recovery is an exact string match to the kn canonical target. The LOVD columns come from `paper/scripts/lovd_head_to_head.py`, which runs the same cases through the LOVD HGVS syntax checker (v1.2.2, local PHP CLI) under the same scoring rule (Methods); "top-1" scores the checker's highest-confidence -suggested correction. Production ops with no string-level injector (structure -reconstruction, empty-version dropping, provider-verified accession-prefix restoration) -are absent by design. Examples are synthesised from public `NM_000059.4` (BRCA2). -Regenerate with the commands in each script's docstring; totals also appear in -`paper/empirical_results/cleaning.csv` and `lovd_comparison.csv`. - -| Injected error (example → target) | n | `clean_hgvs()` | LOVD top-1 | -|---|---|---|---| -| Whitespace (` NM_000059.4: c.68del`) | 200 | 100% | 100% | -| Lowercased bases (`c.316g>a`) | 200 | 100% | 100% | -| Trailing protein suffix (`c.68del p.Arg100Ter`) | 200 | 100% | 91.5% | -| Gene wrapper with colon (`NM_000059.4:(BRCA2):c.68del`) | 200 | 100% | 0% | -| Gene/transcript swapped (`BRCA2(NM_000059.4):c.68del`) | 200 | 100% | 49.5% | -| Surrounding quotes (`"NM_000059.4:c.68del"`) | 200 | 100% | 0% | -| Doubled colon (`NM_000059.4::c.68del`) | 200 | 100% | 0% | -| Unbalanced bracket (`(NM_000059.4:c.68del`) | 200 | 100% | 0% | -| Separator typo (`NM_000059.4:c,68del`) | 200 | 100% | 100% | -| Doubled version dot (`NM_000059..4:c.68del`) | 200 | 100% | 0% | -| Leading junk (`GRCh38.p2 NM_000059.4:c.68del`) | 200 | 100% | 0% | -| Doubled kind (`NM_000059.4:c.c.68del`) | 200 | 100% | 0% | -| Lowercased accession (`nm_000059.4:c.68del`) | 200 | 100% | 52.0% | -| Redundant del/dup count (`c.68_69del23`) | 77 | 100% | 0% | -| Missing accession underscore (`NM000059.4:c.68del`) | 200 | 100% | 0% | -| Colon in accession prefix (`NM:_000059.4:c.68del`) | 200 | 100% | 0% | -| Uppercased mutation type (`c.68DEL`) | 142 | 100% | 100% | -| Dropped accession letter (`M_000059.4:c.68del`) | 200 | 100% | 0% | -| **Total** | **{{ lovd_comparison.n_cases | commas }}** | **{{ lovd_comparison.cdot_pct | dp(1) }}%** | **{{ lovd_comparison.lovd_top1_pct | dp(1) }}%** | +suggested correction. The VariantValidator (VV) and Mutalyzer columns come from +`paper/scripts/vv_mutalyzer_head_to_head.py`, which runs the same cases through the two +services' public REST APIs ({{ vv_mutalyzer_comparison.service_date }}; Methods) under +the same rule, scoring the validated (VV) or corrected/normalized (Mutalyzer) +description. Production ops with no string-level injector (structure reconstruction, +empty-version dropping, provider-verified accession-prefix restoration) are absent by +design. Examples are synthesised from public `NM_000059.4` (BRCA2). Regenerate with the +commands in each script's docstring; totals also appear in +`paper/empirical_results/cleaning.csv`, `lovd_comparison.csv` and +`vv_mutalyzer_comparison.csv`. + +| Injected error (example → target) | n | `clean_hgvs()` | LOVD top-1 | VV | Mutalyzer | +|---|---|---|---|---|---| +| Whitespace (` NM_000059.4: c.68del`) | 200 | 100% | 100% | 96.5% | 89.5% | +| Lowercased bases (`c.316g>a`) | 200 | 100% | 100% | 96.0% | 89.5% | +| Trailing protein suffix (`c.68del p.Arg100Ter`) | 200 | 100% | 91.5% | 0% | 0% | +| Gene wrapper with colon (`NM_000059.4:(BRCA2):c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Gene/transcript swapped (`BRCA2(NM_000059.4):c.68del`) | 200 | 100% | 49.5% | 94.5% | 0% | +| Surrounding quotes (`"NM_000059.4:c.68del"`) | 200 | 100% | 0% | 97.5% | 0% | +| Doubled colon (`NM_000059.4::c.68del`) | 200 | 100% | 0% | 99.0% | 0% | +| Unbalanced bracket (`(NM_000059.4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Separator typo (`NM_000059.4:c,68del`) | 200 | 100% | 100% | 0% | 0% | +| Doubled version dot (`NM_000059..4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Leading junk (`GRCh38.p2 NM_000059.4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Doubled kind (`NM_000059.4:c.c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Lowercased accession (`nm_000059.4:c.68del`) | 200 | 100% | 52.0% | 94.0% | 35.5% | +| Redundant del/dup count (`c.68_69del23`) | 77 | 100% | 0% | 0% | 0% | +| Missing accession underscore (`NM000059.4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Colon in accession prefix (`NM:_000059.4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| Uppercased mutation type (`c.68DEL`) | 142 | 100% | 100% | 94.4% | 0% | +| Dropped accession letter (`M_000059.4:c.68del`) | 200 | 100% | 0% | 0% | 0% | +| **Total** | **{{ lovd_comparison.n_cases | commas }}** | **{{ lovd_comparison.cdot_pct | dp(1) }}%** | **{{ lovd_comparison.lovd_top1_pct | dp(1) }}%** | **{{ vv_mutalyzer_comparison.vv_pct | dp(1) }}%** | **{{ vv_mutalyzer_comparison.mut_pct | dp(1) }}%** | Weighted by the production rescue-op distribution (Results Table 2) the totals are -{{ lovd_comparison.cdot_weighted_pct | dp(1) }}% for `clean_hgvs()` and -{{ lovd_comparison.lovd_top1_weighted_pct | dp(1) }}% for LOVD top-1; accepting the +{{ lovd_comparison.cdot_weighted_pct | dp(1) }}% for `clean_hgvs()`, +{{ lovd_comparison.lovd_top1_weighted_pct | dp(1) }}% for LOVD top-1, +{{ vv_mutalyzer_comparison.vv_weighted_pct | dp(1) }}% for VariantValidator and +{{ vv_mutalyzer_comparison.mut_weighted_pct | dp(1) }}% for Mutalyzer; accepting the target anywhere in LOVD's ranked correction list changes no case. On the -{{ lovd_comparison.originals_n | commas }} uncorrupted originals neither tool falsely +{{ lovd_comparison.originals_n | commas }} uncorrupted originals no tool falsely corrects any input ({{ lovd_comparison.cdot_false_corrections | int }} for `clean_hgvs()`, -{{ lovd_comparison.lovd_false_corrections | int }} for LOVD); LOVD flags +{{ lovd_comparison.lovd_false_corrections | int }} for LOVD, +{{ vv_mutalyzer_comparison.vv_false_corrections | int }} for VariantValidator, +{{ vv_mutalyzer_comparison.mut_false_corrections | int }} for Mutalyzer); LOVD flags {{ lovd_comparison.lovd_flagged_invalid_pct | dp(1) }}% of them (intronic positions on a transcript reference) as requiring a genomic reference while leaving the string unchanged. The LOVD partial rates are one-sided accession-family support: the @@ -143,6 +152,20 @@ re-casing for Ensembl but not RefSeq, and the two rates track the corpus's rough even RefSeq/Ensembl split. The protein-suffix misses absorb the stray `p` into the edit (`c.68del p.` becomes `c.68delP`). +The services' sub-100% ceilings in the categories they otherwise repair are +service-level rejections independent of the injected error, measured on the same +uncorrupted originals: VariantValidator rejects +{{ vv_mutalyzer_comparison.vv_rejected_valid | int }} of the +{{ vv_mutalyzer_comparison.originals_n | commas }} valid originals (all recent Ensembl +transcripts absent from its vvta_2025_02 database), and Mutalyzer rejects +{{ vv_mutalyzer_comparison.mut_rejected_valid | int }} (`EINTRONIC`: intronic positions +on a RefSeq transcript reference; it resolves intronic Ensembl descriptions). +Mutalyzer's lowercased-accession rate is one-sided in the opposite direction to LOVD's +(it re-cases `enst...` but not `nm_...`). VariantValidator's responses for input it +could not parse embedded the LOVD syntax checker's ranked suggestions in +{{ vv_mutalyzer_comparison.vv_lovd_suggestion_cases | commas }} of the +{{ vv_mutalyzer_comparison.n_cases | commas }} cases without applying them (Discussion). + ### Table S6: Residual error classes after cleaning **[Tier 2].** Single-label classification of the 1,075 production queries (826 unique