Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 36 additions & 0 deletions claude/plans/new_analyses_plan.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion paper/Snakefile
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
104 changes: 65 additions & 39 deletions paper/empirical_results/PROVENANCE.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)

Expand Down
6 changes: 6 additions & 0 deletions paper/empirical_results/benchmark_fullscale.csv
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
field,value
n_pairs,4423358
resolved_pct,99.4
cold_tps,632
hot_tps,640
wall_min,115
4 changes: 2 additions & 2 deletions paper/empirical_results/clinvar_submitted.csv
Original file line number Diff line number Diff line change
@@ -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
4 changes: 2 additions & 2 deletions paper/empirical_results/clinvar_submitted_residual.csv
Original file line number Diff line number Diff line change
@@ -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
2 changes: 2 additions & 0 deletions paper/empirical_results/submitter_attribution.csv
Original file line number Diff line number Diff line change
@@ -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
9 changes: 4 additions & 5 deletions paper/methods.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading