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 948bc2d..bb1f54d 100644 --- a/paper/discussion.md +++ b/paper/discussion.md @@ -39,6 +39,24 @@ error mix, LOVD's top-ranked correction restored common single-token errors but not on free-text damage or cross-accession-family repairs; neither tool altered any valid input. +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 the per-release JSON published on the GitHub releases page is a faster-loading drop-in for the corresponding GTF/GFF. 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 eac588b..79b35a8 100644 --- a/paper/methods.md +++ b/paper/methods.md @@ -152,7 +152,23 @@ clean ClinVar strings. `lovd_head_to_head.py` runs the same injected cases throu ({{ lovd_comparison.lovd_version }}, run locally as a PHP CLI), scoring a case as recovered when the tool's output (for LOVD, its top-ranked correction) exactly matches the canonical target; the same false-correction check runs over the uncorrupted -originals. Version-fallback safety is measured by `compute_version_stability.py` on +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, over a seeded {{ version_stability.sample_n | commas }}-accession sample of accessions cdot holds at two or more versions; the same run bins preserved coding bases by relative CDS position (Supplementary Figure S1). 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 ad0e745..2c9913f 100644 --- a/paper/supplementary.md +++ b/paper/supplementary.md @@ -185,41 +185,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 @@ -228,6 +237,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