Paper: reframe around resolution goal, fair comparisons, terminology

claude/plans/new_analyses_plan.md

@@ -0,0 +1,191 @@
+# New analyses to do for the paper (deferred)
+
+From Dave's feedback (`claude/paper_feedbac_2026-08-19.md`, 2026-08-19). These need
+re-running scripts / new data work, so they are parked here. The prose/framing edits
+that do **not** need new data have already been applied to the paper.
+
+Constraint reminder: never run against full datasets (CLAUDE.md). Use samples / local
+UTA / committed test data. No `cdot_private` example strings in-repo.
+
+---
+
+## 1. Fair time-bucketed sampling of the ClinVar submitted (SCV) corpus  [R2]
+
+**Why:** the current submitted-corpus sample is a single fixed-seed flat draw over all
+unique (AlleleID, string) pairs, which has recency bias (ClinVar grows over time, so a
+flat draw over-represents recent submissions).
+
+**Do:**
+- In `build_clinvar_submitted_pairs.py`, add time-bucketed sampling: bucket SCV records
+  by submission date (or release/version era) and draw evenly across buckets, so
+  historical submissions are represented fairly.
+- Keep a **second** run that is a whole-file random sample, reported alongside and
+  explicitly labelled as recency-biased (reflects the live distribution).
+- Report both. The time-bucketed one is the fair "what labs submitted over the years"
+  number; the random one is "what the current file looks like".
+
+**Touches:** R2 numbers in `results.md`, `clinvar_submitted.csv`, Methods description of
+the sampling.
+
+**Code DONE 2026-08-19 (run pending):** `build_clinvar_submitted_pairs.py` now captures
+each pair's earliest SCV `SubmissionDate` from the VCV XML and writes it as a trailing
+`submit_date` corpus column (downstream `iter_vcf_pairs` reads by name, so it is
+ignored there). Two draws: `--sample-out` (the existing whole-file random draw,
+recency-biased) and `--time-bucketed-out` (new, even across submission-year eras,
+`--bucket-years` width, short eras drawn fully with the deficit redistributed). Both can
+be emitted from one build/`--from-pairs`. Validated on synthetic skewed corpora + a tiny
+synthetic VCV XML (earliest-date tracking, SCV collapse, even allocation, undated-drop,
+and the missing-column error path). STILL TO RUN (needs full VCV XML + VCF, a big job):
+rebuild the corpus with dates, emit both samples, resolve each via
+`resolve_clinvar_pass.py`, then update the R2 numbers + Methods prose to report both.
+Note: the committed `clinvar_submitted_500.tsv` predates the date column and should be
+regenerated when the corpus is rebuilt.
+
+**RUN DONE 2026-08-19 (results pending prose):** rebuilt from `ClinVarVCVRelease_2026-06`.
+New corpus = 3,198,528 unique (AlleleID, string) pairs, 100% RefSeq / 0 Ensembl, all with
+a submission date, 19 eras 2008-2026. (Corpus is larger than the committed 2,933,667: the
+XML path scans to the first *transcript* HGVS per assertion, capturing SCVs whose first
+HGVS attribute is genomic/protein, which the CSV-extraction path dropped. More complete
+and self-contained.) Raw distribution is heavily recency-skewed: 2024 = 1.14M, 2025 =
+1.21M pairs vs a few thousand per year in 2010-2013.
+
+cdot resolution (0.2.34 refseq GRCh38, FastaSeqFetcher, --with-fixes) on 3,000-pair
+seed-42 samples, VCF-coordinate scored:
+
+| | correct (after fix) | residual | 2008-2015 | 2016-2020 | 2021-2026 |
+|---|---|---|---|---|---|
+| RANDOM (recency-biased) | 2969 (99.0%) | 31 (1.0%) | 96.9% (n=32) | 99.2% | 99.0% |
+| TIME-BUCKETED (fair) | 2939 (98.0%) | 61 (2.0%) | 95.6% (n=1072) | 99.2% | 99.3% |
+
+The fair draw is 2x harder (2.0% vs 1.0% residual) because it up-weights old submissions:
+the random draw holds only 32 pre-2016 pairs (1%), the fair draw 1072 (36%), and pre-2016
+strings resolve at ~95.6% vs ~99% for recent ones. Residual composition also differs:
+random residual is mostly `no_data` (recent strings citing versions cdot lacks), fair
+residual is mostly `incorrect`/`error` (old versions projecting to drifted/invalid coords).
+Pass CSVs: scratchpad `random_pass_cdot.csv`, `bucketed_pass_cdot.csv`. Corpus + samples
+under `/data/clinvar/`.
+
+Dave's decisions 2026-08-19: ADOPT the 3.20M corpus; run UTA head-to-head on the samples
+(done); HOLD item 1.
+
+Full measured numbers (new 3.20M corpus, cdot 0.2.34 refseq GRCh38, local uta_20241220):
+
+cdot-vs-UTA head-to-head (VCF-scored):
+| draw | cdot correct | UTA correct | UTA no_data | cdot-only | uta-only |
+|---|---|---|---|---|---|
+| RANDOM (recency-biased) | 2969 (99.0%) | 2466 (82.2%) | 530 (17.7%) | 504 (16.8pts) | 1 |
+| TIME-BUCKETED (fair)    | 2939 (98.0%) | 2396 (79.9%) | 586 (19.5%) | 544 (18.1pts) | 1 |
+
+The fair draw WIDENS the cdot advantage (16.8 -> 18.1 pts): historical submissions cite
+more superseded versions UTA holds no GRCh38 alignment for.
+
+Version-age on the new corpus (`version_age_new.csv`): version_not_current 75.1% (was
+81.8% on the old 2.93M CSV-path corpus - the drop is the corpus redefinition),
+scv_weighted_not_current 69.8%, base_retired 0.7%, not_current_in_cdot 99.3%,
+absent_cdot_grch38 0.6%.
+
+STILL TODO (R2 prose/CSV rewrite): report BOTH samples (bucketed = fair headline, random
+= recency-biased); restructure `clinvar_submitted.csv` to carry both draws; pick which
+sample the abstract reports; recompute the residual taxonomy
+(`clinvar_submitted_residual.csv`) for the headline sample; regenerate committed
+`clinvar_submitted_500.tsv` from the new corpus. Pass CSVs in scratchpad:
+{random,bucketed}_pass_{cdot,uta}.csv.
+
+## PR #128 review comments (Dave, 2026-08-19) - addressed in parallel
+
+- DONE: renamed R2 heading "What laboratories submitted" -> "ClinVar submissions as a
+  historical record of transcripts used" (+ reworded lead-in).
+- DONE: turned the hardcoded R4 cleaning-corpus constants into a vibepaper fact
+  (`cleaning_corpus.csv`, wired into Snakefile FROZEN_FACTS); templated N, rescued,
+  rates, residual across results.md Table 2 + supplementary S6. Render verified.
+- DONE (draft, for Dave review): abstract motivation now leads with the Shariant origin
+  story. Intro already had the Shariant paragraph (lines 33-38).
+
+## 2. Run local UTA over (all of) ClinVar  [R2 / R3]
+
+**Why:** we currently only sample the cdot-vs-UTA comparison "because UTA is slow", but
+that caveat is for the *remote* server. The **local** `uta_20241220` may be able to do
+the whole ClinVar set in a few hours.
+
+**Do:** run the full ClinVar resolution through local UTA (background job), so the
+head-to-head is full-scale, not sampled. Removes the "gated to a sample by UTA
+throughput" caveat for the local case.
+
+**Touches:** R2 (drop/reduce the sampling caveat), possibly R3.
+
+## 3. Why is version substitution needed at all? (submitter attribution)  [new analysis]
+
+**Why:** we ingest every RefSeq/Ensembl release off the FTP sites, so in principle we
+should have every published version. Yet labs cite versions we don't hold. Dave's
+hypothesis: those transcripts come from a small number of labs that align transcripts
+to the genome themselves (so the version never appeared in an official annotation
+release).
+
+**Do (exploratory first, time-boxed):**
+- Take the cited transcript versions that are absent from cdot's data.
+- Group by submitter (SCV submitter / lab) in the ClinVar VCV XML.
+- Test whether the absent-version citations concentrate in a small number of submitters,
+  or whether specific labs dominate.
+- Characterise: are these self-aligned transcripts, pre-release versions, or genuine
+  gaps in our ingest?
+
+**Scope risk:** this is a genuinely new result and could balloon. Decide whether it
+earns a place in the paper (probably a short paragraph in R2 or Discussion) before
+investing heavily.
+
+## 4. Warm-cache benchmark rerun  [R3]
+
+**Why:** the full-scale R3 paragraph currently says the local-JSON and REST runs' "cache
+conditions differed, so their small difference is not evidence of REST outrunning local
+JSON." That is a hand-wave. Replace it with a hot-cache protocol.
+
+**Do:** for the full-scale runs, do one throwaway pass to warm the sequence-layer cache,
+discard it, then time the next pass. Report the hot-cache time and note it is only
+marginally slower than cold. Removes the "cache conditions differed" sentence entirely.
+
+**Touches:** R3 full-scale paragraph, `benchmark.csv` (if regenerated), Methods.
+
+## 5. Drop non-HGVS input from the cleaning residual corpus  [R4]  — DONE 2026-08-19
+
+**Why:** the R4 residual currently includes "non-HGVS input (pasted URLs, report
+templates, prose)". Dave's point: that is a data-collection artifact, not a measure of
+tool correctness, and should be removed from the denominator before counting the
+residual, since it does not reflect what the cleaner could ever be expected to fix.
+
+**Do:** filter non-HGVS input out of the production cleaning corpus before computing the
+residual rate, so the residual reflects genuine HGVS-repair ceiling only. Recompute the
+R4 residual number and the Table S6 taxonomy denominators.
+
+**Touches:** R4 residual paragraph + Table 2 context, `cleaning.csv`, Supplementary
+Table S6.
+
+**Done:** the R4/S6 residual and cleaning-rate numbers are literal frozen constants in
+the prose (not a facts CSV; `cleaning.csv` is the unrelated injection benchmark and was
+not touched). Kept the frozen LLM per-class taxonomy counts and moved the 81 non-HGVS
+queries out of the corpus. Reran `cdot_private/analyze_cleaning.py` to confirm the
+frozen baseline reproduces (as-is 29,956, after 31,676, residual 1,076 vs frozen 1,075,
+±1 code drift). Recomputed against N = 32,752 − 81 = 32,671:
+
+| quantity | was | now |
+|---|---|---|
+| corpus N | 32,752 | 32,671 |
+| parseable as-submitted | 91.5% | 91.7% (29,956/32,671) |
+| parseable after cleaning | 96.7% | 97.0% (31,677/32,671) |
+| absolute gain | +5.3% | +5.3% (1,721/32,671 = 5.27%) |
+| share of as-submitted failures rescued | ~62% of 8.5 pp | ~63% of 8.3 pp (1,721/2,715) |
+| residual | 3.3% (1,075 q) | 3.0% (994 q) |
+| S6 class %s | of 1,075 | of 994 (counts unchanged) |
+
+The "826 unique strings" figure was dropped (the unique count *within* the 81 non-HGVS
+needs the lost per-string LLM labels; residual is now reported in queries only).
+
+Edited: `results.md` (main R4 para, Table 2 header, residual para), `supplementary.md`
+(S6 header + table %s), `paper/README.md` (Tier-2 rate note). Paper re-rendered OK.
+
+---
+
+## Also consider (from feedback, smaller)
+
+- ~~**VEP + RefSeq alignment gaps citation:**~~ DONE. R1 now cites Ensembl/ensembl-vep
+  issue #1053 (`@VepHgvsGaps` in references.bib) for the claim that VEP does not apply
+  transcript-to-genome alignment gaps when converting c. to genomic coordinates.

