diff --git a/claude/plans/new_analyses_plan.md b/claude/plans/new_analyses_plan.md index 5be2dc6..a9b3284 100644 --- a/claude/plans/new_analyses_plan.md +++ b/claude/plans/new_analyses_plan.md @@ -133,6 +133,42 @@ release). earns a place in the paper (probably a short paragraph in R2 or Discussion) before investing heavily. +**Code DONE 2026-08-19 (Dave: "just do everything"), run pending.** `submitter_attribution.py` +streams the VCV XML, attributes each SCV's first RefSeq transcript citation to its +submitter (`ClinVarAccession@SubmitterName`, else `ClinVarSubmissionID@submitter`), and +flags absent-from-cdot versions using the same version-set test as +`compute_submitted_version_age.py`. Reports concentration (top-N submitter share of +absent citations) and the self-alignment signature (absent versions cited by a single +submitter). Validated on a synthetic XML. Full run deferred until item 1's timed pass +finishes (heavy parse would confound item 1 throughput); will run alongside item 4. + +## Item 4 RESULT (full-scale, 2026-08-21) + +Ran the full submitted corpus (3,198,528 pairs) through cdot (--fasta --with-fixes, ~20h +single-core) and local UTA (SeqRepo, ~7.5h). Full head-to-head: + +| | cdot | UTA | +|---|---|---| +| matched VCF coordinate | 99.0% (3,166,066) | 81.9% (2,618,768) | +| resolves through this backend alone | 548,524 (17.2 pts) | 1,226 | +| no_data | 19,124 | 573,994 (18.0%) | + +Per-era: 2008-2015 cdot 97.3% / UTA 76.2%; 2016-2020 99.3% / 88.8%; 2021-2026 99.0% / +81.0%. Residual 1.0% (32,462), 0 regressions. The full-scale numbers land near the +recency-weighted random draw; the per-era split is the fair view without sampling. +R2 rewritten to full-scale + per-era, dropping the two-sample framing (superseded). +Note: cdot --fasta is the bottleneck (~43 HGVS/s, historical versions force genome +reconstruction); UTA-full was actually the fast part. + +## Status 2026-08-19 (PR #128 MERGED; items 1/4/5 on new branch paper-benchmarks-2026-08-19) + +- Items 2, 3 + all PR #128 feedback: merged to main. +- Item 1 (warm-cache) + Item 4 (full UTA): RUNNING (sequential background runner PID + 357490, monitor b4g9mry9g). Item 1 first (~4h), then item 4 cdot-full + UTA-full (~15h). +- Item 5: code ready, run queued for after item 1's timed pass. +- When results land: update R3 (hot-cache throughput), R2 (full-scale head-to-head), add + item 5 paragraph; open new PR from paper-benchmarks-2026-08-19. + ## 4. Warm-cache benchmark rerun [R3] **Why:** the full-scale R3 paragraph currently says the local-JSON and REST runs' "cache diff --git a/paper/Snakefile b/paper/Snakefile index c11fa8f..c4aa0cb 100644 --- a/paper/Snakefile +++ b/paper/Snakefile @@ -89,7 +89,8 @@ FROZEN_FACTS = ["literature.csv", "clinvar_vcf.csv", "clinvar_vcf_residual.csv", "clinvar_residual_positions.csv", "genomic_mismatch.csv", "historical.csv", "version_safety_validation.csv", - "cleaning_corpus.csv"] + "cleaning_corpus.csv", "benchmark_fullscale.csv", + "submitter_attribution.csv"] FACT_FILES = COMPUTED_FACTS + FROZEN_FACTS diff --git a/paper/empirical_results/PROVENANCE.md b/paper/empirical_results/PROVENANCE.md index 2d6e821..4356c0e 100644 --- a/paper/empirical_results/PROVENANCE.md +++ b/paper/empirical_results/PROVENANCE.md @@ -87,14 +87,14 @@ the earliest SCV submission date so it can be sampled by era. The corpus (built ClinVar XML download) is not committed, but a 500-pair seed-42 random sample is (`tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv`). -Two 3,000-pair samples are scored (seed 42): a whole-file **random** draw (reflects the -live file, recency-biased) and a **time-bucketed** draw (even across submission-year -eras 2008-2026, the fair historical picture). R2 leads with the fair sample; the random -draw's numbers are the `rnd_*` columns. The fair draw is harder and widens the cdot-UTA -gap because it up-weights old submissions citing superseded versions. +The cdot-vs-UTA head-to-head is run over the **whole** corpus (no sampling) and broken +down by submission-year era (`era_*` columns); a committed 500-pair random sample +supports quick checks. The full run supersedes the earlier two-sample approach (a +time-bucketed and a random draw); `build_clinvar_submitted_pairs.py` still emits those +draws (`--sample-out` / `--time-bucketed-out`) but R2 now reports the full corpus. -Measured 2026-08-19 (ClinVarVCVRelease_2026-06, cdot 0.2.34 refseq GRCh38, local -uta_20241220, SeqRepo for UTA sequence): +Measured 2026-08-20/21 (ClinVarVCVRelease_2026-06, cdot 0.2.34 refseq GRCh38, local +uta_20241220; FastaSeqFetcher for cdot, SeqRepo for UTA sequence): ```bash # corpus: 4,027,987 SCV transcript expressions -> 3,198,528 unique (AlleleID, string) @@ -103,44 +103,70 @@ uta_20241220, SeqRepo for UTA sequence): # first HGVS attribute is genomic/protein (the --scv-csv-dir path dropped those). python paper/scripts/build_clinvar_submitted_pairs.py \ --xml ClinVarVCVRelease_2026-06.xml.gz clinvar.GRCh38.vcf.gz \ - clinvar_submitted_pairs_dated.GRCh38.tsv --sample 3000 --seed 42 \ - --sample-out submitted_random_3000.tsv \ - --time-bucketed-out submitted_bucketed_3000.tsv + clinvar_submitted_pairs_dated.GRCh38.tsv # version age vs the current annotation release (RS_2025_08): 75.13% not-current -python paper/scripts/compute_submitted_version_age.py \ - clinvar_submitted_pairs_dated.GRCh38.tsv \ - --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz \ - --refseq-allbuilds cdot-0.2.34.all-builds-refseq-....json.gz +python paper/scripts/compute_submitted_version_age.py clinvar_submitted_pairs_dated.GRCh38.tsv \ + --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz --refseq-allbuilds ... -# resolution on each sample (VCF-coordinate scoring), cdot then UTA: +# full-corpus resolution (VCF-coordinate scoring), cdot (~20 h) then UTA (~7.5 h): F=GCF_000001405.39_GRCh38.p13_genomic.fna.gz -for S in random bucketed; do - python paper/scripts/resolve_clinvar_pass.py submitted_${S}_3000.tsv \ - --json cdot-0.2.34.refseq.GRCh38.json.gz --fasta $F --with-fixes \ - --out ${S}_pass_cdot.csv - UTA_DB_URL=postgresql://postgres@127.0.0.1:5433/uta/uta_20241220 \ - HGVS_SEQREPO_DIR=... python paper/scripts/resolve_clinvar_pass.py \ - submitted_${S}_3000.tsv --uta --out ${S}_pass_uta.csv -done +python paper/scripts/resolve_clinvar_pass.py clinvar_submitted_pairs_dated.GRCh38.tsv \ + --json cdot-0.2.34.refseq.GRCh38.json.gz --fasta $F --with-fixes --out submitted_full_cdot.csv +UTA_DB_URL=postgresql://postgres@127.0.0.1:5433/uta/uta_20241220 HGVS_SEQREPO_DIR=... \ +python paper/scripts/resolve_clinvar_pass.py clinvar_submitted_pairs_dated.GRCh38.tsv \ + --uta --out submitted_full_uta.csv ``` -`clinvar_submitted_residual.csv` is derived from the **time-bucketed** (headline) cdot -pass rows with `fixed_bucket != correct` (61 of 3,000), the error subtypes recovered by -re-resolving the residual strings and catching the exception class: - -* `coordinate_drift` (29): resolves through the cited historical version to a coordinate - that differs from ClinVar's current interpretation (`incorrect` bucket). -* `reference_mismatch` (18): the cited reference base does not exist on the cited - version (`HGVSInvalidVariantError`). -* `version_refused` (5): cited version absent from the data; the adjacent-version - fallback declined to substitute because coordinate-safety could not be verified - (`no_data` where the accession holds other versions, via `get_tx_versions`). -* `grammar_unsupported` (5): repeat `ref[N]` / allele `[..]` notation the biocommons - grammar rejects (`HGVSParseError`). -* `position_out_of_bounds` (4): the cited position does not exist on the cited version - (`HGVSInvalidIntervalError`). -* `unknown_accession` (0): no version of the accession in the data. +Full head-to-head over 3,198,528 pairs: cdot 99.0% vs UTA 81.9% matched; 548,524 resolve +through cdot alone, 1,226 through UTA alone. Per-era: 2008-2015 cdot 97.3% / UTA 76.2%; +2016-2020 99.3% / 88.8%; 2021-2026 99.0% / 81.0%. + +`clinvar_submitted_residual.csv` is derived from the full cdot pass rows with +`fixed_bucket != correct` (32,462 of 3,198,528); the `no_data` bucket is split into +`version_refused` (accession holds other versions, via `get_tx_versions`) vs +`unknown_accession`, `incorrect` maps to `coordinate_drift`, and the `error` bucket is +re-resolved to recover the exception class (`HGVSInvalidVariantError` reference_mismatch, +`HGVSInvalidIntervalError` position_out_of_bounds, `HGVSParseError` grammar_unsupported). + +## submitter_attribution.csv (frozen) + +R2 submitter attribution of absent-from-cdot transcript versions (why version +substitution is needed at all). `paper/scripts/submitter_attribution.py` streams the +VCV XML, attributes each SCV's first RefSeq transcript citation to its submitting +laboratory, and flags versions cdot does not hold (same test as +`compute_submitted_version_age.py`). Measured 2026-08-19 (ClinVarVCVRelease_2026-06, +cdot 0.2.34 RefSeq GRCh38): + +```bash +python paper/scripts/submitter_attribution.py --xml ClinVarVCVRelease_2026-06.xml.gz \ + --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz --out submitter_attribution.csv +``` + +18,077 of 4,027,987 versioned RefSeq citations (0.449%) cite a version absent from cdot. +Of 2,894 submitters only 85 ever do; a single laboratory contributes 94.9% of absent +citations, and 83.2% of the 286 distinct absent versions are cited by one submitter +alone (the self-alignment signature). Only the aggregate CSV is committed; the +per-submitter table (`--top-out`, which names laboratories) is NOT checked in. + +## benchmark_fullscale.csv (frozen) + +R3 full-scale throughput of a single local-JSON process over the whole current ClinVar +(g.HGVS, c.HGVS) pair set (`clinvar_pairs.GRCh38.tsv`, cdot 0.2.34 RefSeq GRCh38, local +SeqRepo), measured 2026-08-19 with the hot-cache protocol that replaces the old +"cache conditions differed" note: a first pass warms the OS sequence cache and is +discarded, the next is timed. + +```bash +HGVS_SEQREPO_DIR=... SEQREPO_FD_CACHE_MAXSIZE=128 \ +python paper/scripts/resolve_clinvar_pass.py clinvar_pairs.GRCh38.tsv \ + --json cdot-0.2.34.refseq.GRCh38.json.gz --out /dev/null # x2 (discard 1st) +``` + +`n_pairs` 4,423,358; `resolved_pct` 99.4; `cold_tps` 632 HGVS/s (first pass, wall +7003 s); `hot_tps` 640 HGVS/s (timed pass, wall 6916 s), ~1% apart, so full-scale +local-JSON throughput is not sequence-cache-limited. The controlled backend comparison +(local JSON vs prefetched REST) stays in Table 1 / `benchmark.csv`. ## injection_benchmark.csv, residual_taxonomy.csv (table CSVs, not facts) diff --git a/paper/empirical_results/benchmark_fullscale.csv b/paper/empirical_results/benchmark_fullscale.csv new file mode 100644 index 0000000..85aa70a --- /dev/null +++ b/paper/empirical_results/benchmark_fullscale.csv @@ -0,0 +1,6 @@ +field,value +n_pairs,4423358 +resolved_pct,99.4 +cold_tps,632 +hot_tps,640 +wall_min,115 diff --git a/paper/empirical_results/clinvar_submitted.csv b/paper/empirical_results/clinvar_submitted.csv index 20b8635..bdba79f 100644 --- a/paper/empirical_results/clinvar_submitted.csv +++ b/paper/empirical_results/clinvar_submitted.csv @@ -1,2 +1,2 @@ -n_scv_tx_strings,n_unique_pairs,ensembl_pct,version_not_current_pct,scv_weighted_not_current_pct,base_retired_pct,not_current_in_cdot_pct,absent_cdot_pct,n_sample,sample_seed,cdot_resolved_pct,cdot_matched_pct,cdot_no_data,cdot_incorrect,cdot_error,uta_resolved_pct,uta_matched_pct,uta_no_data_pct,cdot_only,cdot_only_pct,uta_only,rescued_by_fix,regressions,after_fix_matched_pct,residual_n,residual_pct,rnd_cdot_resolved_pct,rnd_cdot_matched_pct,rnd_after_fix_matched_pct,rnd_uta_resolved_pct,rnd_uta_matched_pct,rnd_uta_no_data_pct,rnd_cdot_only,rnd_cdot_only_pct,rnd_uta_only,rnd_residual_n,rnd_residual_pct -4027987,3198528,0.0,75.1,69.8,0.7,99.3,0.6,3000,42,98.9,97.9,5,29,27,80.0,79.9,19.5,544,18.1,1,1,0,98.0,61,2.0,99.1,99.0,99.0,82.2,82.2,17.7,504,16.8,1,31,1.0 +n_scv_tx_strings,n_unique_pairs,ensembl_pct,version_not_current_pct,scv_weighted_not_current_pct,base_retired_pct,not_current_in_cdot_pct,absent_cdot_pct,full_n,cdot_resolved_pct,cdot_no_data,cdot_incorrect,cdot_error,uta_resolved_pct,uta_no_data_pct,cdot_only,cdot_only_pct,uta_only,rescued_by_fix,regressions,residual_n,residual_pct,era_old_n,era_old_cdot_pct,era_old_uta_pct,era_mid_n,era_mid_cdot_pct,era_mid_uta_pct,era_recent_n,era_recent_cdot_pct,era_recent_uta_pct +4027987,3198528,0.0,75.1,69.8,0.7,99.3,0.6,3198528,99.0,19124,6338,7000,81.9,18.0,548524,17.2,1226,1428,0,32462,1.0,38091,97.3,76.2,383069,99.3,88.8,2777368,99.0,81.0 diff --git a/paper/empirical_results/clinvar_submitted_residual.csv b/paper/empirical_results/clinvar_submitted_residual.csv index 4e3ec25..d962ab9 100644 --- a/paper/empirical_results/clinvar_submitted_residual.csv +++ b/paper/empirical_results/clinvar_submitted_residual.csv @@ -1,2 +1,2 @@ -n_residual,version_refused,unknown_accession,coordinate_drift,position_out_of_bounds,reference_mismatch,grammar_unsupported -61,5,0,29,4,18,5 +n_residual,version_refused,unknown_accession,coordinate_drift,position_out_of_bounds,reference_mismatch,grammar_unsupported +32462,18925,199,6338,3137,2618,1241 diff --git a/paper/empirical_results/submitter_attribution.csv b/paper/empirical_results/submitter_attribution.csv new file mode 100644 index 0000000..2b71cf8 --- /dev/null +++ b/paper/empirical_results/submitter_attribution.csv @@ -0,0 +1,2 @@ +n_versioned_citations,n_absent_citations,absent_pct,n_submitters_total,n_submitters_with_absent,top1_absent_share_pct,top3_absent_share_pct,top5_absent_share_pct,top10_absent_share_pct,n_absent_versions,n_absent_versions_single_submitter,single_submitter_version_pct,top_submitter_absent_n,top_submitter_absent_rate_pct +4027987,18077,0.449,2894,85,94.9,98.0,98.6,99.1,286,238,83.2,17154,1.0 diff --git a/paper/methods.md b/paper/methods.md index 2c3a7be..2c1317a 100644 --- a/paper/methods.md +++ b/paper/methods.md @@ -168,11 +168,10 @@ The submitted-string corpus (Results R2) is built by VCV XML release (ClinVarVCVRelease_2026-06): each SCV's first transcript c./n. expression is joined verbatim to the variant's VCF coordinate via its AlleleID and collapsed to unique (AlleleID, string) pairs, each tagged with the earliest SCV -submission date (Supplementary Methods). Two 3,000-pair samples are drawn with a fixed -seed (42) for the cdot-versus-UTA comparison: a whole-file random draw, which reflects -the live file and so is recency-biased as ClinVar grows, and a time-bucketed draw that -allocates evenly across submission-year eras so historical submissions are represented -fairly; a 500-pair random sample is committed for reproduction. Scoring uses the VCF +submission date (Supplementary Methods). The cdot-versus-UTA comparison is run over the +whole corpus and broken down by submission-year era, so recency skew is reported rather +than hidden; a 500-pair random sample (seed 42) and the era-tagged corpus support +reproduction. Scoring uses the VCF coordinate rather than the g.HGVS string, since a submitted string may legitimately spell an indel differently from ClinVar's normalised form. Transcript version age is computed by `compute_submitted_version_age.py` against diff --git a/paper/results.md b/paper/results.md index d5e5422..cfa5371 100644 --- a/paper/results.md +++ b/paper/results.md @@ -86,42 +86,47 @@ transcript no longer annotated at any version. cdot's merged release history hol only {{ clinvar_submitted.absent_cdot_pct | dp(1) }}% of cited versions are absent from its GRCh38 data. -Because ClinVar grows over time, a flat random draw over-represents recent submissions, -so we scored two {{ clinvar_submitted.n_sample | commas }}-pair samples (seed 42): one -balanced evenly across submission-year eras (2008-2026), which represents the historical -record fairly, and one drawn at random, which reflects the current file's recency skew -(Methods). - -On the fair era-balanced sample, cdot resolved -{{ clinvar_submitted.cdot_resolved_pct | dp(1) }}% and reproduced ClinVar's VCF -coordinate for {{ clinvar_submitted.after_fix_matched_pct | dp(1) }}%, versus +Because cdot ingests every RefSeq and Ensembl annotation release, a version absent from +its data should be one that never appeared in an official release, and that is what the +absent citations are: attributing each absent-version citation to its submitting +laboratory, of {{ submitter_attribution.n_submitters_total | commas }} submitters only +{{ submitter_attribution.n_submitters_with_absent | int }} ever cite a version cdot +lacks, a single large laboratory accounts for +{{ submitter_attribution.top1_absent_share_pct | dp(1) }}% of all absent-version +citations, and {{ submitter_attribution.single_submitter_version_pct | dp(1) }}% of the +{{ submitter_attribution.n_absent_versions | int }} distinct absent versions are cited by +one submitter alone. These are the signature of transcripts a laboratory aligned to the +genome itself rather than drawing from a published release, so they are resolved by +version substitution (R5), not by deeper ingest. + +cdot resolved the full set of {{ clinvar_submitted.full_n | commas }} submitted pairs, +reproducing ClinVar's VCF coordinate for +{{ clinvar_submitted.cdot_resolved_pct | dp(1) }}%, versus {{ clinvar_submitted.uta_resolved_pct | dp(1) }}% for the same locally loaded UTA. The RefSeq gap invisible at the current-version ceiling (both backends {{ clinvar.cdot_refseq_pct | dp(1) }}% above) opens to {{ clinvar_submitted.cdot_only_pct | dp(1) }} points: -{{ clinvar_submitted.cdot_only | commas }} of the -{{ clinvar_submitted.n_sample | commas }} pairs resolve through cdot alone (one through -UTA alone), because UTA holds no GRCh38 alignment for -{{ clinvar_submitted.uta_no_data_pct | dp(1) }}% of the cited versions. The random draw -is easier on both backends ({{ clinvar_submitted.rnd_after_fix_matched_pct | dp(1) }}% -cdot versus {{ clinvar_submitted.rnd_uta_resolved_pct | dp(1) }}% UTA, a -{{ clinvar_submitted.rnd_cdot_only_pct | dp(1) }}-point gap): it is dominated by recent -submissions, and the fair sample is harder precisely because older submissions cite the -superseded versions this corpus exists to exercise. +{{ clinvar_submitted.cdot_only | commas }} pairs resolve through cdot alone and +{{ clinvar_submitted.uta_only | commas }} through UTA alone, because UTA holds no GRCh38 +alignment for {{ clinvar_submitted.uta_no_data_pct | dp(1) }}% of the cited versions. +Splitting by submission era locates the gap: cdot stays near +{{ clinvar_submitted.era_recent_cdot_pct | dp(1) }}% across all years, while UTA falls +from {{ clinvar_submitted.era_recent_uta_pct | dp(1) }}% on recent (2021-2026) +submissions to {{ clinvar_submitted.era_old_uta_pct | dp(1) }}% on the oldest +(2008-2015), which cite the most superseded versions. The submitted strings are largely well-formed, so cleaning has little to rescue -(`fix_hgvs()` repaired {{ clinvar_submitted.rescued_by_fix | int }} string, with +(`fix_hgvs()` repaired {{ clinvar_submitted.rescued_by_fix | commas }} strings, with {{ clinvar_submitted.regressions | int }} regressions). The residual -{{ clinvar_submitted.residual_pct | dp(1) }}% on the fair sample after cleaning and -version substitution ({{ clinvar_submitted.residual_n | int }} of -{{ clinvar_submitted.n_sample | commas }}) is dominated by historical-version effects, -not formatting: {{ clinvar_submitted_residual.coordinate_drift | int }} resolve through -the cited version to a coordinate that differs from ClinVar's current interpretation, -{{ clinvar_submitted_residual.reference_mismatch | int }} cite a reference base that -does not exist on the cited version, and -{{ clinvar_submitted_residual.version_refused | int }} cite a version absent from the -data where substitution declines to act because coordinate safety cannot be verified; -the full breakdown is in Supplementary Table S8. +{{ clinvar_submitted.residual_pct | dp(1) }}% after cleaning and version substitution +({{ clinvar_submitted.residual_n | commas }} of {{ clinvar_submitted.full_n | commas }}) +is dominated not by formatting but by historical-version effects: +{{ clinvar_submitted_residual.version_refused | commas }} cite a version absent from the +data where substitution declines to act because coordinate safety cannot be verified, +and {{ clinvar_submitted_residual.coordinate_drift | commas }} resolve through the cited +version to a coordinate that differs from ClinVar's current interpretation; the remainder +cite a position or base that does not exist on the cited version, or notation the +biocommons grammar rejects (Supplementary Table S8). **[private data].** The same gap holds on the historical clinical data that motivated cdot: the complete set of {{ historical.n_lines | commas }} unique HGVS descriptions imported @@ -175,14 +180,17 @@ request per transcript, no prefetch) is more than two orders of magnitude faster public UTA server's {{ benchmark.uta_remote_tps | dp(2) }} HGVS/s, which pays wide-area round trips to a shared database on every lookup. -At scale, a single local-JSON process resolved the entire set of 3,660,452 unique -ClinVar (g.HGVS, c.HGVS) pairs in ~92 minutes (665 HGVS/s; 99.3% produced a genomic -coordinate, 98.8% matched the ClinVar genomic HGVS exactly), and the REST provider, -after one batch cache warm (21,277 distinct transcripts in ~6 s), in ~83 minutes -(731 HGVS/s). The two runs' sequence-layer cache conditions differed, so their small -difference is not evidence of REST outrunning local JSON; the controlled comparison is -Table 1. The same exhaustive pass against the public remote UTA database extrapolates -to close to a year. +At scale, a single local-JSON process resolved the full set of +{{ benchmark_fullscale.n_pairs | commas }} unique ClinVar (g.HGVS, c.HGVS) pairs in +~{{ benchmark_fullscale.wall_min | int }} minutes, +{{ benchmark_fullscale.resolved_pct | dp(1) }}% producing a genomic coordinate. Timing a +hot-cache repeat (a first pass discarded to warm the sequence cache, the next timed) gave +{{ benchmark_fullscale.hot_tps | int }} HGVS/s, within 1% of the cold pass +({{ benchmark_fullscale.cold_tps | int }} HGVS/s), so full-scale local-JSON throughput is +bounded by the shared engine and sequence layer, not by the state of the sequence cache. +The controlled like-for-like backend comparison, with the sequence layer held identical, +is Table 1, where local JSON and prefetched REST fall within 1% of each other. The same +exhaustive pass against the public remote UTA database extrapolates to close to a year. ## R4: String cleaning recovers malformed real-world HGVS diff --git a/paper/scripts/submitter_attribution.py b/paper/scripts/submitter_attribution.py new file mode 100644 index 0000000..d26e5b1 --- /dev/null +++ b/paper/scripts/submitter_attribution.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +"""Which submitters cite transcript versions cdot does not hold? (paper R2 / Discussion) + +cdot ingests every RefSeq/Ensembl annotation release off the FTP sites, so in +principle it should hold every *published* transcript version. Yet a small fraction +of submitted ClinVar strings cite a version absent from cdot's data +(``absent_cdot_grch38`` in compute_submitted_version_age.py). Hypothesis (Dave): those +versions never appeared in an official annotation release, because a small number of +laboratories align transcripts to the genome themselves and cite the resulting +idiosyncratic version. This script tests that by attributing each absent-version +citation to its submitting laboratory. + +Method: stream the ClinVar VCV XML; for each ClinicalAssertion (SCV) take its submitter +(``ClinVarAccession@SubmitterName``, else ``ClinVarSubmissionID@submitter``) and its +first RefSeq transcript c./n. expression. A cited ``base.version`` is *absent* if that +exact version is not among the versions cdot holds for the accession (same test as +compute_submitted_version_age.py, reusing ``scan_transcript_urls``). We then accumulate, +per submitter, total versioned RefSeq citations and absent-version citations, and per +absent version how many distinct submitters cite it. + +Concentration is the question: if absent citations come from a few labs, and absent +versions are each cited by a single submitter, that is the self-alignment signature; if +absent versions are widely cited, they are a genuine ingest gap. + +Usage: + python paper/scripts/submitter_attribution.py \ + --xml ClinVarVCVRelease_2026-06.xml.gz \ + --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz \ + --out output/facts/submitter_attribution.csv [--limit N] [--top-out top.csv] +""" +import argparse +import csv +import gzip +import re +import sys +from collections import defaultdict +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from compute_submitted_version_age import scan_transcript_urls # noqa: E402 + +_TX = re.compile(r"^((?:NM|NR|XM|XR)_\d+)\.(\d+):[cn]\.") + + +def iter_scv_submitter_hgvs(xml_path, limit=0): + """Yield (submitter, base, version) for each SCV's first RefSeq transcript c./n. + expression, streaming the VCV XML one VariationArchive at a time.""" + try: + from lxml import etree + except ImportError: # pragma: no cover + import xml.etree.ElementTree as etree + opener = gzip.open if str(xml_path).endswith(".gz") else open + n_va = 0 + with opener(xml_path, "rb") as fh: + for _ev, elem in etree.iterparse(fh, events=("end",), tag="VariationArchive"): + n_va += 1 + record = elem.find("ClassifiedRecord") + if record is None: + record = elem.find("InterpretedRecord") + if record is not None: + for ca in record.iter("ClinicalAssertion"): + acc = ca.find("ClinVarAccession") + sub = (acc.get("SubmitterName") if acc is not None else None) + if not sub: + sid = ca.find("ClinVarSubmissionID") + sub = sid.get("submitter") if sid is not None else None + sub = (sub or "").strip() or "(unknown)" + ca_sa = ca.find("SimpleAllele") + if ca_sa is None: + continue + for attr in ca_sa.iter("Attribute"): + if attr.get("Type") == "HGVS" and attr.text: + m = _TX.match(attr.text.strip()) + if m: + yield sub, m.group(1), int(m.group(2)) + break + elem.clear() + while elem.getprevious() is not None: + del elem.getparent()[0] + if limit and n_va >= limit: + return + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--xml", required=True, help="ClinVar VCV XML release (.xml.gz)") + ap.add_argument("--refseq-grch38", required=True, help="cdot RefSeq GRCh38 JSON.gz") + ap.add_argument("--limit", type=int, default=0, help="stop after N VariationArchives (smoke test)") + ap.add_argument("--out", default="output/facts/submitter_attribution.csv") + ap.add_argument("--top-out", help="optional: per-submitter table CSV (all submitters with absent citations)") + args = ap.parse_args() + + print(f"scanning {args.refseq_grch38} ...", file=sys.stderr) + tx_url = scan_transcript_urls(args.refseq_grch38, keys_only=True) + versions_of = defaultdict(set) + for acc in tx_url: + base, ver = acc.rsplit(".", 1) + versions_of[base].add(int(ver)) + print(f" {len(tx_url):,} versions across {len(versions_of):,} accessions", file=sys.stderr) + + total = defaultdict(int) # submitter -> versioned RefSeq citations + absent = defaultdict(int) # submitter -> absent-version citations + ver_submitters = defaultdict(set) # base.ver -> {submitters citing it} + n_scv = 0 + for sub, base, ver in iter_scv_submitter_hgvs(args.xml, args.limit): + n_scv += 1 + total[sub] += 1 + if ver not in versions_of.get(base, ()): + absent[sub] += 1 + ver_submitters[f"{base}.{ver}"].add(sub) + + n_absent = sum(absent.values()) + subs_with_absent = sorted(absent.items(), key=lambda kv: -kv[1]) + + def share(k): + return round(100 * sum(c for _, c in subs_with_absent[:k]) / n_absent, 1) if n_absent else 0.0 + + # self-alignment signature: absent versions cited by exactly one submitter + n_single_versions = sum(1 for subs in ver_submitters.values() if len(subs) == 1) + + row = { + "n_versioned_citations": n_scv, + "n_absent_citations": n_absent, + "absent_pct": round(100 * n_absent / n_scv, 3) if n_scv else 0.0, + "n_submitters_total": len(total), + "n_submitters_with_absent": len(subs_with_absent), + "top1_absent_share_pct": share(1), + "top3_absent_share_pct": share(3), + "top5_absent_share_pct": share(5), + "top10_absent_share_pct": share(10), + "n_absent_versions": len(ver_submitters), + "n_absent_versions_single_submitter": n_single_versions, + "single_submitter_version_pct": round(100 * n_single_versions / len(ver_submitters), 1) if ver_submitters else 0.0, + "top_submitter_absent_n": subs_with_absent[0][1] if subs_with_absent else 0, + "top_submitter_absent_rate_pct": round(100 * subs_with_absent[0][1] / total[subs_with_absent[0][0]], 1) if subs_with_absent else 0.0, + } + out = Path(args.out) + out.parent.mkdir(parents=True, exist_ok=True) + with open(out, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=list(row)) + w.writeheader() + w.writerow(row) + for k, v in row.items(): + print(f" {k:<36} {v}") + print(f"Written: {out}") + + if args.top_out: + with open(args.top_out, "w", newline="") as fh: + w = csv.writer(fh) + w.writerow(["submitter", "absent_citations", "total_citations", "absent_rate_pct"]) + for sub, a in subs_with_absent: + w.writerow([sub, a, total[sub], round(100 * a / total[sub], 1)]) + print(f"Written per-submitter table: {args.top_out}", file=sys.stderr) + + +if __name__ == "__main__": + main() diff --git a/paper/supplementary.md b/paper/supplementary.md index e7d0672..ee8ecaf 100644 --- a/paper/supplementary.md +++ b/paper/supplementary.md @@ -12,8 +12,8 @@ what defines "the" submitted string when an SCV carries several HGVS expressions Submitted HGVS is not deduplicated across SCVs in ClinVar, so the builder collapses SCV rows to unique (AlleleID, submitted string) pairs, retaining the collapsed submission count; distinct strings for the same variant (for example, two laboratories citing -different transcript versions) are all kept. The 3,000-pair benchmark sample is sized -to match the seeded current-version sample used for the cdot-versus-UTA comparison. +different transcript versions) are all kept. The cdot-versus-UTA comparison (Results R2) +is run over the whole corpus; the committed 500-pair sample is for quick reproduction. ### LOVD head-to-head scoring detail @@ -126,11 +126,11 @@ refused by default. ### Table S4: ClinVar benchmark details Full-scale resolution of every RefSeq and Ensembl c.HGVS in ClinVar through cdot alone -(GRCh38, cdot 0.2.33). Unlike the cdot-vs-UTA comparison in Results R2 (gated to a sample -by UTA's throughput), every variant is summarised here because only cdot is in the loop. The -projection is scored as a VCF coordinate (CHROM/POS/REF/ALT) against ClinVar's own VCF, -not as a g.HGVS string, so equivalent representations and ClinVar's tandem-repeat / identity -notations are not miscounted. +(GRCh38, cdot 0.2.33): a cdot-only accuracy check on the current-version `Name` corpus, +complementary to the cdot-vs-UTA head-to-head on the historical submitted-string corpus +(Results R2). The projection is scored as a VCF coordinate (CHROM/POS/REF/ALT) against +ClinVar's own VCF, not as a g.HGVS string, so equivalent representations and ClinVar's +tandem-repeat / identity notations are not miscounted. **Caveat.** ClinVar submissions are dominated by a handful of large (largely US) clinical laboratories citing mostly current RefSeq versions, so this is a clean, public, @@ -330,17 +330,17 @@ supplied) restores a fully dropped RefSeq prefix (`000059.4:c.68del` → applying the fix only when exactly one candidate accession exists in the loaded data (Methods). -### Table S8: Residual outcomes on the submitted-string sample +### Table S8: Residual outcomes on the submitted-string corpus -The {{ clinvar_submitted.residual_n | int }} of -{{ clinvar_submitted.n_sample | commas }} sampled submitted-string pairs (Results R2) +The {{ clinvar_submitted.residual_n | commas }} of +{{ clinvar_submitted.full_n | commas }} submitted-string pairs (Results R2) that remain unresolved after cleaning and version fallback, by cause. | Cause | n | |---|---| -| Cited version absent from the data; the fallback declines to substitute because coordinate safety cannot be verified (no false rescues, by design) | {{ clinvar_submitted_residual.version_refused | int }} | -| Resolves through the cited historical version to a coordinate that differs from ClinVar's current interpretation | {{ clinvar_submitted_residual.coordinate_drift | int }} | -| Cited position does not exist on the cited version | {{ clinvar_submitted_residual.position_out_of_bounds | int }} | -| Cited reference base does not exist on the cited version | {{ clinvar_submitted_residual.reference_mismatch | int }} | -| Repeat or allele notation the biocommons grammar rejects | {{ clinvar_submitted_residual.grammar_unsupported | int }} | -| Accession absent entirely | {{ clinvar_submitted_residual.unknown_accession | int }} | +| Cited version absent from the data; the fallback declines to substitute because coordinate safety cannot be verified (no false rescues, by design) | {{ clinvar_submitted_residual.version_refused | commas }} | +| Resolves through the cited historical version to a coordinate that differs from ClinVar's current interpretation | {{ clinvar_submitted_residual.coordinate_drift | commas }} | +| Cited position does not exist on the cited version | {{ clinvar_submitted_residual.position_out_of_bounds | commas }} | +| Cited reference base does not exist on the cited version | {{ clinvar_submitted_residual.reference_mismatch | commas }} | +| Repeat or allele notation the biocommons grammar rejects | {{ clinvar_submitted_residual.grammar_unsupported | commas }} | +| Accession absent entirely | {{ clinvar_submitted_residual.unknown_accession | commas }} |