paper/README.md

@@ -82,7 +82,7 @@ snakemake -s paper/Snakefile full --config pdf=true --cores 1 \
 
 ### Full build — regenerates every fact, then renders
 
-Recomputes the reproducible (Tier-1) facts from scratch, refreshes the
+Recomputes the reproducible (public-data) facts from scratch, refreshes the
 `paper/empirical_results/` snapshot, then renders. Slow, and needs the inputs each fact requires.
 
 ```bash
@@ -103,7 +103,7 @@ needs are absent, it copies the committed CSV instead, so a build always renders
 | `coverage.csv` | `paper/scripts/compute_coverage.py` over the release JSON.gz files | `data_dir` |
 | `benchmark.csv` | `paper/scripts/compute_benchmark.py` (Table 1 throughput) | `data_dir`, `uta_uri`, SeqRepo |
 | `version_stability.csv`, `positional_drift.csv` | `paper/scripts/compute_version_stability.py` | `data_dir` |
-| `cleaning.csv` | `paper/scripts/inject_and_clean.py` (Tier-1 injection benchmark) | committed test data (none) |
+| `cleaning.csv` | `paper/scripts/inject_and_clean.py` (public-data injection benchmark) | committed test data (none) |
 | `lovd_comparison.csv` | `paper/scripts/lovd_head_to_head.py` | `lovd_checker` (local PHP CLI) |
 | `sources.csv` | `paper/scripts/compute_sources.py` over `cdot_transcripts.yaml` | committed (none) |
 
@@ -119,17 +119,17 @@ copies them into the facts dir. To change one, edit the CSV and update its entry
 
 The full-scale ClinVar throughput runs take ~1.5 h each — see `claude/benchmark_plan.md`.
 
-### Two-tier facts and the private corpus
+### Public and private facts
 
-- **Tier 1 (reproducible)** lives in the fact CSVs above and regenerates from public
-  data committed here.
-- **Tier 2 (production validation, not reproducible)** — the cleaning rescue rate
-  (91.5% → 96.7%), the per-fix rescue distribution (Results Table 2), and the residual
-  error taxonomy — comes from the private `cdot_private` corpus and is written into
+- **Public data (reproducible)** lives in the fact CSVs above and regenerates from
+  public data committed here.
+- **Private data (production validation, not reproducible):** the cleaning rescue rate
+  (91.7% → 97.0%), the per-fix rescue distribution (Results Table 2), and the residual
+  error taxonomy come from the private `cdot_private` corpus and are written into
   `results.md` as **literal frozen constants**, not regenerable facts. No corpus string
   ever enters this repo. When the corpus is re-analysed (`cdot_private/analyze_cleaning.py`),
   transcribe the new aggregate numbers into `results.md` by hand. These are flagged
-  `[Tier 2]` in the text.
+  `[private data]` in the text.
 
 To refresh the committed snapshot without a full data run (e.g. after editing one
 analysis script), regenerate that one CSV into `output/facts/` and copy it into

paper/Snakefile

@@ -88,7 +88,8 @@ FROZEN_FACTS = ["literature.csv",
                 "clinvar_submitted.csv", "clinvar_submitted_residual.csv",
                 "clinvar_vcf.csv", "clinvar_vcf_residual.csv",
                 "clinvar_residual_positions.csv", "genomic_mismatch.csv",
-                "historical.csv", "version_safety_validation.csv"]
+                "historical.csv", "version_safety_validation.csv",
+                "cleaning_corpus.csv"]
 
 FACT_FILES = COMPUTED_FACTS + FROZEN_FACTS
 

paper/abstract.md

@@ -4,12 +4,16 @@
 
 **Motivation:**
 
-Resolving a transcript-level HGVS variant description to genomic coordinates requires
-accurate, versioned transcript-to-genome alignments. The standard source for this,
-UTA, requires a PostgreSQL database, covers only
-~{{ literature.uta_count | commas }} alignments, retains limited transcript history,
-and omits Ensembl, while the descriptions real pipelines receive are often malformed or
-cite retired transcript versions.
+cdot was built for the Australian Genomics Shariant project [@Tudini2022], which pools
+variant classifications from clinical laboratories across the country, each recorded
+against whatever transcript version the submitting lab used at the time. Making that
+shared history usable meant resolving as many real-world transcript-level HGVS
+descriptions as possible, including the malformed strings and long-retired transcript
+versions that accumulate in variant databases and clinical records. The standard
+transcript source for the Python HGVS libraries, UTA, forces a tradeoff between a
+locally installed PostgreSQL database and a slow public server that many clinical
+networks firewall; it covers only ~{{ literature.uta_count | commas }} RefSeq
+alignments, retains limited transcript history, and omits Ensembl.
 
 **Results:**
 
@@ -21,16 +25,18 @@ the Python HGVS libraries, T2T-CHM13v2.0, as a single gzipped JSON file or a RES
 version that is no longer current, cdot resolved
 {{ clinvar_submitted.cdot_resolved_pct | dp(1) }}% versus
 {{ clinvar_submitted.uta_resolved_pct | dp(1) }}% for UTA. Loaded in memory, cdot
-resolves ~{{ benchmark.cdot_local_tps | commas }} HGVS/second, nearly four orders of
-magnitude above the public UTA server. A parser-independent cleaning step
+resolves ~{{ benchmark.cdot_local_tps | commas }} HGVS/second, about four times a
+locally installed UTA, and compared remote-to-remote its REST API is roughly two orders
+of magnitude faster than the public UTA server. A parser-independent cleaning step
 (`clean_hgvs()`) repairs common formatting errors before resolution, and an opt-in
-fallback substitutes a retired transcript version only when a coordinate-safety check
-confirms the substitution does not move the variant.
+version-substitution step supplies a retired transcript version only when a
+coordinate-safety check confirms it does not change the variant's coordinate.
 
 **Availability and Implementation:**
 
-https://github.com/SACGF/cdot; `pip install cdot`; MIT licence. Data files (JSON.gz) at
-cdotlib.org and [Zenodo DOI].
+https://github.com/SACGF/cdot; `pip install cdot`; MIT licence. Data files (JSON.gz) are
+published on the GitHub releases page and archived at [Zenodo DOI]; a REST API is served
+at cdotlib.org.
 
 **Contact:** [email]