From 6e2d163b151be6c9bef562aa70b35cf9e80a86ab Mon Sep 17 00:00:00 2001 From: Dave Lawrence Date: Mon, 17 Aug 2026 15:04:51 +0000 Subject: [PATCH 1/3] clean_hgvs repair ops and get_tx_exons caching - #112 #5 --- CHANGELOG.md | 6 + cdot/hgvs/__init__.py | 1 + cdot/hgvs/clean.py | 72 ++++++- cdot/hgvs/dataproviders/json_data_provider.py | 27 ++- cdot/hgvs/gene_hgvs.py | 115 ++++++++++- claude/20260808_paper_feedback_plan.md | 195 ------------------ tests/test_clean_hgvs.py | 43 ++++ tests/test_gene_hgvs.py | 114 +++++++++- 8 files changed, 367 insertions(+), 206 deletions(-) delete mode 100644 claude/20260808_paper_feedback_plan.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ff0bf5..75b2183 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,8 +2,14 @@ ### Added +- #112 - HGVS cleaning: `resolve_missing_accession_prefix(hgvs_string, data_provider, genome_build=None)` restores a RefSeq prefix the user dropped entirely, eg `000059.4:c.68del` → `NM_000059.4:c.68del`. The kind letter picks the candidate prefixes (`c.` → `NM_`/`XM_`, `n.` → `NR_`/`XR_`; a 7-8 digit field is zero-padded to the 9-digit form) and each candidate is checked against the data provider (via `get_tx_versions`), so the fix is applied only when exactly one candidate accession exists there - with zero or several matches the string is unchanged. Reported as a WARNING `HGVSFix` with the new code `ADDED_ACCESSION_PREFIX`, never silent. `fix_hgvs` runs it automatically whenever a `data_provider` is supplied (before gene resolution); without a provider the string is left alone +- #112 - HGVS cleaning: new string-level ops driven by the search-log corpus - a semicolon or underscore in place of the accession:allele colon (`NM_000059.4;c.68del` / `NM_001754.5_c.749G>A` → `…:c.…`, reported under `FIXED_SEPARATOR_TYPO`); a space in place of the reference:kind colon after a bare gene symbol (`BRCA2 c.68del` → `BRCA2:c.68del`, new op `HGVSCleanOp.FIX_SPACE_BEFORE_KIND`, which also repairs symbols that resemble an accession, eg `NR2E3 c.349G>A`); a dangling dot left by a dropped transcript version (`NM_001754.:c.749G>A` → `NM_001754:c.749G>A`, new op `HGVSCleanOp.DROP_EMPTY_VERSION` / fix code `DROPPED_EMPTY_VERSION`); and leading-junk stripping now also covers junk separated from the accession by whitespace or a stray dash (`": NM_000059.4(BRCA2):c.68del"`) - #51 - RefSeq GRCh38 data now includes NCBI's historical transcript alignments (`RefSeq_historical_alignments`, RS_2024_08 set): alignments of replaced and suppressed NM_/NR_ versions, many of which never appeared in any annotation release cdot ingests. They are merged at the lowest priority above UTA, so a transcript version from any official annotation release is unchanged. Note some historical alignments are partial: a few transcript bases at the ends (eg an unaligned poly-A tail) or occasionally the first bases do not align to the genome, and the stored exon cDNA coordinates reflect that, so positions in an unaligned region cannot be projected. Data only (no client code change), from the next data release +### Changed + +- #5 - Performance: the biocommons HGVS data providers now keep a small per-instance LRU cache (10,000 entries) of `get_tx_exons()` results, so repeated lookups of the same transcript/contig no longer rebuild the per-exon dict list (the per-variant hot path when resolving a stream of variants). Applies to `JSONDataProvider` and `RESTDataProvider`; results are cached only after a successful build, so lazily fetched providers still pick up transcripts that appear later (eg after a REST `prefetch()`) + ### Fixed - #95 - Data fix: the genomic `cds_start`/`cds_end` of UTA-sourced minus-strand transcripts were both 1 too high (a 0-based/1-based conversion error in the UTA import; plus-strand transcripts were unaffected). These fields are only used by the PyHGVS integration, so biocommons HGVS resolution was never affected. The UTA import also no longer requires the SACGF PyHGVS fork: when the pipeline ran with standard pyhgvs installed instead, every coding UTA transcript failed conversion and was silently skipped, which is why release 0.2.33 has only 3,954 UTA-sourced GRCh37 records (all non-coding) against 46,659 in 0.2.28. Coding UTA transcripts are restored from the next data release, and the pipeline's UTA fetch now fails loudly instead of keeping a truncated download. Data only (no client code change) diff --git a/cdot/hgvs/__init__.py b/cdot/hgvs/__init__.py index 2ff9c30..850d700 100644 --- a/cdot/hgvs/__init__.py +++ b/cdot/hgvs/__init__.py @@ -22,6 +22,7 @@ fix_hgvs, rank_transcripts_for_gene, resolve_gene_hgvs, + resolve_missing_accession_prefix, resolve_transcript_version, UnsafeVersionPolicy, ) diff --git a/cdot/hgvs/clean.py b/cdot/hgvs/clean.py index 2e29502..f5e1085 100644 --- a/cdot/hgvs/clean.py +++ b/cdot/hgvs/clean.py @@ -54,6 +54,7 @@ class HGVSFixCode(Enum): DROPPED_GENOMIC_REF_IN_PARENS = "dropped_genomic_ref_in_parens" ADDED_TRANSCRIPT_UNDERSCORE = "added_transcript_underscore" FIXED_MULTIPLE_UNDERSCORE = "fixed_multiple_underscore" + DROPPED_EMPTY_VERSION = "dropped_empty_version" UPPERCASED_BASES = "uppercased_bases" RECONSTRUCTED_STRUCTURE = "reconstructed_structure" ADDED_MISSING_KIND = "added_missing_kind" @@ -74,6 +75,9 @@ class HGVSFixCode(Enum): # Gene → transcript resolution (WARNING on success, ERROR on failure) RESOLVED_GENE_TO_TRANSCRIPT = "resolved_gene_to_transcript" NO_TRANSCRIPT_FOR_GENE = "no_transcript_for_gene" + # Accession prefix restoration for a bare-number accession, verified against + # a data provider (see gene_hgvs.resolve_missing_accession_prefix; WARNING) + ADDED_ACCESSION_PREFIX = "added_accession_prefix" # Validation errors (ERROR) NO_COLON = "no_colon" MISSING_REFERENCE_SEQUENCE = "missing_reference_sequence" @@ -92,6 +96,7 @@ class HGVSCleanOp(Enum): """ STRIP_PROTEIN_SUFFIX = "strip_protein_suffix" STRIP_LEADING_JUNK = "strip_leading_junk" + FIX_SPACE_BEFORE_KIND = "fix_space_before_kind" STRIP_WHITESPACE = "strip_whitespace" STRIP_SURROUNDING_PUNCTUATION = "strip_surrounding_punctuation" FIX_MULTIPLE_COLON = "fix_multiple_colon" @@ -101,6 +106,7 @@ class HGVSCleanOp(Enum): FIX_PREFIX_COLON = "fix_prefix_colon" FIX_MULTIPLE_KIND = "fix_multiple_kind" FIX_SEPARATOR_TYPO = "fix_separator_typo" + DROP_EMPTY_VERSION = "drop_empty_version" ADD_N_PREFIX = "add_n_prefix" LOWERCASE_MUTATION_TYPE = "lowercase_mutation_type" DROP_DEL_DUP_COUNT = "drop_del_dup_count" @@ -244,9 +250,10 @@ def _looks_like_hgvs_prefix(s: str) -> bool: # Transcript accession prefixes (used to locate the HGVS core within a string) _TX_PREFIX = r"(?:NM_|NR_|XM_|XR_|NC_|NG_|ENST|LRG_)" -# Leading junk before a transcript accession, e.g. "#", ":", a literal "\t", -# or a genome-build prefix like "GRCh38.p2 " / "GRCh37.p12#" -_LEADING_JUNK = re.compile(r"^(?:GRCh\d+(?:\.p\d+)?[#\s]*|\\t|[#:])+", re.IGNORECASE) +# Leading junk before a transcript accession, e.g. "#", ":", "-", a literal +# "\t", whitespace (": NM_..."), or a genome-build prefix like "GRCh38.p2 " / +# "GRCh37.p12#". Only stripped when a transcript accession actually follows. +_LEADING_JUNK = re.compile(r"^(?:GRCh\d+(?:\.p\d+)?[#\s]*|\\t|[#:\-\s])+", re.IGNORECASE) # Gene symbol wedged between extra colons, e.g. "tx:(GENE):c." / "tx:(GENE)c." / # "tx:GENE:c." — should be "tx(GENE):c." @@ -270,6 +277,26 @@ def _looks_like_hgvs_prefix(s: str) -> bool: # Period used in place of the substitution '>', e.g. "1030C.T" -> "1030C>T" _SUB_PERIOD = re.compile(r"([ACGT])\.([ACGT])$", re.IGNORECASE) +# Semicolon or underscore used in place of the accession:allele colon, e.g. +# "NM_000059.4;c.68del" or "NM_000059.4_c.68del". Anchored to the start so a +# ';' inside a multi-variant allele is never touched. LRG is excluded (its +# accessions legitimately contain underscores). +_ACCESSION_KIND_SEPARATOR = re.compile( + r"^((?:NM_|NR_|XM_|XR_|NC_|NG_|ENST)\d+(?:\.\d+)?)[;_](?=[cgnmp]\.)", re.IGNORECASE) + +# Whitespace used in place of the reference:kind colon after a bare gene symbol, +# e.g. "BRCA2 c.68del" -> "BRCA2:c.68del". Restricted to a plain symbol (no +# '_'/'.'/':') so accessions are left to the whitespace/reconstruction steps, +# and runs before whitespace stripping (which would otherwise glue the symbol to +# the kind letter, mangling symbols that resemble an accession, eg NR2E3). +_SPACE_BEFORE_KIND = re.compile(r"^([A-Za-z][A-Za-z0-9]{1,9})\s+(?=[cgnmp]\.)", re.IGNORECASE) + +# An empty version field: accession, a dot, then no digits before the colon, +# e.g. "NM_001754.:c.749G>A". The version digits were dropped, so drop the dot +# and leave a valid unversioned accession. +_EMPTY_VERSION = re.compile( + r"^((?:NM_|NR_|XM_|XR_|NC_|NG_|ENST)\d+)\.(?=:)", re.IGNORECASE) + # del/dup with an explicit range followed by a redundant base count, e.g. # "c.1315_1337dup23" -> "c.1315_1337dup" (the range already gives the length; # normalise to plain del/dup). The range guard avoids changing the meaning of a @@ -723,12 +750,45 @@ def _op_fix_gene_wrapper(s: str) -> tuple[str, list[HGVSFix]]: return s, [] +def _op_fix_space_before_kind(s: str) -> tuple[str, list[HGVSFix]]: + # Whitespace in place of the reference:kind colon after a bare gene symbol, + # e.g. "BRCA2 c.68del" -> "BRCA2:c.68del". Must run before whitespace + # stripping, which would glue the symbol to the kind letter and leave + # accession-lookalike symbols (eg NR2E3) unrepairable. + s2 = _SPACE_BEFORE_KIND.sub(r"\1:", s) + if s2 != s: + return s2, [HGVSFix( + severity=HGVSFixSeverity.WARNING, + code=HGVSFixCode.FIXED_SEPARATOR_TYPO, + message="Replaced space between reference and kind with ':'", + original=s, fixed=s2, + )] + return s, [] + + +def _op_drop_empty_version(s: str) -> tuple[str, list[HGVSFix]]: + # An accession with a dot but no version digits, e.g. "NM_001754.:c.749G>A" + # -> "NM_001754:c.749G>A". The version was dropped; an unversioned accession + # is valid HGVS, so drop the dangling dot rather than invent a version. + s2 = _EMPTY_VERSION.sub(r"\1", s) + if s2 != s: + return s2, [HGVSFix( + severity=HGVSFixSeverity.WARNING, + code=HGVSFixCode.DROPPED_EMPTY_VERSION, + message="Dropped dangling dot of an empty transcript version", + original=s, fixed=s2, + )] + return s, [] + + def _op_fix_separator_typo(s: str) -> tuple[str, list[HGVSFix]]: - # Comma ("c,1811") or colon ("c:1811") in place of the kind dot, and period - # in place of the substitution '>' ("1030C.T" -> "1030C>T"). + # Comma ("c,1811") or colon ("c:1811") in place of the kind dot, period + # in place of the substitution '>' ("1030C.T" -> "1030C>T"), and semicolon + # or underscore in place of the accession:allele colon ("NM_000059.4;c.68del"). s2 = _KIND_COMMA.sub(r":\1.", s) s2 = _KIND_COLON.sub(r":\1.", s2) s2 = _SUB_PERIOD.sub(r"\1>\2", s2) + s2 = _ACCESSION_KIND_SEPARATOR.sub(r"\1:", s2) if s2 != s: return s2, [HGVSFix( severity=HGVSFixSeverity.WARNING, @@ -777,6 +837,7 @@ def _op_drop_genomic_ref_in_parens(s: str) -> tuple[str, list[HGVSFix]]: _PIPELINE: list[tuple[HGVSCleanOp, "callable"]] = [ (HGVSCleanOp.STRIP_PROTEIN_SUFFIX, _op_strip_protein_suffix), (HGVSCleanOp.STRIP_LEADING_JUNK, _op_strip_leading_junk), + (HGVSCleanOp.FIX_SPACE_BEFORE_KIND, _op_fix_space_before_kind), (HGVSCleanOp.STRIP_WHITESPACE, _op_strip_whitespace), (HGVSCleanOp.STRIP_SURROUNDING_PUNCTUATION, _op_strip_surrounding_punctuation), (HGVSCleanOp.FIX_MULTIPLE_COLON, _op_fix_multiple_colon), @@ -786,6 +847,7 @@ def _op_drop_genomic_ref_in_parens(s: str) -> tuple[str, list[HGVSFix]]: (HGVSCleanOp.FIX_PREFIX_COLON, _op_fix_prefix_colon), (HGVSCleanOp.FIX_MULTIPLE_KIND, _op_fix_multiple_kind), (HGVSCleanOp.FIX_SEPARATOR_TYPO, _op_fix_separator_typo), + (HGVSCleanOp.DROP_EMPTY_VERSION, _op_drop_empty_version), (HGVSCleanOp.ADD_N_PREFIX, _op_add_n_prefix), (HGVSCleanOp.LOWERCASE_MUTATION_TYPE, _op_lowercase_mutation_type), (HGVSCleanOp.DROP_DEL_DUP_COUNT, _op_drop_del_dup_count), diff --git a/cdot/hgvs/dataproviders/json_data_provider.py b/cdot/hgvs/dataproviders/json_data_provider.py index d8cf80e..3cfdc0b 100644 --- a/cdot/hgvs/dataproviders/json_data_provider.py +++ b/cdot/hgvs/dataproviders/json_data_provider.py @@ -2,7 +2,7 @@ import gzip import requests -from collections import defaultdict +from collections import defaultdict, OrderedDict from lazy import lazy from hgvs.dataproviders.interface import Interface from hgvs.dataproviders.seqfetcher import SeqFetcher @@ -31,6 +31,12 @@ class AbstractJSONDataProvider(Interface): # All cdot data is 'splign', it's the method used in NCBI/Ensembl GTFs, and we also only pull out 'splign' from UTA NCBI_ALN_METHOD = "splign" required_version = "1.1" + # get_tx_exons() converts cdot's stored exon arrays into the list of per-exon dicts + # (with CIGARs) that biocommons hgvs expects. That conversion is the per-variant hot + # path, and its input is immutable once loaded/fetched, so each provider instance + # keeps a small LRU of built results (see get_tx_exons). Size is a per-instance + # bound on distinct (tx_ac, alt_ac) entries; each entry is just the exon dict list. + TX_EXONS_CACHE_SIZE = 10_000 def __init__(self, assemblies: List[str] = None, mode=None, cache=None, seqfetcher=None): """ assemblies: defaults to ["GRCh37", "GRCh38"] @@ -54,6 +60,10 @@ def __init__(self, assemblies: List[str] = None, mode=None, cache=None, seqfetch self.assembly_by_contig = {} for assembly_name, contig_map in self.assembly_maps.items(): self.assembly_by_contig.update({contig: assembly_name for contig in contig_map.keys()}) + # Per-instance LRU for get_tx_exons results, keyed (tx_ac, alt_ac). Instance-level + # (not functools.lru_cache on the method) so the cache dies with the instance and + # is never shared across providers holding different data. + self._tx_exons_cache = OrderedDict() @abc.abstractmethod @@ -283,7 +293,17 @@ def _convert_gap_to_cigar(gap): return "".join(cigar_ops) def get_tx_exons(self, tx_ac, alt_ac, alt_aln_method): - self._check_alt_aln_method(alt_aln_method) + self._check_alt_aln_method(alt_aln_method) # raises unless alt_aln_method == splign + # alt_aln_method is validated above (always splign) so it is not part of the key. + # Only successful builds are cached: a None result can become available later on + # providers that fetch lazily (eg RESTDataProvider after prefetch()). + # Callers get the cached list itself and must not mutate it (biocommons hgvs + # only reads it; sorted() takes a copy). + cache_key = (tx_ac, alt_ac) + if (cached := self._tx_exons_cache.get(cache_key)) is not None: + self._tx_exons_cache.move_to_end(cache_key) + return cached + transcript = self._get_transcript(tx_ac) if not transcript: return None @@ -320,6 +340,9 @@ def get_tx_exons(self, tx_ac, alt_ac, alt_aln_method): } tx_exons.append(exon_data) + self._tx_exons_cache[cache_key] = tx_exons + if len(self._tx_exons_cache) > self.TX_EXONS_CACHE_SIZE: + self._tx_exons_cache.popitem(last=False) return tx_exons def get_tx_identity_info(self, tx_ac): diff --git a/cdot/hgvs/gene_hgvs.py b/cdot/hgvs/gene_hgvs.py index 38f5af4..2f18854 100644 --- a/cdot/hgvs/gene_hgvs.py +++ b/cdot/hgvs/gene_hgvs.py @@ -11,6 +11,7 @@ """ import inspect +import re from enum import Enum from typing import Optional @@ -343,6 +344,105 @@ def resolve_gene_hgvs( return resolved, fixes +# --------------------------------------------------------------------------- +# Public: resolve_missing_accession_prefix — restore a dropped RefSeq prefix +# --------------------------------------------------------------------------- + +# A bare-number accession: the user dropped the RefSeq prefix entirely (or left +# only its underscore), e.g. "000059.4:c.68del" / "_001128425.1:c.667A>G". +# RefSeq digit fields are 6 or 9 digits, but 7-8 are accepted as a 9-digit +# accession missing its leading zeros (see _digit_field_candidates). +_BARE_NUMBER_ACCESSION = re.compile(r"^_?(\d{6,9})((?:\.\d+)?):([cn])\.", re.IGNORECASE) + +# The kind letter constrains which RefSeq prefixes could have been dropped: +# a coding "c." implies an mRNA accession, a non-coding "n." an RNA accession. +_PREFIXES_FOR_KIND = { + "c": ("NM_", "XM_"), + "n": ("NR_", "XR_"), +} + + +def _digit_field_candidates(digits: str) -> list[str]: + """Candidate digit fields for a bare-number accession. + + RefSeq accessions have 6 or 9 digits after the prefix. A 6 or 9 digit + input is taken as written; a 7-8 digit input can only be a 9-digit + accession whose leading zeros were also dropped, so it is zero-padded. + """ + if len(digits) in (6, 9): + return [digits] + return [digits.zfill(9)] + + +def resolve_missing_accession_prefix( + hgvs_string: str, + data_provider, + genome_build: Optional[str] = None, +) -> tuple[str, list[HGVSFix]]: + """ + Restore a RefSeq prefix that was dropped from the transcript accession, + e.g. "000059.4:c.68del" → "NM_000059.4:c.68del". + + The kind letter narrows the possibilities ("c." implies NM_/XM_, "n." + implies NR_/XR_) but cannot fully disambiguate, so the candidates are + checked against the data provider (via ``get_tx_versions``) and the fix is + applied only when exactly one candidate accession exists there. With zero + or several matches the string is returned unchanged - this function never + guesses. + + Returns: + (resolved_string, fixes) + fixes contains a WARNING ADDED_ACCESSION_PREFIX HGVSFix when the prefix + was restored, and is empty otherwise (no bare-number accession, no + provider match, ambiguous match, or a provider that cannot enumerate + versions). + + ``genome_build`` is optional and only used to restrict the existence check + to accessions placeable in that build (for providers whose + ``get_tx_versions`` accepts it). + + Example:: + + resolved, fixes = resolve_missing_accession_prefix( + "000059.4:c.68del", data_provider) + # resolved = "NM_000059.4:c.68del" (if only NM_000059 exists) + # fixes[0].code = HGVSFixCode.ADDED_ACCESSION_PREFIX + """ + m = _BARE_NUMBER_ACCESSION.match(hgvs_string) + if m is None: + return hgvs_string, [] + + digits, version, kind = m.group(1), m.group(2), m.group(3).lower() + candidates = [] + for prefix in _PREFIXES_FOR_KIND[kind]: + for digit_field in _digit_field_candidates(digits): + accession = f"{prefix}{digit_field}" + try: + versions = _get_tx_versions(data_provider, accession, genome_build) + except NotImplementedError: + return hgvs_string, [] # provider can't enumerate, don't guess + if versions: + candidates.append(accession) + + if len(candidates) != 1: + return hgvs_string, [] # unknown or ambiguous, never guess + + accession = candidates[0] + original_head = hgvs_string[:m.end(2)] + fixed_head = f"{accession}{version}" + resolved = fixed_head + hgvs_string[m.end(2):] + return resolved, [HGVSFix( + severity=HGVSFixSeverity.WARNING, + code=HGVSFixCode.ADDED_ACCESSION_PREFIX, + message=( + f"Restored missing accession prefix: '{original_head}' → " + f"'{fixed_head}' (unique match in data provider)" + ), + original=original_head, + fixed=fixed_head, + )] + + # --------------------------------------------------------------------------- # Public: resolve_transcript_version — adjacent-version fallback # --------------------------------------------------------------------------- @@ -540,9 +640,10 @@ def fix_hgvs( """ Clean and resolve an HGVS string in one call. - Runs clean_hgvs() first (always), then resolve_gene_hgvs() if a - data_provider and genome_build are supplied. All fixes from both steps - are returned in a single list. + Runs clean_hgvs() first (always), then resolve_missing_accession_prefix() + if a data_provider is supplied, then resolve_gene_hgvs() if a data_provider + and genome_build are supplied. All fixes from every step are returned in a + single list. If data_provider/genome_build are omitted, only string cleaning is performed. This is useful when the caller knows the input already contains a transcript. @@ -586,6 +687,14 @@ def fix_hgvs( """ result, fixes = clean_hgvs(hgvs_string, ops=ops) + if data_provider is not None: + # Restore a dropped RefSeq prefix (eg "000059.4:c.68del") before gene + # resolution, so the bare number is not mistaken for a gene symbol. + result, prefix_fixes = resolve_missing_accession_prefix( + result, data_provider, genome_build=genome_build, + ) + fixes.extend(prefix_fixes) + if data_provider is not None and genome_build is not None: result, resolution_fixes = resolve_gene_hgvs( result, data_provider, genome_build, diff --git a/claude/20260808_paper_feedback_plan.md b/claude/20260808_paper_feedback_plan.md deleted file mode 100644 index 35eb228..0000000 --- a/claude/20260808_paper_feedback_plan.md +++ /dev/null @@ -1,195 +0,0 @@ -# Paper feedback plan - -Response to the presentation feedback in `claude/feedback.md`. Each item below records what -the paper currently says, what to do about it, and what else it touches. Ordered by item -number; a suggested working order is at the end. - -## 1. The Mutalyzer "50%" claim (introduction) - -`paper/introduction.md` says Mutalyzer "found ~50% error rates in submitted HGVS -descriptions over five years, many attributable to missing transcript data". Lefter 2021 -reports the opposite split: ~50% of submitted descriptions were *correct*, ~41% had a -syntactic or semantic error, and only ~7% could be automatically corrected by Mutalyzer. -The "many attributable to missing transcript data" clause appears unsupported. - -To do: - -- Rewrite the sentence with the real breakdown (50% correct, 41% error, 7% auto-corrected). - Change `literature.csv` and `paper/Snakefile` (the `literature` rule) from the single - `hgvs_error_rate_pct: 50` fact to the three separate facts. -- Expand it into the motivation for R4. An independent, large-scale production stream says - roughly 4 in 10 human-entered HGVS strings are broken, and the leading checker repairs - only ~7% of submissions. That is the gap `clean_hgvs()` fills: we rescue 5.1 points of - the 8.5 that fail, about 60% of the failures. Much stronger framing than the current line. -- Take the numbers and denominators from the paper body (PMC8479679), not the abstract, and - check whether Mutalyzer's rate is per submission or per unique string, since we report both. - -## 2. Figure 1: show UTA in the data providers - -UTA appears once in Figure 1, as an input in panel A (alignments absent from any annotation -file). Panel B shows only `JSONDataProvider`, `RESTDataProvider` and `EnsemblTarkDataProvider`. -That undersells the drop-in-replacement claim and leaves the R2/R3 comparator invisible. - -To do: - -- Add a greyed `UTADataProvider` to `PostgreSQL` node in panel B alongside cdot's providers, - with both plugging into the same `hgvs.dataproviders.interface`. Grey already means - external in the existing legend, so no legend change is needed. -- Update the figure legend in `paper/figures.md` to name UTA as the alternative backend. - -This makes the "swap only the transcript-data layer" experimental design visible at a glance. - -## 3. Literature search on existing HGVS cleaning - -The discussion currently describes VariantValidator and Mutalyzer as validators only. Both -in fact perform repair. VariantValidator auto-corrects where it can, including re-mapping -incorrectly reported intron/exon boundary coordinates. Mutalyzer 2 returns a corrected -description for ~7% of submissions. A third tool is missing from the paper entirely: the -LOVD HGVS syntax checker (LOVDnl/HGVS-syntax-checker), which offers ranked, -confidence-scored corrections of invalid descriptions. - -To do: - -- Search properly: VariantValidator, Mutalyzer 2 and 3, the LOVD syntax checker, - hgvs-weaver, the ClinGen Allele Registry, and the leniency of the biocommons parser - itself. Add a short related-work paragraph to `paper/discussion.md`. -- State the differentiator honestly. Those tools are web services and validators that repair - as a side effect of validation. cdot's contribution is an offline Python function that - repairs before parsing, reports every change as an inspectable `HGVSFix`, and guarantees - no regressions. Position on repair scope, auditability and offline use, not on "they do - not do this". -- Optional and stronger: run the injection corpus (`paper/scripts/inject_and_clean.py`) - through the LOVD syntax checker for a head-to-head recovery comparison. This is the kind - of comparison a Bioinformatics reviewer is likely to ask for, so it may be worth doing - pre-emptively. - -Reading: - -- Mutalyzer 2 (Lefter 2021): https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8479679/ -- VariantValidator (Freeman 2018): https://onlinelibrary.wiley.com/doi/full/10.1002/humu.23348 -- LOVD HGVS syntax checker: https://github.com/LOVDnl/HGVS-syntax-checker - -## 4. ClinVar resolution is too high (highest value item) - -Confirmed from `paper/scripts/build_clinvar_pairs.py`: the c.HGVS comes from the `Name` -column of `variant_summary.txt.gz`, which is ClinVar's own recomputed preferred-transcript -name and is always at the current version. So the 99.4% resolved / 98.8% matched figures -measure the easy case, and cannot exercise the historical-depth claim at all. That is why -the interesting version story currently sits in the Tier 2 (private, non-reproducible) -Shariant corpus. - -The XML parser exists, but not in `too-many-transcripts`. It is -`../clinvar-hidden-structures/scripts/extract_xml_to_csv.py`, which already walks the VCV -XML, pulls `ClinicalAssertion .//AttributeSet/Attribute[@Type="HGVS"]` values, and has -`_extract_coding_transcript()` to grab the versioned accession the submitting lab actually -used. - -To do: - -- Add `paper/scripts/build_clinvar_submitted_pairs.py`: SCV-submitted HGVS (as written) - joined to the VCV VCF coordinate as ground truth, via VariationID / AlleleID. -- Report the version-age distribution. How many submitted strings cite a transcript version - that is no longer current, and how many cite a version retired from the current RefSeq - release. This is the number that makes the historical-depth argument concrete. -- Re-run R2 as cdot vs UTA on this corpus. Expect the gap to widen and the absolute numbers - to drop, which is the point. -- Biggest win in the whole list: submitted SCV strings are also messy, so this yields a - public, reproducible (Tier 1) counterpart to the Tier 2 cleaning result in R4, plus a - public residual-error taxonomy. That would demote the private corpus to a confirmation - rather than the primary evidence for two separate results. -- Keep the existing variant_summary pass as an explicit "current-version ceiling" row, so - the comparison is stated rather than silently swapped. - -Caveats to flag in Methods: the full VCV XML is large, and submitted HGVS is not -deduplicated across SCVs, so the dedup and sampling rules need documenting. - -## 5. REST faster than local is a red flag, not a result - -`empirical_results/benchmark.csv` has local JSON at 540 to 665 HGVS/s and REST after -prefetch at 731. R3 currently explains this away as run-to-run variance. A reviewer is -unlikely to accept a benchmark whose headline configuration loses to the network. - -To do: - -- Replace the single-shot numbers with N=5 repeats per configuration over an identical fixed - HGVS set, reporting median and IQR (or a boxplot as a supplementary figure). This means - updating `paper/scripts/compute_benchmark.py`, which currently does one run of - `N_BENCHMARK=500`, and carrying dispersion into the Table 1 cells. Note also that the - 540 to 665 range and the 731 REST figure came from differently sized runs, which is an - apples-to-oranges problem on its own. -- Investigate two concrete mechanisms before settling on variance as the explanation: - - `AbstractJSONDataProvider.get_tx_exons` (`cdot/hgvs/dataproviders/json_data_provider.py`) - rebuilds the full exon dict list on every call, with no memoisation. Both providers pay - this, and it is the per-variant hot path. An `lru_cache` here is a real, client-visible - speedup and a genuine changelog entry. - - GC pressure. Local JSON holds hundreds of thousands of transcript dicts resident, while - a warmed REST cache holds a few thousand. `gc.freeze()` after load (or `gc.disable()` - for the duration of the benchmark) is a cheap test. If it explains the inversion, that - is both a better sentence for the paper and a real fix. -- If REST genuinely matches local after repeats, say so plainly with the dispersion shown, - and keep the existing explanation (both resolve from an in-memory dict, so the shared - sequence layer bounds both). - -## 6. Fix the dropped `NM_` prefix in `clean_hgvs()` - -Agreed, this is the most tractable slice of the residual. The "bad accession" class is 167 -queries (14.9% of the residual) and dropped prefix is its headline example. `clean.py` -already has `_op_add_transcript_underscore` (`NM000059` to `NM_000059`), so the missing -sibling is the bare-number case. - -To do: - -- Measure first. Run the sub-breakdown of the 167 bad-accession residuals through - `../cdot_private/analyze_cleaning.py` to see how much is genuinely dropped prefix versus - misplaced or truncated version. Implement rules for whatever dominates. -- The dropped-prefix rule can be disambiguated by the kind letter: a bare - `\d{6,9}\.\d+:c.` implies `NM_` or `XM_`, and `:n.` implies `NR_` or `XR_`. The NM/XM - ambiguity is real, so make the rule data-provider aware: try the candidate accessions - against the loaded provider and apply only when exactly one exists. This fits the existing - design (reported as an `HGVSFix`, never silent) and demonstrates something no pure - string-level checker can do, namely that having the transcript data locally enables better - cleaning. -- Add a new `HGVSFixCode`, plus tests using synthesised public examples only (BRCA2 - `NM_000059.4`, RUNX1 `NM_001754.5`). Nothing sourced from `cdot_private`. -- Add a `CHANGELOG.md` entry under `[unreleased]` / `### Added` referencing #112. -- Re-run the production corpus and regenerate Table 2 and Table S6, since both the rescue - count and the residual taxonomy shift. - -## 7. Quantify positional drift along the transcript - -Good hypothesis, and directly computable from data we already have. -`paper/scripts/compute_version_stability.py::_preserved_fraction` already walks the CDS -breakpoint by breakpoint, so binning preserved and total bases by relative CDS position -(deciles) is a small change to an existing loop. - -To do: - -- Emit a per-decile preservation curve across consecutive version bumps, for RefSeq and - Ensembl. Prediction: monotonic decline toward the 3' end. -- State an important nuance up front, because it partly contradicts the intuition. The - current facts show drift is overwhelmingly whole-CDS (RefSeq 3.2% full versus 0.7% - partial), which is a relocation and therefore position independent. The positional effect - can only live inside the partial-drift bucket, so the analysis must be conditioned on - partial drift. Reported honestly this is a better result: most version risk is - all-or-nothing and detectable, and the position-dependent part is the small partial-drift - tail. -- Give the 1,766 incorrect ClinVar projections (`empirical_results/clinvar_vcf_residual.csv`) - the same treatment: bin by relative position in the transcript to see whether errors - concentrate at the 3' end. -- If the effect is clear, it is a supplementary figure plus one sentence in R5. It also gives - the version fallback a usable rule of thumb: a 5' coding variant is safer to substitute - than a 3' one. - -## Suggested order - -1. Item 4 (submitted-HGVS ClinVar corpus). Largest gain, and it feeds items 6 and 7 with a - public corpus. -2. Item 5 (benchmark repeats plus the two mechanism checks). Reviewer critical, and may - produce a genuine performance fix. -3. Items 1 and 3 together (fix the Mutalyzer number, add the repair-tool related work). Both - cheap, both touch the same motivation paragraph, and item 1 is currently a factual error. -4. Item 6 (`NM_` prefix rule), then regenerate the R4 tables. -5. Item 7 (positional drift), then item 2 (figure). Both small. - -Items 1, 2, 3 and 7 are paper only. Items 5 and 6 touch shipped code and need `CHANGELOG.md` -entries. Item 4 is analysis tooling only and does not. diff --git a/tests/test_clean_hgvs.py b/tests/test_clean_hgvs.py index 3219462..50f7ac9 100644 --- a/tests/test_clean_hgvs.py +++ b/tests/test_clean_hgvs.py @@ -213,6 +213,21 @@ def test_all_clean_ops_covers_every_op(): # #112 — genomic accession wedged into the gene-symbol parenthetical slot ("NM_000059.4(NC_000013.11):c.68del", "NM_000059.4:c.68del", {C.DROPPED_GENOMIC_REF_IN_PARENS}), ("NM_001754.5(NG_042763.1):c.1415T>C", "NM_001754.5:c.1415T>C", {C.DROPPED_GENOMIC_REF_IN_PARENS}), + # #112: semicolon/underscore in place of the accession:allele colon + ("NM_000059.4;c.68del", "NM_000059.4:c.68del", {C.FIXED_SEPARATOR_TYPO}), + ("NM_001754.5_c.749G>A", "NM_001754.5:c.749G>A", {C.FIXED_SEPARATOR_TYPO}), + ("NM_000059.4; c.68del", "NM_000059.4:c.68del", {C.FIXED_SEPARATOR_TYPO}), + # #112: dangling dot of an empty (dropped) transcript version + ("NM_001754.:c.749G>A", "NM_001754:c.749G>A", {C.DROPPED_EMPTY_VERSION}), + ("NM_001754.::c..749G>A", "NM_001754:c.749G>A", {C.DROPPED_EMPTY_VERSION}), + # #112: leading colon/dash junk separated from the accession by whitespace + (": NM_000059.4(BRCA2):c.68del", "NM_000059.4(BRCA2):c.68del", {C.STRIPPED_LEADING_JUNK}), + ("- NM_001754.5:c.749G>A", "NM_001754.5:c.749G>A", {C.STRIPPED_LEADING_JUNK}), + # #112: space in place of the reference:kind colon after a gene symbol. + # Must also work for symbols that resemble an accession (2-letter RefSeq + # prefix + digit, eg NR2E3), which reconstruction deliberately refuses. + ("BRCA2 c.68del", "BRCA2:c.68del", {C.FIXED_SEPARATOR_TYPO}), + ("NR2E3 c.349G>A", "NR2E3:c.349G>A", {C.FIXED_SEPARATOR_TYPO}), ] @@ -295,6 +310,34 @@ def test_drop_genomic_ref_in_parens_collapses_genomic_selector_form(): assert C.DROPPED_GENOMIC_REF_IN_PARENS in codes(fixes) +def test_accession_separator_semicolon_inside_allele_untouched(): + # A ';' inside a (balanced-bracket) multi-variant allele is real HGVS + # punctuation, not a separator typo, so the anchored rule must not fire. + s = "NM_000059.4:c.[68del;70A>G]" + cleaned, fixes = clean_hgvs(s, validate=False) + assert cleaned == s + assert not fixes + + +def test_empty_version_requires_colon_after_dot(): + # A real version must never be dropped; the rule only fires on a dangling + # dot immediately before the colon. + s = "NM_001754.5:c.749G>A" + cleaned, fixes = clean_hgvs(s, validate=False) + assert cleaned == s + assert not fixes + + +def test_space_before_kind_leaves_accessions_to_whitespace_strip(): + # An accession head (contains '_'/'.') is glued by the whitespace strip and + # then repaired by reconstruction. The space-separator rule must not fire, + # so the reported fixes stay the same as before the rule existed. + cleaned, fixes = clean_hgvs("nm_000059.4 c.316+5G>A") + assert cleaned == "NM_000059.4:c.316+5G>A" + assert C.STRIPPED_WHITESPACE in codes(fixes) + assert C.FIXED_SEPARATOR_TYPO not in codes(fixes) + + def test_clean_hgvs_protein_suffix(): cleaned, fixes = clean_hgvs("NM_000059.4:c.316+5G>A p.Arg106*") assert cleaned == "NM_000059.4:c.316+5G>A" diff --git a/tests/test_gene_hgvs.py b/tests/test_gene_hgvs.py index d269b81..232c2c9 100644 --- a/tests/test_gene_hgvs.py +++ b/tests/test_gene_hgvs.py @@ -21,6 +21,7 @@ fix_hgvs, rank_transcripts_for_gene, resolve_gene_hgvs, + resolve_missing_accession_prefix, resolve_transcript_version, UnsafeVersionPolicy, ) @@ -436,7 +437,8 @@ def test_parse_versioned_transcript_no_colon(): # --------------------------------------------------------------------------- -# resolve_transcript_version — adjacent-version fallback +# resolve_missing_accession_prefix (#112): restore a dropped RefSeq prefix, +# verified against the data provider (fires only on a unique match) # --------------------------------------------------------------------------- class _VersionStubProvider: @@ -448,6 +450,116 @@ def get_tx_versions(self, accession): return self._version_map.get(accession, []) +def test_missing_prefix_unambiguous_nm(version_provider): + # "c." implies NM_/XM_; only NM_000059 exists, so the fix applies. + resolved, fixes = resolve_missing_accession_prefix( + "000059.4:c.36del", version_provider) + assert resolved == "NM_000059.4:c.36del" + assert len(fixes) == 1 + assert fixes[0].severity == W + assert fixes[0].code == C.ADDED_ACCESSION_PREFIX + assert fixes[0].original == "000059.4" + assert fixes[0].fixed == "NM_000059.4" + + +def test_missing_prefix_kept_underscore(version_provider): + # The user dropped the letters but kept the underscore. + resolved, fixes = resolve_missing_accession_prefix( + "_000059.4:c.36del", version_provider) + assert resolved == "NM_000059.4:c.36del" + assert fixes[0].code == C.ADDED_ACCESSION_PREFIX + + +def test_missing_prefix_ambiguous_candidates_no_fix(): + # Both NM_ and XM_ exist for these digits: never guess. + provider = _VersionStubProvider({"NM_001754": [5], "XM_001754": [1]}) + resolved, fixes = resolve_missing_accession_prefix( + "001754.5:c.749G>A", provider) + assert resolved == "001754.5:c.749G>A" + assert fixes == [] + + +def test_missing_prefix_unknown_accession_no_fix(version_provider): + # No candidate exists in the provider: leave the string unchanged. + resolved, fixes = resolve_missing_accession_prefix( + "001754.5:c.749G>A", version_provider) + assert resolved == "001754.5:c.749G>A" + assert fixes == [] + + +def test_missing_prefix_kind_letter_routes_to_rna_prefixes(): + # "n." implies NR_/XR_, so an NM_-only provider must not match ... + provider = _VersionStubProvider({"NM_001754": [5]}) + resolved, fixes = resolve_missing_accession_prefix( + "001754.5:n.100A>G", provider) + assert resolved == "001754.5:n.100A>G" + assert fixes == [] + # ... and an NR_ match resolves the "n." form only. + provider = _VersionStubProvider({"NR_038196": [1]}) + resolved, fixes = resolve_missing_accession_prefix( + "038196.1:n.100A>G", provider) + assert resolved == "NR_038196.1:n.100A>G" + assert fixes[0].code == C.ADDED_ACCESSION_PREFIX + resolved, fixes = resolve_missing_accession_prefix( + "038196.1:c.100A>G", provider) + assert resolved == "038196.1:c.100A>G" + assert fixes == [] + + +def test_missing_prefix_zero_pads_seven_digit_field(): + # RefSeq digit fields are 6 or 9 digits; a 7-digit input is a 9-digit + # accession that also lost its leading zeros. + provider = _VersionStubProvider({"NM_001128425": [1]}) + resolved, fixes = resolve_missing_accession_prefix( + "1128425.1:c.667A>G", provider) + assert resolved == "NM_001128425.1:c.667A>G" + assert fixes[0].code == C.ADDED_ACCESSION_PREFIX + + +def test_missing_prefix_unversioned_accession(version_provider): + resolved, fixes = resolve_missing_accession_prefix( + "000059:c.36del", version_provider) + assert resolved == "NM_000059:c.36del" + assert fixes[0].code == C.ADDED_ACCESSION_PREFIX + + +def test_missing_prefix_provider_cannot_enumerate_no_fix(): + class _NoVersions: + def get_tx_versions(self, accession): + raise NotImplementedError("nope") + + resolved, fixes = resolve_missing_accession_prefix( + "000059.4:c.36del", _NoVersions()) + assert resolved == "000059.4:c.36del" + assert fixes == [] + + +def test_missing_prefix_noop_without_bare_number(version_provider): + # Complete references and gene symbols are untouched. + for s in ("NM_000059.4:c.36del", "BRCA2:c.36del", "12345.1:c.36del"): + assert resolve_missing_accession_prefix(s, version_provider) == (s, []) + + +def test_fix_hgvs_restores_missing_prefix(version_provider): + # fix_hgvs applies the provider-aware prefix restoration (no genome_build + # needed) after string cleaning. + result, fixes = fix_hgvs(" 000059.4:c.36DEL", data_provider=version_provider) + assert result == "NM_000059.4:c.36del" + fix_codes = {f.code for f in fixes} + assert C.ADDED_ACCESSION_PREFIX in fix_codes + assert C.STRIPPED_WHITESPACE in fix_codes + + +def test_fix_hgvs_no_provider_leaves_bare_number_unchanged(): + result, _fixes = fix_hgvs("000059.4:c.36del") + assert result == "000059.4:c.36del" + + +# --------------------------------------------------------------------------- +# resolve_transcript_version — adjacent-version fallback +# --------------------------------------------------------------------------- + + @pytest.fixture def version_provider(): return _VersionStubProvider({"NM_000059": [2, 3, 4]}) From b9ab0e22cf874b2c5b7cdd056bb288d80e0edbc6 Mon Sep 17 00:00:00 2001 From: Dave Lawrence Date: Mon, 17 Aug 2026 15:04:51 +0000 Subject: [PATCH 2/3] Paper: address presentation feedback (submitted-SCV corpus, benchmarks, positional drift) --- claude/clinvar_diff_coordinates.md | 16 + paper/README.md | 2 +- paper/Snakefile | 175 +++++- paper/discussion.md | 21 +- paper/empirical_results/benchmark.csv | 4 +- .../clinvar_residual_positions.csv | 2 + paper/empirical_results/clinvar_submitted.csv | 2 + .../clinvar_submitted_residual.csv | 2 + paper/empirical_results/literature.csv | 4 +- paper/empirical_results/positional_drift.csv | 2 + paper/figures.md | 5 +- paper/figures/figure1.svg | 61 ++- paper/figures/figure_s1_positional_drift.svg | 101 ++++ paper/introduction.md | 10 +- paper/methods.md | 61 ++- paper/references.bib | 17 + paper/results.md | 257 ++++++--- paper/scripts/bin_residual_positions.py | 122 +++++ .../scripts/build_clinvar_submitted_pairs.py | 241 +++++++++ paper/scripts/compute_benchmark.py | 269 ++++++---- .../scripts/compute_submitted_version_age.py | 203 +++++++ paper/scripts/compute_version_stability.py | 99 +++- paper/scripts/inject_and_clean.py | 31 +- paper/scripts/make_positional_figure.py | 155 ++++++ paper/scripts/resolve_clinvar_pass.py | 40 +- paper/supplementary.md | 56 +- .../clinvar_hgvs/clinvar_submitted_500.tsv | 501 ++++++++++++++++++ 27 files changed, 2198 insertions(+), 261 deletions(-) create mode 100644 paper/empirical_results/clinvar_residual_positions.csv create mode 100644 paper/empirical_results/clinvar_submitted.csv create mode 100644 paper/empirical_results/clinvar_submitted_residual.csv create mode 100644 paper/empirical_results/positional_drift.csv create mode 100644 paper/figures/figure_s1_positional_drift.svg create mode 100644 paper/scripts/bin_residual_positions.py create mode 100644 paper/scripts/build_clinvar_submitted_pairs.py create mode 100644 paper/scripts/compute_submitted_version_age.py create mode 100644 paper/scripts/make_positional_figure.py create mode 100644 tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv diff --git a/claude/clinvar_diff_coordinates.md b/claude/clinvar_diff_coordinates.md index 0e7ad7e..cf4acba 100644 --- a/claude/clinvar_diff_coordinates.md +++ b/claude/clinvar_diff_coordinates.md @@ -105,8 +105,24 @@ artifacts. `clinvar_latest.vcf.gz` 2026-03-21 by ALLELEID; the slight release skew contributes a few wrong-ground-truth rows in `incorrect`. +## Positional distribution (2026-08-17 re-run) + +The original per-variant table was deleted, so the pass was re-run for the R5 +positional-drift analysis (same code/settings, cdot 0.2.33 RefSeq + genome FASTA) on the +Jun 26 rebuild of the pairs file (4,423,358 pairs, slightly newer ClinVar): 4,395,535 +correct / 1,772 incorrect / 3,228 no_data / 22,823 error, i.e. the committed totals +within six variants. `paper/scripts/bin_residual_positions.py` bins each correct and +incorrect row's cited c. position by relative CDS position (deciles). Result: correct +projections are nearly uniform (9.6-10.4% per decile); the 1,718 binned incorrect rows +show no 3' concentration, only a mild mid-CDS excess (14.3% in decile 6) and a depleted +3'-most decile (5.5% vs 9.6% baseline), consistent with the representation/multi-mapping +composition above. Facts committed to +`paper/empirical_results/clinvar_residual_positions.csv`. + ## Files - Full per-variant table (not in source control): `output/clinvar_pass/refseq_full_vcf.csv` + (2026-08-17 re-run: `/data/clinvar/clinvar_pass_rerun/refseq_full_vcf.csv`) - Residual / errors (not in source control): `output/clinvar_pass/vcf_errors.csv` - Headline facts (committed): `paper/empirical_results/clinvar_vcf.csv` +- Positional binning (committed): `paper/empirical_results/clinvar_residual_positions.csv` diff --git a/paper/README.md b/paper/README.md index fc3c37c..5772257 100644 --- a/paper/README.md +++ b/paper/README.md @@ -109,7 +109,7 @@ The full-scale ClinVar throughput runs take ~1.5 h each — see `claude/benchmar - **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.4%), the per-fix rescue distribution (Results Table 1), and the residual + (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 `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`), diff --git a/paper/Snakefile b/paper/Snakefile index 2063668..99f5fdd 100644 --- a/paper/Snakefile +++ b/paper/Snakefile @@ -65,8 +65,9 @@ PDF_FLAG = "--pdf" if config.get("pdf", False) else "" EMPIRICAL = "paper/empirical_results" GEN_FACTS = "output/facts" FACT_FILES = ["literature.csv", "coverage.csv", "benchmark.csv", - "clinvar.csv", "cleaning.csv", "sources.csv", "historical.csv", - "version_stability.csv"] + "clinvar.csv", "clinvar_submitted.csv", "clinvar_submitted_residual.csv", + "cleaning.csv", "sources.csv", "historical.csv", + "version_stability.csv", "positional_drift.csv"] PAPER_SOURCES = [ "paper/abstract.md", @@ -76,6 +77,7 @@ PAPER_SOURCES = [ "paper/discussion.md", "paper/figures.md", "paper/figures/figure1.svg", + "paper/figures/figure_s1_positional_drift.svg", "paper/supplementary.md", "paper/references.bib", "paper/bioinformatics.csl", @@ -149,7 +151,13 @@ rule literature: os.makedirs("output/facts", exist_ok=True) rows = [{ "clinvar_variants": 3000000, - "hgvs_error_rate_pct": 50, + # Lefter 2021 (Mutalyzer): breakdown of the ~26 million unique + # descriptions from five years of production logs (per unique + # description, not per submission). Error = 15.8% syntactic + + # 25.4% semantic failures of the Name Checker. + "mutalyzer_correct_pct": 50.4, + "mutalyzer_error_pct": 41.2, + "mutalyzer_autocorrect_pct": 7.1, "lof_agreement_pct": 44, "brca2_accuracy_pct": 32, "uta_count": 141000, @@ -197,12 +205,15 @@ rule coverage: rule benchmark: - """Benchmark local JSON throughput and load time. - - Requires --config data_dir=/path/to/data pointing to a GRCh38 RefSeq JSON.gz. - Without data, writes the measured throughput facts below (HGVS/s, GRCh38 RefSeq, - identical biocommons engine + shared local SeqRepo; see paper/scripts/benchmark_resolution.py - and paper/scripts/compare_providers.py). These are frozen measurements, not placeholders. + """Benchmark HGVS resolution throughput per transcript backend (Table 1). + + Requires --config data_dir=/path/to/data pointing to a GRCh38 RefSeq JSON.gz + (plus uta_uri=... for the local-UTA row) and HGVS_SEQREPO_DIR in the environment. + paper/scripts/compute_benchmark.py resolves the identical committed 500-pair + ClinVar set through every backend, 5 timed passes each, one shared fd-cached + SeqRepo sequence layer, reporting median and IQR. The public remote UTA rows take + ~20 min (~0.1 HGVS/s); skip them during development with an edited invocation. + Without data, writes the frozen measured facts below (not placeholders). """ output: "output/facts/benchmark.csv" run: @@ -214,17 +225,30 @@ rule benchmark: shell(f"{PYTHON} paper/scripts/compute_benchmark.py --refseq-grch38 {grch38} {uta_flag}") else: print("No data_dir — writing the frozen measured benchmark facts to benchmark.csv") - # Measured HGVS/s (GRCh38 RefSeq, biocommons engine, shared local SeqRepo): - # local JSON 540 (500-set) -> 665 (full 3.66M ClinVar); REST 39 cold, 731 after a - # batch prefetch; local UTA 24; public remote UTA 0.1. See R5 / Table 2. + # Measured 2026-08-17 (cdot 0.2.34 GRCh38 RefSeq, uta_20241220 local, public + # uta.biocommons.org remote, cdotlib.org REST, biocommons engine, one shared + # local SeqRepo with fd caching): median (q1, q3) HGVS/s over 5 passes of the + # committed 500-pair set (remote UTA: first 25 pairs). See R3 / Table 1. row = { - "cdot_local_min_tps": 540, - "cdot_local_max_tps": 665, - "cdot_rest_tps": 39, - "cdot_rest_prefetch_tps": 731, - "uta_local_tps": 24, - "uta_remote_tps": 0.1, - "grch38_load_time_s": 10.6, + "n_pairs": 500, + "n_repeats": 5, + "cdot_local_tps": 748, + "cdot_local_tps_q1": 743, + "cdot_local_tps_q3": 752, + "grch38_load_time_s": 4.8, + "cdot_rest_tps": 16, + "cdot_rest_tps_q1": 15, + "cdot_rest_tps_q3": 16, + "cdot_rest_prefetch_tps": 753, + "cdot_rest_prefetch_tps_q1": 752, + "cdot_rest_prefetch_tps_q3": 753, + "uta_local_tps": 187, + "uta_local_tps_q1": 187, + "uta_local_tps_q3": 188, + "uta_remote_tps": 0.12, + "uta_remote_tps_q1": 0.12, + "uta_remote_tps_q3": 0.12, + "uta_remote_n": 25, } with open(output[0], "w", newline="") as fh: writer = csv.DictWriter(fh, fieldnames=list(row)) @@ -292,6 +316,105 @@ rule clinvar: print(f"Written: {output[0]}") +rule clinvar_submitted: + """R2 submitted-HGVS ClinVar corpus facts: version-age distribution plus the + cdot-vs-UTA comparison on strings as the submitting labs wrote them. + + The variant_summary-based corpora take c.HGVS from ClinVar's recomputed Name + column, always at the current transcript version, so they are a current-version + ceiling. This corpus instead uses the per-SCV HGVS attributes of the VCV XML + (submitted strings verbatim), joined to the ClinVar VCF coordinate via AlleleID. + Frozen measured constants; the corpus (built from a ClinVar XML download) is not + committed, but a 500-pair seed-42 sample is + (tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv). Provenance + (measured 2026-08-17, ClinVarVCVRelease_2026-06, cdot 0.2.34 refseq GRCh38, + local uta_20241220): + # corpus: 5,652,560 SCV HGVS values -> 3,495,275 transcript c./n. strings + # -> 2,933,667 unique (AlleleID, string) pairs; 100.00% RefSeq, 0 ENST + python paper/scripts/build_clinvar_submitted_pairs.py \\ + --xml ClinVarVCVRelease_2026-06.xml.gz clinvar.GRCh38.vcf.gz \\ + clinvar_submitted_pairs.GRCh38.tsv # (or --scv-csv-dir extraction) + # version age vs the current annotation release (RS_2025_08, auto-detected + # from per-transcript source URLs in the cdot JSON): + python paper/scripts/compute_submitted_version_age.py \\ + clinvar_submitted_pairs.GRCh38.tsv \\ + --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz \\ + --refseq-allbuilds cdot-0.2.34.all-builds-refseq-....json.gz + # resolution on the seed-42 3,000-pair sample (VCF-coordinate scoring): + F=GCF_000001405.39_GRCh38.p13_genomic.fna.gz + python paper/scripts/resolve_clinvar_pass.py clinvar_submitted_sample3000.tsv \\ + --json cdot-0.2.34.refseq.GRCh38.json.gz --fasta $F --with-fixes \\ + --out submitted_pass_cdot_3000.csv + UTA_DB_URL=... HGVS_SEQREPO_DIR=... python paper/scripts/resolve_clinvar_pass.py \\ + clinvar_submitted_sample3000.tsv --uta --out submitted_pass_uta_3000.csv + Residual taxonomy (clinvar_submitted_residual.csv), derived from the cdot pass + rows with fixed_bucket != correct (37 of 3,000): + * version_refused (26): cited version absent from the data; the adjacent-version + fallback declined to substitute because coordinate-safety could not be + verified (REFUSED_UNSAFE_VERSION; no false rescues by design). Split from + no_data via summarize_clinvar_pass.py --split-no-data (26 unknown-version). + * unknown_accession (1): no version of the accession in the data. + * coordinate_drift (4): resolves through the cited historical version to a + coordinate that differs from ClinVar's current interpretation. + * position_out_of_bounds (3) / reference_mismatch (1): the cited position or + base does not exist on the cited version (raises HGVSInvalidIntervalError / + HGVSInvalidVariantError). + * grammar_unsupported (2): repeat ref[N] and allele [..] notation the + biocommons grammar rejects. + """ + output: + facts = "output/facts/clinvar_submitted.csv", + residual = "output/facts/clinvar_submitted_residual.csv", + run: + import os, csv + os.makedirs("output/facts", exist_ok=True) + row = { + "n_scv_tx_strings": 3495275, # transcript c./n. strings across SCVs + "n_unique_pairs": 2933667, # unique (AlleleID, submitted string) + "ensembl_pct": 0.0, # not one submitted ENST string + "version_not_current_pct": 81.8, # cited version != current (RS_2025_08) + "scv_weighted_not_current_pct": 80.3, # same, weighted by scv_count + "base_retired_pct": 0.7, # accession absent from current release + "not_current_in_cdot_pct": 99.3, # cdot GRCh38 holds the superseded version + "absent_cdot_pct": 0.6, # cited version absent from cdot GRCh38 + "n_sample": 3000, # seed-42 uniform sample + "sample_seed": 42, + "cdot_resolved_pct": 98.9, # 2966/3000 + "cdot_matched_pct": 98.7, # 2962/3000 reproduce the VCF coordinate + "cdot_no_data": 27, # 26 unknown-version + 1 unknown-accession + "cdot_incorrect": 4, + "cdot_error": 7, + "uta_resolved_pct": 80.1, # 2404/3000 (all matched) + "uta_no_data_pct": 19.7, # 592/3000: no alignment for the version + "cdot_only": 563, # resolved by cdot, not UTA + "cdot_only_pct": 18.8, + "uta_only": 1, + "rescued_by_fix": 1, # dup with length suffix (DROPPED_DEL_DUP_COUNT) + "regressions": 0, + "after_fix_matched_pct": 98.8, # 2963/3000 + "residual_n": 37, + "residual_pct": 1.2, + } + with open(output.facts, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=list(row)) + writer.writeheader() + writer.writerow(row) + resid = { + "n_residual": 37, + "version_refused": 26, + "unknown_accession": 1, + "coordinate_drift": 4, + "position_out_of_bounds": 3, + "reference_mismatch": 1, + "grammar_unsupported": 2, + } + with open(output.residual, "w", newline="") as fh: + writer = csv.DictWriter(fh, fieldnames=list(resid)) + writer.writeheader() + writer.writerow(resid) + print(f"Written: {output.facts}, {output.residual}") + + rule cleaning: """Tier-1 reproducible injection cleaning benchmark (issue #112). @@ -388,10 +511,18 @@ rule version_stability: stability improvement from Ensembl 116 entering the data (verified by re-running on the 0.2.33 files: RefSeq matched the old snapshot exactly, Ensembl already matched the new numbers). + + The same script run also bins preserved coding bases by relative CDS position + (deciles, 5'->3') into positional_drift.csv, conditioned on the partial-drift + pairs and unconditioned (see the R5 positional-drift figure); its frozen + constants live in paper/empirical_results/positional_drift.csv and are copied + here when no data is present. """ - output: "output/facts/version_stability.csv" + output: + "output/facts/version_stability.csv", + "output/facts/positional_drift.csv", run: - import os, csv + import os, csv, shutil os.makedirs("output/facts", exist_ok=True) rs38 = FILES.get("refseq_grch38", "") en38 = FILES.get("ensembl_grch38", "") @@ -406,6 +537,8 @@ rule version_stability: shell(f"{PYTHON} paper/scripts/compute_version_stability.py {' '.join(flags)}") else: print("No data_dir — writing the frozen measured version-stability facts") + shutil.copy(f"{EMPIRICAL}/positional_drift.csv", + "output/facts/positional_drift.csv") row = { "sample_n": 12000, "build": "GRCh38", diff --git a/paper/discussion.md b/paper/discussion.md index a8427d5..971d352 100644 --- a/paper/discussion.md +++ b/paper/discussion.md @@ -29,8 +29,25 @@ adds RefSeq coverage, fully offline operation, and support for T2T-CHM13v2.0. To as VariantValidator [@Freeman2018], built on the biocommons/hgvs library with a self-hosted copy of UTA, and Mutalyzer [@Lefter2021], which uses its own independent normalisation stack and retrieves transcripts directly from NCBI and Ensembl, are widely -used to check HGVS correctness; cdot is complementary to them, supplying the -transcript-coordinate layer rather than validating descriptions. +used to check and correct HGVS descriptions; cdot is complementary to them, supplying +the transcript-coordinate layer rather than a validation service. + +cdot is also not the first tool to repair broken HGVS rather than merely reject it. +VariantValidator automatically corrects the mistakes it can interpret, including +intronic positions described relative to the wrong exon boundary [@Freeman2018]. +Mutalyzer's Name Checker returned a corrected description for +~{{ literature.mutalyzer_autocorrect_pct | dp(0) }}% of the unique descriptions in its +production logs [@Lefter2021]. The LOVD HGVS syntax checker [@LovdHgvsChecker] checks +syntax without needing a reference sequence and suggests corrections for invalid +descriptions, ranked by likelihood; the ClinGen Allele Registry [@Pawliczek2018] +normalises the descriptions it registers to canonical allele identifiers; and the +biocommons/hgvs parser itself tolerates a few common deviations, such as a gene symbol +in parentheses after the accession. For these tools repair happens as part of validation +or registration, usually behind a web service. `clean_hgvs()` differs in role and +placement: it is an offline Python function that repairs the string before parsing, so +anything that consumes the string downstream benefits; every change is returned as an +`HGVSFix` the caller can audit; and it guarantees no regressions, never breaking a +description that already parsed. Beyond HGVS resolution, the JSON representation is useful in its own right. It parses far faster than the GTF/GFF files it is built from and loads trivially over HTTP, so cdot diff --git a/paper/empirical_results/benchmark.csv b/paper/empirical_results/benchmark.csv index 7487b56..449cb5b 100644 --- a/paper/empirical_results/benchmark.csv +++ b/paper/empirical_results/benchmark.csv @@ -1,2 +1,2 @@ -cdot_local_min_tps,cdot_local_max_tps,cdot_rest_tps,cdot_rest_prefetch_tps,uta_local_tps,uta_remote_tps,grch38_load_time_s -540,665,39,731,24,0.1,10.6 +n_pairs,n_repeats,cdot_local_tps,cdot_local_tps_q1,cdot_local_tps_q3,grch38_load_time_s,cdot_rest_tps,cdot_rest_tps_q1,cdot_rest_tps_q3,cdot_rest_prefetch_tps,cdot_rest_prefetch_tps_q1,cdot_rest_prefetch_tps_q3,uta_local_tps,uta_local_tps_q1,uta_local_tps_q3,uta_remote_tps,uta_remote_tps_q1,uta_remote_tps_q3,uta_remote_n +500,5,748,743,752,4.8,16,15,16,753,752,753,187,187,188,0.12,0.12,0.12,25 diff --git a/paper/empirical_results/clinvar_residual_positions.csv b/paper/empirical_results/clinvar_residual_positions.csv new file mode 100644 index 0000000..592c010 --- /dev/null +++ b/paper/empirical_results/clinvar_residual_positions.csv @@ -0,0 +1,2 @@ +correct_binned,correct_decile1_pct,correct_decile2_pct,correct_decile3_pct,correct_decile4_pct,correct_decile5_pct,correct_decile6_pct,correct_decile7_pct,correct_decile8_pct,correct_decile9_pct,correct_decile10_pct,correct_5putr,correct_3putr,correct_noncoding_n,correct_no_cds_len,correct_unparsed,incorrect_binned,incorrect_decile1_pct,incorrect_decile2_pct,incorrect_decile3_pct,incorrect_decile4_pct,incorrect_decile5_pct,incorrect_decile6_pct,incorrect_decile7_pct,incorrect_decile8_pct,incorrect_decile9_pct,incorrect_decile10_pct,incorrect_5putr,incorrect_3putr,incorrect_noncoding_n,incorrect_no_cds_len,incorrect_unparsed +4285669,10.4,10.0,9.9,9.9,10.0,10.0,10.0,10.1,10.1,9.6,36536,68642,4688,0,0,1718,9.9,9.3,8.0,8.5,11.2,14.3,10.5,12.9,9.9,5.5,28,25,1,0,0 diff --git a/paper/empirical_results/clinvar_submitted.csv b/paper/empirical_results/clinvar_submitted.csv new file mode 100644 index 0000000..a552f88 --- /dev/null +++ b/paper/empirical_results/clinvar_submitted.csv @@ -0,0 +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_no_data_pct,cdot_only,cdot_only_pct,uta_only,rescued_by_fix,regressions,after_fix_matched_pct,residual_n,residual_pct +3495275,2933667,0.0,81.8,80.3,0.7,99.3,0.6,3000,42,98.9,98.7,27,4,7,80.1,19.7,563,18.8,1,1,0,98.8,37,1.2 diff --git a/paper/empirical_results/clinvar_submitted_residual.csv b/paper/empirical_results/clinvar_submitted_residual.csv new file mode 100644 index 0000000..f227239 --- /dev/null +++ b/paper/empirical_results/clinvar_submitted_residual.csv @@ -0,0 +1,2 @@ +n_residual,version_refused,unknown_accession,coordinate_drift,position_out_of_bounds,reference_mismatch,grammar_unsupported +37,26,1,4,3,1,2 diff --git a/paper/empirical_results/literature.csv b/paper/empirical_results/literature.csv index 78907b7..1f6366b 100644 --- a/paper/empirical_results/literature.csv +++ b/paper/empirical_results/literature.csv @@ -1,2 +1,2 @@ -clinvar_variants,hgvs_error_rate_pct,lof_agreement_pct,brca2_accuracy_pct,uta_count,uta_remote_tps,mane_coverage_pct,seqrepo_speedup_fold -3000000,50,44,32,141000,1,97,1300 +clinvar_variants,mutalyzer_correct_pct,mutalyzer_error_pct,mutalyzer_autocorrect_pct,lof_agreement_pct,brca2_accuracy_pct,uta_count,uta_remote_tps,mane_coverage_pct,seqrepo_speedup_fold +3000000,50.4,41.2,7.1,44,32,141000,1,97,1300 diff --git a/paper/empirical_results/positional_drift.csv b/paper/empirical_results/positional_drift.csv new file mode 100644 index 0000000..b844c06 --- /dev/null +++ b/paper/empirical_results/positional_drift.csv @@ -0,0 +1,2 @@ +sample_n,build,refseq_partial_pairs,refseq_partial_bases,refseq_partial_decile1_pct,refseq_partial_decile2_pct,refseq_partial_decile3_pct,refseq_partial_decile4_pct,refseq_partial_decile5_pct,refseq_partial_decile6_pct,refseq_partial_decile7_pct,refseq_partial_decile8_pct,refseq_partial_decile9_pct,refseq_partial_decile10_pct,refseq_partial_5p_half_pct,refseq_partial_3p_half_pct,refseq_all_decile1_pct,refseq_all_decile2_pct,refseq_all_decile3_pct,refseq_all_decile4_pct,refseq_all_decile5_pct,refseq_all_decile6_pct,refseq_all_decile7_pct,refseq_all_decile8_pct,refseq_all_decile9_pct,refseq_all_decile10_pct,refseq_all_5p_half_pct,refseq_all_3p_half_pct,ensembl_partial_pairs,ensembl_partial_bases,ensembl_partial_decile1_pct,ensembl_partial_decile2_pct,ensembl_partial_decile3_pct,ensembl_partial_decile4_pct,ensembl_partial_decile5_pct,ensembl_partial_decile6_pct,ensembl_partial_decile7_pct,ensembl_partial_decile8_pct,ensembl_partial_decile9_pct,ensembl_partial_decile10_pct,ensembl_partial_5p_half_pct,ensembl_partial_3p_half_pct,ensembl_all_decile1_pct,ensembl_all_decile2_pct,ensembl_all_decile3_pct,ensembl_all_decile4_pct,ensembl_all_decile5_pct,ensembl_all_decile6_pct,ensembl_all_decile7_pct,ensembl_all_decile8_pct,ensembl_all_decile9_pct,ensembl_all_decile10_pct,ensembl_all_5p_half_pct,ensembl_all_3p_half_pct +12000,GRCh38,138,306855,98.9,87.9,77.5,69.5,65.3,63.6,58.8,53.1,49.0,37.5,79.8,52.4,97.1,97.0,96.9,96.9,96.8,96.8,96.8,96.7,96.7,96.6,96.9,96.7,49,92270,93.2,75.3,66.0,58.1,50.3,44.3,36.1,30.4,24.5,17.1,68.6,30.5,98.9,98.8,98.7,98.7,98.7,98.6,98.6,98.6,98.6,98.5,98.8,98.6 diff --git a/paper/figures.md b/paper/figures.md index 73c7ded..6ab679b 100644 --- a/paper/figures.md +++ b/paper/figures.md @@ -11,7 +11,10 @@ JSON, hosted on GitHub and also loaded into the cdot_rest REST API possibly malformed, is repaired by `clean_hgvs()`, then parsed and mapped between c. and g. coordinates by biocommons/hgvs, which draws exons and alignments from one of cdot's data providers (in-memory `JSONDataProvider`, -`RESTDataProvider`, or `EnsemblTarkDataProvider`) and transcript sequence from +`RESTDataProvider`, or `EnsemblTarkDataProvider`); the PostgreSQL-backed +`UTADataProvider` built into biocommons/hgvs plugs into the same interface and +is the alternative backend cdot is compared against in Results. Transcript +sequence comes from the sequence layer (SeqRepo and/or `FastaSeqFetcher`, tried in order via `ChainedSeqFetcher`). Opt-in helpers rewrite the input string: `resolve_gene_hgvs()` maps a gene symbol to its MANE transcript, and diff --git a/paper/figures/figure1.svg b/paper/figures/figure1.svg index 7f984f9..523cf06 100644 --- a/paper/figures/figure1.svg +++ b/paper/figures/figure1.svg @@ -104,27 +104,38 @@ - cdot data providers (biocommons/hgvs Interface) - - - JSONDataProvider - in-memory, - interval trees - 540-665 var/s - - - EnsemblTark - DataProvider - Ensembl TARK - web service - - - RESTDataProvider - per-transcript - requests + - prefetch() - - also: JSONPyHGVSTranscriptFactory → PyHGVS (legacy) + data providers (biocommons/hgvs Interface) + + + JSONDataProvider + in-memory, + interval trees + 540-665 var/s + + + EnsemblTark + DataProvider + Ensembl TARK + web service + + + UTADataProvider + hgvs built-in, + SQL queries + + + + + PostgreSQL + UTA database + + + RESTDataProvider + per-transcript + requests + + prefetch() + + also: JSONPyHGVSTranscriptFactory → PyHGVS (legacy) @@ -161,11 +172,11 @@ c. → g. or g. → c. - + download, load ~10 s - - HTTPS, - prefetch() + + HTTPS, + prefetch() diff --git a/paper/figures/figure_s1_positional_drift.svg b/paper/figures/figure_s1_positional_drift.svg new file mode 100644 index 0000000..e93d9ac --- /dev/null +++ b/paper/figures/figure_s1_positional_drift.svg @@ -0,0 +1,101 @@ + + +RefSeqEnsembl +Coding bases preserved (%) +A  Partial-drift bumps only +n = 138 RefSeq / 49 Ensembl pairs + +0 + +25 + +50 + +75 + +100 + +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +CDS position decile (5′ → 3′) + + + + + + + + + + + +RefSeq 38% + + + + + + + + + + + +Ensembl 17% +B  All version bumps +unconditioned; whole-CDS drift dominates + +0 + +25 + +50 + +75 + +100 + +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +CDS position decile (5′ → 3′) + + + + + + + + + + + + + + + + + + + + + + + diff --git a/paper/introduction.md b/paper/introduction.md index 4a790c7..edf7e73 100644 --- a/paper/introduction.md +++ b/paper/introduction.md @@ -11,10 +11,12 @@ the matching transcript-version alignment for the build in question. The scale i ClinVar alone provides HGVS descriptions for >{{ literature.clinvar_variants | commas }} variants [@Landrum2025]. Resolving HGVS descriptions reliably is hard, especially for the human-entered strings that reach search boxes and importers: -Mutalyzer, whose users type descriptions into a free-text box, found -~{{ literature.hgvs_error_rate_pct | dp(0) }}% error rates in -submitted HGVS descriptions over five years [@Lefter2021], many attributable to missing -transcript data. Transcript choice also has downstream consequences: only +of the ~26 million unique descriptions submitted to Mutalyzer over five years of +production logs, only ~{{ literature.mutalyzer_correct_pct | dp(0) }}% were correct as +written, ~{{ literature.mutalyzer_error_pct | dp(0) }}% contained a syntactic or +semantic error, and Mutalyzer could automatically correct only +~{{ literature.mutalyzer_autocorrect_pct | dp(0) }}% (all rates per unique +description, not per submission) [@Lefter2021]. Transcript choice also has downstream consequences: only {{ literature.lof_agreement_pct | dp(0) }}% of putative loss-of-function variants were classified as loss-of-function by both RefSeq and Ensembl annotation sets in ANNOVAR [@McCarthy2014], so which annotation source you pick changes the answer and covering both diff --git a/paper/methods.md b/paper/methods.md index 7fb4935..bce594c 100644 --- a/paper/methods.md +++ b/paper/methods.md @@ -122,6 +122,17 @@ integer length instead of the inserted sequence); these mark incomplete input ra than formatting noise. By default cleaning never raises and always returns its best attempt, though callers may opt into raising on the first error. +One repair goes beyond string manipulation. A bare-number accession, where the user +dropped the RefSeq prefix entirely (`000059.4:c.68del`), cannot be repaired from the +string alone: the kind letter narrows the possibilities (`c.` implies `NM_` or `XM_`, +`n.` implies `NR_` or `XR_`, with 7-8 digit fields zero-padded to the 9-digit form) but +cannot fully disambiguate them. `resolve_missing_accession_prefix()`, applied by +`fix_hgvs()` whenever a data provider is supplied, generates the candidate accessions +and checks each against the loaded transcript data, restoring the prefix only when +exactly one candidate exists there; with zero or several matches the string is left +unchanged. Like every other repair it is reported as an `HGVSFix`, never applied +silently. + A separate, opt-in helper, `get_best_transcript_version()`, addresses transcript-version drift. Unlike the cleaning operations above, which are unambiguous formatting corrections, substituting a different transcript version is a heuristic that @@ -142,9 +153,9 @@ cdot's client stack (Figure 1B) offers three data providers behind the same inte **Local JSON**: `JSONDataProvider` loads a JSON.gz file into memory on initialisation (typically ~{{ benchmark.grch38_load_time_s | dp(0) }} seconds for GRCh38 RefSeq), lazily building interval trees for region queries on first use and dictionaries for -transcript and gene lookup. Transcript retrieval is then O(1), giving throughput of -{{ benchmark.cdot_local_min_tps | commas }}–{{ benchmark.cdot_local_max_tps | commas }} -transcripts/second (Results, Table 1). +transcript and gene lookup. Transcript retrieval is then O(1), giving end-to-end +resolution throughput of ~{{ benchmark.cdot_local_tps | commas }} HGVS/second +(Results, Table 1). **REST API**: `cdot_rest` (https://github.com/SACGF/cdot_rest) serves the same JSON data at cdotlib.org. `RESTDataProvider` fetches transcripts one request per transcript @@ -216,5 +227,45 @@ reproducible control, with `inject_and_clean.py`, which injects each fix categor clean ClinVar strings. Version-fallback safety is measured by `compute_version_stability.py` on GRCh38, using a seeded {{ version_stability.sample_n | commas }}-accession sample drawn from accessions cdot holds at two or more versions (the only accessions where a version -bump can be assessed). Throughput comparisons hold the sequence layer constant (a shared -local SeqRepo) so that only the transcript-data layer varies (Results R3). +bump can be assessed). The same run bins preserved coding bases by relative CDS position +(deciles, 5' to 3') to test whether drift within a version bump is positional +(Supplementary Figure S1). + +The submitted-string ClinVar corpus (Results R2) is built by +`build_clinvar_submitted_pairs.py` from the per-submission HGVS attributes of a ClinVar +VCV XML release (ClinVarVCVRelease_2026-06): for each ClinicalAssertion (SCV), the +first RefSeq/Ensembl transcript c./n. expression among its `Attribute[@Type="HGVS"]` +values is kept verbatim and joined to the variant's VCF coordinate through the record's +AlleleID, the same ground-truth join as `build_clinvar_pairs.py`. 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. Benchmark samples are drawn uniformly at random with a fixed +seed (42): 3,000 pairs for the cdot-versus-UTA comparison (sized to match the seeded +current-version sample), and a 500-pair sample committed to +`tests/test_data/clinvar_hgvs/`. 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. Two practical caveats: the full VCV XML is tens of +gigabytes uncompressed, so the builder streams it (and also accepts a pre-extracted +per-SCV table); and the first-attribute selection is what defines "the" submitted +string when an SCV carries several HGVS expressions. Transcript version age is +computed by `compute_submitted_version_age.py` against the released cdot RefSeq JSON, +whose per-transcript source URL identifies whether an entry comes from the current +annotation release (RS_2025_08) or survives only through cdot's historical releases. + +Throughput (Table 1) is measured by `compute_benchmark.py`. Every configuration resolves +the identical committed set of {{ benchmark.n_pairs | commas }} ClinVar (g.HGVS, c.HGVS) +pairs through the same biocommons/hgvs engine, with the sequence layer held constant: +a single shared local SeqRepo instance, with its file-descriptor cache enabled, serves +every configuration, so the only thing that varies between rows is the transcript-data +layer. (This matters in practice: without file-descriptor caching, SeqRepo re-opens a +compressed FASTA file on every fetch and the sequence layer, not the transcript layer, +bounds every backend.) Each configuration is timed over {{ benchmark.n_repeats | int }} passes +of the set and reported as median (IQR). The timed portion is steady-state resolution +throughput (HGVS parse plus coordinate projection per string); provider setup (JSON load, +database connection), the REST cache-warming `prefetch()`, and one untimed warm-up pass +per configuration are excluded. The REST rows run against the production server at +cdotlib.org over the public internet, so they include real network conditions. The public +remote UTA row is measured on the first {{ benchmark.uta_remote_n | int }} pairs of the same +set, because at ~{{ benchmark.uta_remote_tps | dp(1) }} HGVS/s a full-size pass is +impractical (Results R3). diff --git a/paper/references.bib b/paper/references.bib index 257fa10..2f9a0a5 100644 --- a/paper/references.bib +++ b/paper/references.bib @@ -36,6 +36,13 @@ @misc{FerroHgvs note = {Accessed 2026}, } +@misc{LovdHgvsChecker, + author = {{LOVD}}, + title = {The {LOVD} {HGVS} syntax checker: a library to validate variant descriptions in the {HGVS} nomenclature syntax}, + howpublished = {\url{https://github.com/LOVDnl/HGVS-syntax-checker}}, + note = {Accessed 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}}, @@ -206,6 +213,16 @@ @article{Park2022 doi = {10.1093/labmed/lmab074}, } +@article{Pawliczek2018, + author = {Pawliczek, Piotr and Patel, Ronak Y. and Ashmore, Lawrence R. and others}, + title = {{ClinGen} {Allele} {Registry} links information about genetic variants}, + journal = {Human Mutation}, + year = {2018}, + volume = {39}, + pages = {1690--1701}, + doi = {10.1002/humu.23637}, +} + @article{Pozo2022, author = {Pozo, Fernando and Martinez-Gomez, Laura and Walsh, Tom A. and others}, title = {Assessing the functional relevance of splice isoforms}, diff --git a/paper/results.md b/paper/results.md index 4072115..752bfef 100644 --- a/paper/results.md +++ b/paper/results.md @@ -68,6 +68,59 @@ dominated by paralog and copy-number transcripts that map to more than one genom and by indel-representation differences, not by coordinate errors; the per-source split and the residual breakdown are in Supplementary Table S4. +Both corpora above take their c.HGVS from the `Name` column of ClinVar's +variant_summary, ClinVar's own recomputed preferred-transcript name, which is always at +the current transcript version. They therefore measure a current-version ceiling and +cannot exercise historical transcript depth. To measure what laboratories actually +write, we built a second public corpus from the per-submission (SCV) HGVS attributes of +the ClinVar VCV XML: each submitted string kept verbatim and joined to the variant's VCF +coordinate as ground truth via its AlleleID +({{ clinvar_submitted.n_unique_pairs | commas }} unique submitted-string/variant pairs +from {{ clinvar_submitted.n_scv_tx_strings | commas }} SCV transcript expressions; +construction, dedup and sampling rules in Methods). The corpus is entirely RefSeq (not +one submitted string cites an Ensembl transcript), and its version profile confirms that +submitted traffic is historical: +{{ clinvar_submitted.version_not_current_pct | dp(1) }}% of submitted strings cite a +transcript version that is no longer the version in the current RefSeq annotation +release ({{ clinvar_submitted.scv_weighted_not_current_pct | dp(1) }}% weighted by +submission count), and {{ clinvar_submitted.base_retired_pct | dp(1) }}% cite a +transcript no longer annotated at any version. cdot's merged release history holds +{{ clinvar_submitted.not_current_in_cdot_pct | dp(1) }}% of those superseded versions; +only {{ clinvar_submitted.absent_cdot_pct | dp(1) }}% of cited versions are absent from +its GRCh38 data. + +On a fixed-seed sample of {{ clinvar_submitted.n_sample | commas }} submitted pairs (the +comparison is again gated by UTA throughput), cdot resolved +{{ clinvar_submitted.cdot_resolved_pct | dp(1) }}% and reproduced ClinVar's VCF +coordinate for {{ clinvar_submitted.cdot_matched_pct | dp(1) }}%, versus +{{ clinvar_submitted.uta_resolved_pct | dp(1) }}% for the same locally loaded UTA. On +submitted rather than recomputed strings, the RefSeq gap that is 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: UTA holds no GRCh38 +alignment for {{ clinvar_submitted.uta_no_data_pct | dp(1) }}% of the cited versions, +and {{ clinvar_submitted.cdot_only | commas }} of the +{{ clinvar_submitted.n_sample | commas }} pairs resolve through cdot alone (one through +UTA alone). The submitted strings are largely well-formed, so string cleaning has little +to rescue here (`fix_hgvs()` repaired {{ clinvar_submitted.rescued_by_fix | int }} +string, a dup carrying a length suffix, with +{{ clinvar_submitted.regressions | int }} regressions); the residual +{{ clinvar_submitted.residual_pct | dp(1) }}% after cleaning and version fallback +({{ clinvar_submitted.residual_n | int }} of {{ clinvar_submitted.n_sample | commas }}) +is dominated by version effects, not formatting: +{{ clinvar_submitted_residual.version_refused | int }} cite a version absent from the +data where the fallback declines to substitute because coordinate safety cannot be +verified (no false rescues, by design), +{{ clinvar_submitted_residual.coordinate_drift | int }} resolve through the cited +historical version to a coordinate that differs from ClinVar's current interpretation, +{{ clinvar_submitted_residual.position_out_of_bounds | int }} + +{{ clinvar_submitted_residual.reference_mismatch | int }} cite a position or base that +does not exist on the cited version, +{{ clinvar_submitted_residual.grammar_unsupported | int }} use repeat or allele notation +the biocommons grammar rejects, and +{{ clinvar_submitted_residual.unknown_accession | int }} cites an accession absent +entirely. The corpus builder, version-age script, resolution harness and a committed +500-pair sample make this benchmark reproducible end to end. + **[Tier 2].** The same gap holds on the historical clinical data that motivated cdot. We resolved the complete set of {{ historical.n_lines | commas }} unique HGVS descriptions imported into the Australian @@ -84,60 +137,78 @@ versus {{ historical.uta_resolved_pct | dp(1) }}% for the same locally loaded UT ({{ historical.cdot_only_pct | dp(1) }}% of the corpus), {{ historical.cdot_only_historical_pct | dp(0) }}% were RefSeq transcript versions for which UTA holds no GRCh38 alignment and {{ historical.cdot_only_ensembl_pct | dp(0) }}% -were Ensembl transcripts (absent from UTA entirely). The ClinVar comparison above shows -the GTF→JSON→biocommons pipeline works at scale on current transcript versions; this -corpus shows the historical depth matters in practice: it is the -older-version traffic a working clinical lab generates. There is no ground-truth genomic +were Ensembl transcripts (absent from UTA entirely). The public submitted-string corpus +above establishes the version-age effect reproducibly; this private corpus confirms it +on the traffic a working clinical lab generates, where the same historical-version gap +appears at comparable size. There is no ground-truth genomic coordinate for the private corpus, so the metric is resolution rate rather than correctness (Methods, data availability). ## R3: Throughput -Backends were compared with the sequence layer held constant, so the only thing that -varies across rows of Table 1 is the transcript-data layer. - -**Table 1. End-to-end HGVS resolution throughput by transcript backend**, sequence layer -held constant (shared local SeqRepo), identical biocommons/hgvs engine. - -| Configuration | Throughput (HGVS/s) | +Backends were compared over the identical committed {{ benchmark.n_pairs | commas }}-pair +ClinVar set with the engine and sequence layer held constant (Methods), so the only thing +that varies across rows of Table 1 is the transcript-data layer. Each configuration was +timed over {{ benchmark.n_repeats | int }} passes of the set. + +**Table 1. End-to-end HGVS resolution throughput by transcript backend**: median (IQR) +HGVS/s over {{ benchmark.n_repeats | int }} timed passes of the identical +{{ benchmark.n_pairs | commas }}-pair ClinVar set, sequence layer held constant (one +shared local SeqRepo instance), identical biocommons/hgvs engine. Timing covers +steady-state resolution only (provider load, REST `prefetch()` and one warm-up pass per +configuration are untimed). The public remote UTA row uses the first +{{ benchmark.uta_remote_n | int }} pairs of the same set (a full-size pass is impractical at +this throughput). + +| Configuration | Throughput (HGVS/s), median (IQR) | |---|---| -| UTA: public remote database | ~{{ benchmark.uta_remote_tps | dp(1) }} | -| UTA: local PostgreSQL | ~{{ benchmark.uta_local_tps | int }} | -| cdot REST (one request per transcript) | ~{{ benchmark.cdot_rest_tps | int }} | -| cdot REST (after one batch `prefetch()`) | ~{{ benchmark.cdot_rest_prefetch_tps | int }} | -| cdot local JSON | {{ benchmark.cdot_local_min_tps | int }}–{{ benchmark.cdot_local_max_tps | int }} | +| UTA: public remote database | {{ benchmark.uta_remote_tps | dp(2) }} ({{ benchmark.uta_remote_tps_q1 | dp(2) }}–{{ benchmark.uta_remote_tps_q3 | dp(2) }}) | +| UTA: local PostgreSQL | {{ benchmark.uta_local_tps | int }} ({{ benchmark.uta_local_tps_q1 | int }}–{{ benchmark.uta_local_tps_q3 | int }}) | +| cdot REST (one request per transcript) | {{ benchmark.cdot_rest_tps | int }} ({{ benchmark.cdot_rest_tps_q1 | int }}–{{ benchmark.cdot_rest_tps_q3 | int }}) | +| cdot REST (after one batch `prefetch()`) | {{ benchmark.cdot_rest_prefetch_tps | int }} ({{ benchmark.cdot_rest_prefetch_tps_q1 | int }}–{{ benchmark.cdot_rest_prefetch_tps_q3 | int }}) | +| cdot local JSON | {{ benchmark.cdot_local_tps | int }} ({{ benchmark.cdot_local_tps_q1 | int }}–{{ benchmark.cdot_local_tps_q3 | int }}) | A GRCh38 RefSeq JSON file loads in ~{{ benchmark.grch38_load_time_s | dp(0) }} s and then -resolves at {{ benchmark.cdot_local_min_tps | int }}–{{ benchmark.cdot_local_max_tps | -int }} HGVS/s through the biocommons engine. The REST provider serves -~{{ benchmark.cdot_rest_tps | int }} HGVS/s when each transcript version is fetched in its -own request; batching those lookups into one `prefetch()` request amortises the -per-request network and processing overhead across the whole set and warms the transcript -cache. Later lookups are then in-memory hits, so REST throughput rises to match local. A -locally loaded UTA reached only ~{{ benchmark.uta_local_tps | int }} HGVS/s, and the -public remote UTA database only ~{{ benchmark.uta_remote_tps | dp(1) }} HGVS/s. cdot's -local data layer is thus roughly 30× faster than a local UTA on the identical engine. -Once the data layer is local, the remaining bottleneck is sequence fetching rather than -transcript lookup. +resolves at {{ benchmark.cdot_local_tps | int }} HGVS/s (median) through the biocommons +engine. The REST provider serves {{ benchmark.cdot_rest_tps | int }} HGVS/s when each +transcript version is fetched in its own request; batching those lookups into one +`prefetch()` request (all transcripts for the set, under a second, untimed) warms the +transcript cache, after which lookups are in-memory hits and REST throughput +({{ benchmark.cdot_rest_prefetch_tps | int }} HGVS/s) is equivalent to local JSON, the +two differing by under 1% across repeats: with the transcript data in process memory, +both configurations are bounded by the shared engine and sequence layer, not by the +transcript backend. A locally loaded UTA on the +identical engine reached {{ benchmark.uta_local_tps | int }} HGVS/s, about a quarter of +local-JSON throughput, because each transcript lookup is a set of SQL queries rather +than a dict hit. The public remote UTA database, at +{{ benchmark.uta_remote_tps | dp(2) }} HGVS/s, is nearly four orders of magnitude slower +than any local configuration: every lookup pays wide-area round trips to a shared +server. Note that the absolute numbers are sensitive to the shared sequence layer: with +SeqRepo's file-descriptor cache disabled (its default), every configuration is bounded +by sequence fetching at well under 200 HGVS/s and the differences between transcript +backends are masked (Methods). 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). The REST provider -matched this over the network: after one batch cache-warming pass (21,277 distinct -transcripts warmed in ~6 s) it resolved the same set in ~83 minutes (731 HGVS/s), -effectively matching local JSON: once the cache is warm both providers resolve from an -in-memory dict with no further network I/O, so throughput is bounded by the shared -sequence layer. The marginal difference is within run-to-run variance, plausibly because -the warmed REST cache holds only the few thousand transcripts the set actually touches -whereas local JSON holds the entire dataset in memory. The same exhaustive pass is -impractical against the public remote UTA database (extrapolated at hundreds of days from -its ~0.1 HGVS/s). (`paper/scripts/build_clinvar_pairs.py` builds the pair set by joining -ClinVar's variant_summary with the ClinVar VCF.) +genomic coordinate, 98.8% matched the ClinVar genomic HGVS exactly), and the REST +provider, after one batch cache-warming pass (21,277 distinct transcripts warmed in +~6 s), resolved the same set in ~83 minutes (731 HGVS/s). Those two at-scale figures +come from separate one-off runs whose sequence-layer cache conditions differed, so +their small difference is not evidence of REST outrunning local JSON; the controlled +same-set comparison in Table 1, which holds those conditions fixed, shows the two +configurations are equivalent once the REST cache is warm. The same exhaustive pass is +impractical against the public remote UTA database (extrapolated at close to a year +from its ~{{ benchmark.uta_remote_tps | dp(2) }} HGVS/s). (`paper/scripts/build_clinvar_pairs.py` +builds the pair set by joining ClinVar's variant_summary with the ClinVar VCF.) ## R4: String cleaning recovers malformed real-world HGVS -The main test of cleaning is its effect on a real production query stream; a -reproducible injection benchmark provides a supporting safety check. +An independent large-scale production stream shows how common broken input is: +roughly 4 in 10 of the unique HGVS descriptions submitted to Mutalyzer carried a +syntactic or semantic error, and its checker automatically repaired only +~{{ literature.mutalyzer_autocorrect_pct | dp(0) }}% of unique descriptions overall +[@Lefter2021]. The main test of cleaning is its effect on a real production query +stream; a reproducible injection benchmark provides a supporting safety check. **[Tier 2].** The corpus is N = 32,752 real queries typed into the HGVS search box of production clinical and @@ -146,31 +217,50 @@ to a variant or its classification, so the strings are whatever a clinician or c happened to paste or type, not curated HGVS. They arrive carrying the damage of their route to the box: stray whitespace and non-printable characters from copying out of Word documents and report PDFs and pasting between systems, lost casing, transposed -punctuation, and trailing protein annotations. Run over this corpus, `clean_hgvs()` raised -the fraction parseable by biocommons/hgvs from 91.5% as-submitted to 96.6% after -cleaning: a +5.1% absolute gain (1,678 additional strings rescued) with zero -regressions (no already-valid string was broken). This measures what cleaning recovers -from messy, human-entered input rather than from synthetic errors. The rescues break down by fix type as shown in -Table 2; they are dominated by whitespace removal and base re-casing, followed by -protein-suffix stripping and gene/transcript-wrapper repair. +punctuation, and trailing protein annotations. Run over this corpus, the cleaning +pipeline (`clean_hgvs()` plus the data-provider-verified accession-prefix restoration +below) raised the fraction parseable by biocommons/hgvs from 91.5% as-submitted to +96.7% after cleaning: a +5.3% absolute gain (1,721 additional strings rescued) with zero +regressions (no already-valid string was broken). That gain is 5.3 of the 8.5 +percentage points that failed as-submitted, about 62% of the failures. This measures +what cleaning recovers from messy, human-entered input rather than from synthetic +errors. The rescues break down by fix type as shown in +Table 2; they are dominated by whitespace removal and structural-punctuation repair, +followed by gene/transcript-wrapper repair and structure reconstruction. One category is +not a pure string operation: a bare-number accession whose RefSeq prefix was dropped +entirely (`000059.4:c.68del`) is repaired by generating the candidate accessions the +kind letter allows (`c.` implies `NM_`/`XM_`, `n.` implies `NR_`/`XR_`), checking each +against the loaded transcript data, and restoring the prefix only when exactly one +candidate exists there. This repair is impossible for a purely string-level checker: +having the transcript data locally is what makes it safe. + +The public submitted-string ClinVar corpus (R2) provides the reproducible complement to +this Tier 2 result, and locates the two failure axes: formal database submissions are +largely well-formed, so cleaning fires rarely there (one rescue in the 3,000-pair +sample) and failures are dominated by transcript-version age, whereas interactive +human-typed input carries the formatting damage that cleaning repairs. cdot addresses +the first with historical transcript depth and the second with `clean_hgvs()`. **Table 2. Fixes applied across the production corpus (N = 32,752).** Each row is a -`clean_hgvs()` fix category, with the number of rescued queries in which it fired and that -number as a share of the 1,678 rescued queries. Categories overlap (a single query may +cleaning fix category, with the number of rescued queries in which it fired and that +number as a share of the 1,721 rescued queries. Categories overlap (a single query may need several fixes), so the counts sum to more than the total. *(Tier 2; counts are -frozen constants from a deterministic run of `clean_hgvs()` over the production corpus.)* +frozen constants from a deterministic run of the cleaning pipeline, `clean_hgvs()` plus +the provider-verified accession-prefix restoration, over the production corpus.)* | Fix category | Example (→ repaired) | Rescued queries | % of rescued | |---|---|---|---| -| Whitespace / non-printable removal | `NM_000059.4: c.1A>G` → `NM_000059.4:c.1A>G` | 940 | 56.0% | -| Base re-casing | `NM_000059.4:c.1delg` → `…delG` | 553 | 33.0% | -| Structural-punctuation repair | `NM_000059..4:c.1A>G` → `NM_000059.4:c.1A>G` | 310 | 18.5% | -| Gene/transcript-wrapper repair | `BRCA2(NM_000059.4):c.1A>G` → `NM_000059.4(BRCA2):c.1A>G` | 209 | 12.5% | -| Protein-suffix stripping | `NM_000059.4:c.1A>G p.(Met1?)` → `NM_000059.4:c.1A>G` | 130 | 7.7% | -| Genomic-ref-in-parens removal | `NM_000059.4(NC_000013.11):c.68del` → `NM_000059.4:c.68del` | 57 | 3.4% | -| Other (del/dup count, mutation-type case, …) | `NM_000059.4:c.1_2del2` → `…del` | 21 | 1.3% | -| Prefix / kind restoration | `NM_000059.4:1A>G` → `NM_000059.4:c.1A>G` | 11 | 0.7% | -| **Total unique queries rescued** | | **1,678** | **100%** | +| Whitespace / non-printable removal | `NM_000059.4: c.1A>G` → `NM_000059.4:c.1A>G` | 664 | 38.6% | +| Structural-punctuation repair | `NM_000059.4;c.1A>G` → `NM_000059.4:c.1A>G` | 640 | 37.2% | +| Gene/transcript-wrapper repair | `BRCA2(NM_000059.4):c.1A>G` → `NM_000059.4(BRCA2):c.1A>G` | 260 | 15.1% | +| Structure reconstruction | `NM_000059.4c.1A>G` → `NM_000059.4:c.1A>G` | 248 | 14.4% | +| Protein-suffix stripping | `NM_000059.4:c.1A>G p.(Met1?)` → `NM_000059.4:c.1A>G` | 130 | 7.6% | +| Genomic-ref-in-parens removal | `NM_000059.4(NC_000013.11):c.68del` → `NM_000059.4:c.68del` | 57 | 3.3% | +| Base re-casing | `NM_000059.4:c.1delg` → `…delG` | 42 | 2.4% | +| Other (del/dup count, mutation-type case, …) | `NM_000059.4:c.1_2del2` → `…del` | 11 | 0.6% | +| Prefix / kind restoration | `NM_000059.4:1A>G` → `NM_000059.4:c.1A>G` | 4 | 0.2% | +| Accession prefix restoration (provider-verified) | `000059.4:c.68del` → `NM_000059.4:c.68del` | 1 | 0.1% | +| **Total unique queries rescued** | | **1,721** | **100%** | As a reproducible control, `paper/scripts/inject_and_clean.py` injects each `clean_hgvs()` fix category into a seeded @@ -184,35 +274,35 @@ depends. ### Residual errors: the ceiling of cleaning *(Table S6)* -**[Tier 2].** The 3.4% of the production corpus (1,118 queries; 860 +**[Tier 2].** The 3.3% of the production corpus (1,075 queries; 826 unique strings) that still fail to parse after cleaning define the ceiling of what pure string repair can achieve. Each residual string was assigned to one single-label error class under a fixed decision-tree taxonomy. Of the eight classes, the seven repair-relevant ones are shown below with synthesised examples; the eighth was non-HGVS -input (81 queries, 7.2%: pasted URLs, report templates, or prose) and is excluded from +input (81 queries, 7.5%: pasted URLs, report templates, or prose) and is excluded from the table, as there is nothing in it for cleaning to repair. -**Residual error classes after cleaning** (counts and % of the 1,118 residual +**Residual error classes after cleaning** (counts and % of the 1,075 residual queries; examples synthesised from public BRCA2 `NM_000059.4`). *(Tier 2; frozen constants from a deterministic run over the production corpus.)* | Class | Queries | What it is (*example*) | |---|---|---| -| Truncated | 284 (25.4%) | cut off before a complete variant: `NM_000059.4:c.68_69` (range, no edit) | -| No reference | 277 (24.8%) | a bare variant body, no transcript/gene/accession: `c.68_69delAG` | -| Bad accession | 167 (14.9%) | missing prefix, or misplaced/truncated version: `000059.4:c.68del` (`NM_` prefix dropped) | -| Edit syntax | 113 (10.1%) | malformed or non-standard edit operation: `NM_000059.4:c.68AG>T` (multi-base reference in a substitution) | -| Trailing / concatenated | 85 (7.6%) | extra characters after a complete variant, or several run together: `NM_000059.4:c.68delAG;c.70A>G` | -| Grammar gap | 81 (7.2%) | legitimate HGVS the biocommons grammar rejects: `NM_000059.4:c.(67+1_68-1)_(70+1_71-1)del` (uncertain-range deletion) | -| Insertion (length only) | 30 (2.7%) | an insertion given as a base count, not a sequence: `NM_000059.4:c.68_69ins5` (position and length recoverable; inserted bases not) | - -The residual falls into three groups. Just over half (~50%: Truncated + No reference) is +| Truncated | 284 (26.4%) | cut off before a complete variant: `NM_000059.4:c.68_69` (range, no edit) | +| No reference | 277 (25.8%) | a bare variant body, no transcript/gene/accession: `c.68_69delAG` | +| Bad accession | 124 (11.5%) | misplaced or truncated version, or a missing prefix with no unique data match: `NM_000059/4:c.68del` (slash in place of the version dot) | +| Edit syntax | 113 (10.5%) | malformed or non-standard edit operation: `NM_000059.4:c.68AG>T` (multi-base reference in a substitution) | +| Trailing / concatenated | 85 (7.9%) | extra characters after a complete variant, or several run together: `NM_000059.4:c.68delAG;c.70A>G` | +| Grammar gap | 81 (7.5%) | legitimate HGVS the biocommons grammar rejects: `NM_000059.4:c.(67+1_68-1)_(70+1_71-1)del` (uncertain-range deletion) | +| Insertion (length only) | 30 (2.8%) | an insertion given as a base count, not a sequence: `NM_000059.4:c.68_69ins5` (position and length recoverable; inserted bases not) | + +The residual falls into three groups. Just over half (~52%: Truncated + No reference) is incomplete or reference-less user input: information the user never supplied, which no -string-level repair can invent; the integer-length insertions (2.7%) belong with these, +string-level repair can invent; the integer-length insertions (2.8%) belong with these, since the inserted bases were never given (only the position and length are recoverable, so -the variant is not properly resolvable). About 33% (Bad accession + Edit syntax + Trailing / +the variant is not properly resolvable). About 30% (Bad accession + Edit syntax + Trailing / concatenated) is in principle fixable and marks the frontier for future cleaning rules. -The remaining ~14% splits into two equal classes of 81 queries (7.2%) each: a grammar gap +The remaining ~15% splits into two equal classes of 81 queries (7.5%) each: a grammar gap (valid HGVS the biocommons grammar rejects rather than the input) and the non-HGVS input excluded above (strings that should not be parsed at all). Most of what remains is therefore either out of scope for any repair or a downstream grammar limitation. @@ -220,8 +310,12 @@ therefore either out of scope for any repair or a downstream grammar limitation. *Method and limitation:* classification was performed by a large language model (Claude Opus 4, Anthropic; 2026-06-17) applying the shared decision tree to each unique string, single-label and single-rater; no second-rater adjudication was done, so no inter-rater -agreement (κ) is reported. The taxonomy is version `v1`. Synthesised examples (from public -NM_000059.4 / NM_001754.5) illustrate each class; no corpus string is reproduced. +agreement (κ) is reported. The taxonomy is version `v1`. Counts were refreshed on +2026-08-17 after the accession-repair cleaning rules landed: the 43 residual queries the +new rules rescue all sat in the Bad accession class (delta re-rated by Claude Fable 5, +Anthropic), which drops from 167 (14.9% of the previous 1,118-query residual) to 124; +the other classes are unchanged. Synthesised examples (from public NM_000059.4 / +NM_001754.5) illustrate each class; no corpus string is reproduced. ## R5: Transcript version fallback and safe substitution @@ -253,7 +347,16 @@ dangerous case, a *partial* bump that mis-places some variants but not others, i ({{ version_stability.refseq_partial_drift_pct | dp(1) }}% of RefSeq bumps, {{ version_stability.ensembl_partial_drift_pct | dp(1) }}% Ensembl); when a coordinate does move it is almost always the whole CDS, driven by a re-annotation of the coding -region. Drift is also highly concentrated: only +region. Within that small partial-drift tail the risk is strongly positional: a partial +drift keeps a 5' prefix intact up to its first alignment change and moves what lies +downstream, so preservation falls monotonically along the CDS, from +{{ positional_drift.refseq_partial_decile1_pct | dp(0) }}% of coding bases in the 5'-most +decile to {{ positional_drift.refseq_partial_decile10_pct | dp(0) }}% in the 3'-most for +RefSeq ({{ positional_drift.ensembl_partial_decile1_pct | dp(0) }}% to +{{ positional_drift.ensembl_partial_decile10_pct | dp(0) }}% for Ensembl), while the +whole-CDS relocations that dominate drift are position independent (Supplementary Figure +S1). As a rule of thumb, when the safety of a substitution cannot be checked, a 5' coding +variant is a safer bet than a 3' one. Drift is also highly concentrated: only {{ version_stability.refseq_accessions_drift_pct | dp(1) }}% of RefSeq accessions ever drift, and within that set the moves cluster in a small minority of accessions. This concentration has a practical consequence: the risky transcripts form a short, diff --git a/paper/scripts/bin_residual_positions.py b/paper/scripts/bin_residual_positions.py new file mode 100644 index 0000000..ed75866 --- /dev/null +++ b/paper/scripts/bin_residual_positions.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +"""Bin the ClinVar VCF-pass projections by relative CDS position (R5 positional drift). + +Companion to the positional-drift deciles in compute_version_stability.py: that +script asks where along the CDS a *version bump* moves coordinates; this one asks +where along the CDS the residual `incorrect` ClinVar projections sit. If +projection errors concentrated toward the 3' end, the incorrect rows' relative +CDS positions would skew high against the correct rows' baseline distribution. + +Consumes the per-variant table from resolve_clinvar_pass.py (VCF mode; the run +documented in claude/clinvar_diff_coordinates.md) plus the same cdot release +JSON.gz the pass used (only start/stop codons are kept, then the JSON is freed). +Every `correct` and `incorrect` row's cited c. position is anchored to a relative +CDS position (intronic offsets anchor to their exon boundary; ranges to their +5'-most end) and binned into deciles (1 = 5'-most tenth, 10 = 3'-most); 5'UTR +(c.-N) and 3'UTR (c.*N) citations are counted separately, as are n. transcripts. + +Writes a single-row facts CSV: per-decile *shares* (percent of that bucket's +binned coding rows) for the incorrect and correct buckets, plus the UTR / +non-coding / unbinnable counts. The committed snapshot lives at +paper/empirical_results/clinvar_residual_positions.csv (no Snakefile rule, like +the other clinvar_vcf facts: the per-variant table is a dedicated multi-hour run). + +Per CLAUDE.md: never run against full datasets in development; verify on a small +table from tests/test_data/clinvar_hgvs/ first. + +Usage: + python paper/scripts/bin_residual_positions.py \ + output/clinvar_pass/refseq_full_vcf.csv \ + --json cdot-0.2.33.refseq.GRCh38.json.gz \ + --out output/facts/clinvar_residual_positions.csv +""" +import argparse +import csv +import gzip +import json +import re +from collections import Counter +from pathlib import Path + +N_BINS = 10 + +# First cited position after the kind letter: "c.100A>G", "c.-49_12del", +# "c.*103del", "c.4072-1234G>T" -> the leading (possibly -/*) anchor number. +_POS = re.compile(r":[cn]\.(\*?)(-?\d+)") + + +def load_cds_lengths(path): + """Return {accession.version: cds_len} for coding transcripts in a release.""" + with gzip.open(path, "rt") as fh: + transcripts = json.load(fh)["transcripts"] + out = {} + for key, t in transcripts.items(): + sc, ec = t.get("start_codon"), t.get("stop_codon") + if sc is not None and ec is not None and ec > sc: + out[key] = ec - sc + return out + + +def bin_row(c_hgvs, cds_lengths): + """Classify one row: decile 1..N_BINS, '5putr', '3putr', 'noncoding_n', + 'no_cds_len' or 'unparsed'.""" + m = _POS.search(c_hgvs) + if m is None: + return "unparsed" + if ":n." in c_hgvs: + return "noncoding_n" + star, pos = m.group(1), int(m.group(2)) + if star: + return "3putr" + if pos < 0: + return "5putr" + key = c_hgvs.split(":", 1)[0] + cds_len = cds_lengths.get(key) + if not cds_len: + return "no_cds_len" + pos = min(max(pos, 1), cds_len) + return min((pos - 1) * N_BINS // cds_len, N_BINS - 1) + 1 + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("pass_csv", help="per-variant table from resolve_clinvar_pass.py") + ap.add_argument("--json", required=True, help="cdot release JSON.gz the pass used") + ap.add_argument("--out", default="output/facts/clinvar_residual_positions.csv") + args = ap.parse_args() + + print(f"loading CDS lengths from {args.json} ...", flush=True) + cds_lengths = load_cds_lengths(args.json) + print(f" {len(cds_lengths)} coding transcript versions") + + bins = {"correct": Counter(), "incorrect": Counter()} + with open(args.pass_csv, newline="") as fh: + for row in csv.DictReader(fh): + counter = bins.get(row["bucket"]) + if counter is not None: + counter[bin_row(row["c_hgvs"], cds_lengths)] += 1 + + facts = {} + for bucket, counter in bins.items(): + binned = sum(counter[i] for i in range(1, N_BINS + 1)) + facts[f"{bucket}_binned"] = binned + for i in range(1, N_BINS + 1): + facts[f"{bucket}_decile{i}_pct"] = ( + round(100 * counter[i] / binned, 1) if binned else 0.0) + for extra in ("5putr", "3putr", "noncoding_n", "no_cds_len", "unparsed"): + facts[f"{bucket}_{extra}"] = counter[extra] + + 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(facts)) + w.writeheader() + w.writerow(facts) + print(f"Written: {out}") + for k, v in facts.items(): + print(f" {k}: {v}") + + +if __name__ == "__main__": + main() diff --git a/paper/scripts/build_clinvar_submitted_pairs.py b/paper/scripts/build_clinvar_submitted_pairs.py new file mode 100644 index 0000000..bda5875 --- /dev/null +++ b/paper/scripts/build_clinvar_submitted_pairs.py @@ -0,0 +1,241 @@ +#!/usr/bin/env python3 +""" +Build the SUBMITTED-HGVS ClinVar benchmark corpus: SCV-submitted transcript HGVS +strings (as the submitting laboratory wrote them) joined to the variant's genomic +ground truth via AlleleID. + +Why this exists: build_clinvar_pairs.py takes c.HGVS from the ``Name`` column of +variant_summary.txt.gz, which is ClinVar's own recomputed preferred-transcript name +and is always at the CURRENT transcript version. Benchmarks on that corpus measure +the easy case. The per-SCV HGVS attributes in the VCV XML preserve what each lab +actually submitted, including transcript versions that were current when the +variant was classified but have since been superseded, so this corpus can exercise +historical transcript depth. + +Sources (two input modes): + + * ``--xml ClinVarVCVRelease_YYYY-MM.xml.gz`` - stream the full VCV XML release + (ncbi.nlm.nih.gov/clinvar; large, tens of GB uncompressed) with iterparse. + Per VariationArchive: AlleleID from ClassifiedRecord/SimpleAllele (or + InterpretedRecord in pre-2025 releases); per ClinicalAssertion: the first + transcript c./n. expression among its ``AttributeSet/Attribute[@Type="HGVS"]`` + values, kept verbatim. Fully self-contained but slow (hours). + * ``--scv-csv-dir DIR`` - consume per-SCV CSVs previously extracted from that + XML (columns ``hgvs`` = first HGVS attribute value, ``allele_id``). Fast path + when an extraction already exists. The two modes differ only in that the CSV + keeps the first HGVS value rather than the first *transcript* value; on a + 50,000-row sample the two never disagreed (no SCV had a transcript expression + preceded by a non-transcript one). + +Ground truth: the ClinVar VCF (clinvar.GRCh38.vcf.gz), joined on INFO/ALLELEID, +exactly as in build_clinvar_pairs.py: CHROM/POS/REF/ALT plus the CLNHGVS g.HGVS +kept for reference. Scoring against the VCF coordinate (not the g.HGVS string) is +essential here: a submitted string may legitimately spell an indel differently +from ClinVar's normalised form. + +Inclusion rule (documented, deliberately loose so submitted messiness survives): +a string is kept if it starts with a RefSeq/Ensembl transcript accession +(NM_/NR_/XM_/XR_/ENST + digits, version optional) and contains ``:c.`` or +``:n.``. Strings are kept verbatim apart from stripping outer whitespace and +internal tab/newline characters (TSV safety). + +Dedup rule (documented): SCV rows are collapsed to unique +(AlleleID, submitted string) pairs; ``scv_count`` records how many SCVs +contributed each pair. Rationale: submitted HGVS is not deduplicated across SCVs +in ClinVar, and repeated identical submissions measure submitter volume, not +string diversity; distinct strings for the same variant (e.g. two labs citing +different transcript versions) are all kept. + +Sampling rule (documented): ``--sample N --seed S`` draws a uniform random sample +of N pairs with ``random.Random(S).sample`` and writes them in original file +order; the committed test sample and the benchmark sample both use seed 42. + +Output TSV: header ``chrom pos ref alt g_hgvs c_hgvs scv_count`` - the VCF-format +pairs layout consumed by resolve_clinvar_pass.py (extra columns are ignored by its +reader). The full corpus is large and data-derived so it lives under a data dir +and is NOT committed; a 500-pair sample is committed to +tests/test_data/clinvar_hgvs/. + +Usage: + # fast path: consume an existing per-SCV CSV extraction + python paper/scripts/build_clinvar_submitted_pairs.py \ + --scv-csv-dir extracted/ClinVarVCVRelease_2026-06 \ + clinvar.GRCh38.vcf.gz clinvar_submitted_pairs.GRCh38.tsv + + # self-contained: stream the VCV XML release (slow; --limit for a smoke test) + python paper/scripts/build_clinvar_submitted_pairs.py \ + --xml ClinVarVCVRelease_2026-06.xml.gz \ + clinvar.GRCh38.vcf.gz clinvar_submitted_pairs.GRCh38.tsv + + # sample an existing corpus (no rebuild) + python paper/scripts/build_clinvar_submitted_pairs.py \ + --from-pairs clinvar_submitted_pairs.GRCh38.tsv \ + --sample 500 --seed 42 --sample-out clinvar_submitted_500.tsv +""" +import argparse +import csv +import glob +import gzip +import random +import re +import sys +from pathlib import Path + +# Reuse the VCF ground-truth loader so the two corpora share one join convention. +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from build_clinvar_pairs import load_vcf_g, source_of # noqa: E402 + +# Loose inclusion filter: transcript accession at the start (version optional, +# wrappers/messiness allowed) and a coding/non-coding kind somewhere after it. +_TX_START = re.compile(r"^(?:NM_|NR_|XM_|XR_|ENST)\d+") +_HAS_KIND = re.compile(r":[cn]\.") + + +def is_submitted_tx_hgvs(s): + return bool(_TX_START.match(s)) and bool(_HAS_KIND.search(s)) + + +def sanitize(s): + """Strip outer whitespace and TSV-breaking characters, keep the rest verbatim.""" + return s.strip().replace("\t", " ").replace("\n", " ").replace("\r", " ") + + +def iter_scv_csv(csv_dir): + """Yield (allele_id, submitted_hgvs) from an extract_xml_to_csv.py CSV dir. + + The extraction's ``hgvs`` column holds the first HGVS attribute value of the + SCV (which may be genomic or protein; the transcript filter is applied here).""" + files = sorted(glob.glob(str(Path(csv_dir) / "*.csv"))) + if not files: + sys.exit(f"no CSVs found in {csv_dir}") + for fn in files: + with open(fn, newline="") as fh: + for row in csv.DictReader(fh): + h = sanitize(row.get("hgvs") or "") + allele_id = row.get("allele_id") or "" + if h and allele_id: + yield allele_id, h + + +def iter_scv_xml(xml_path, limit=0): + """Yield (allele_id, submitted_hgvs) by streaming a ClinVar VCV XML release. + + One VariationArchive at a time (iterparse + clear), so memory stays flat. + Yields the FIRST transcript c./n. expression per ClinicalAssertion.""" + try: + from lxml import etree + except ImportError: # pragma: no cover - stdlib fallback + 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 _event, elem in etree.iterparse(fh, events=("end",), tag="VariationArchive"): + n_va += 1 + record = elem.find("ClassifiedRecord") + if record is None: # pre-Aug-2025 releases + record = elem.find("InterpretedRecord") + if record is not None: + sa = record.find("SimpleAllele") + allele_id = sa.get("AlleleID") if sa is not None else None + if allele_id: + for ca in record.iter("ClinicalAssertion"): + 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: + h = sanitize(attr.text) + if is_submitted_tx_hgvs(h): + yield allele_id, h + break + elem.clear() + while elem.getprevious() is not None: # lxml only; frees siblings + del elem.getparent()[0] + if limit and n_va >= limit: + return + + +def write_sample(pairs_path, n, seed, out_path): + """Uniform random sample of n data lines, written in original order.""" + with open(pairs_path) as fh: + header = fh.readline() + lines = fh.readlines() + if n > len(lines): + sys.exit(f"--sample {n} > corpus size {len(lines)}") + idx = sorted(random.Random(seed).sample(range(len(lines)), n)) + with open(out_path, "w") as out: + out.write(header) + for i in idx: + out.write(lines[i]) + print(f" sampled {n}/{len(lines):,} (seed {seed}) -> {out_path}", file=sys.stderr) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + src = ap.add_mutually_exclusive_group(required=True) + src.add_argument("--scv-csv-dir", help="dir of per-SCV CSVs (extract_xml_to_csv.py output)") + src.add_argument("--xml", help="ClinVar VCV XML release (.xml.gz); slow streaming parse") + src.add_argument("--from-pairs", help="existing corpus TSV; skip the build and only sample") + ap.add_argument("vcf", nargs="?", help="ClinVar VCF (ground truth, joined on ALLELEID)") + ap.add_argument("out", nargs="?", help="output corpus TSV") + ap.add_argument("--limit", type=int, default=0, help="(--xml) stop after N VariationArchives (smoke test)") + ap.add_argument("--sample", type=int, default=0, help="also write a random sample of N pairs") + ap.add_argument("--seed", type=int, default=42) + ap.add_argument("--sample-out", help="path for the sample TSV") + args = ap.parse_args() + + if args.from_pairs: + if not (args.sample and args.sample_out): + sys.exit("--from-pairs needs --sample and --sample-out") + write_sample(args.from_pairs, args.sample, args.seed, args.sample_out) + return + if not (args.vcf and args.out): + sys.exit("vcf and out are required when building") + + print("reading VCF ground truth...", file=sys.stderr) + g_by_allele = load_vcf_g(args.vcf) + print(f" {len(g_by_allele):,} alleles with genomic coordinate", file=sys.stderr) + + scv_iter = (iter_scv_xml(args.xml, args.limit) if args.xml + else iter_scv_csv(args.scv_csv_dir)) + + n_scv = n_tx = n_joined = 0 + comp = {"refseq": 0, "ensembl": 0} + pairs = {} # (allele_id, c_hgvs) -> scv_count + for allele_id, h in scv_iter: + n_scv += 1 + if not is_submitted_tx_hgvs(h): + continue + n_tx += 1 + if allele_id not in g_by_allele: + continue + n_joined += 1 + key = (allele_id, h) + if key not in pairs: + comp[source_of(h.split(":", 1)[0])] += 1 + pairs[key] = pairs.get(key, 0) + 1 + + with open(args.out, "w") as out: + out.write("chrom\tpos\tref\talt\tg_hgvs\tc_hgvs\tscv_count\n") + for (allele_id, c_hgvs), count in pairs.items(): + chrom, pos, ref, alt, g_hgvs = g_by_allele[allele_id] + out.write(f"{chrom}\t{pos}\t{ref}\t{alt}\t{g_hgvs}\t{c_hgvs}\t{count}\n") + + n_pairs = len(pairs) + print(f"scanned {n_scv:,} SCV HGVS values -> {n_tx:,} transcript c./n. strings " + f"-> {n_joined:,} with VCF ground truth -> {n_pairs:,} unique " + f"(AlleleID, string) pairs", file=sys.stderr) + if n_pairs: + print(f" source mix: refseq {comp['refseq']:,} ({100*comp['refseq']/n_pairs:.2f}%) " + f"ensembl {comp['ensembl']:,} ({100*comp['ensembl']/n_pairs:.2f}%)", file=sys.stderr) + print(f"Written: {args.out}", file=sys.stderr) + + if args.sample: + if not args.sample_out: + sys.exit("--sample needs --sample-out") + write_sample(args.out, args.sample, args.seed, args.sample_out) + + +if __name__ == "__main__": + main() diff --git a/paper/scripts/compute_benchmark.py b/paper/scripts/compute_benchmark.py index 965c2cb..97296b9 100644 --- a/paper/scripts/compute_benchmark.py +++ b/paper/scripts/compute_benchmark.py @@ -1,123 +1,212 @@ #!/usr/bin/env python3 -"""Benchmark cdot local JSON throughput and load time vs UTA remote access. +"""Throughput benchmark for Table 1: HGVS resolution speed by transcript backend. + +Every configuration resolves the SAME fixed set of committed ClinVar (g.HGVS, c.HGVS) +pairs (tests/test_data/clinvar_hgvs/clinvar_hgvs_500.tsv) through the identical +biocommons/hgvs engine, with the sequence layer held constant: one shared local SeqRepo +SeqFetcher instance (file-descriptor caching enabled) is passed to every provider, so +the only thing that varies between configurations is the transcript-data layer. + +Methodology (Results R3 / Methods): + * N_REPEATS timed passes per configuration over the identical pair set. + * Timed: the resolution loop only (parse + c_to_g per string, fresh AssemblyMapper + per pass). Untimed: provider construction/load, REST prefetch(), and one warm-up + pass per configuration (OS page cache, SeqRepo file handles). + * The engine-side memo (the 100-entry LRU biocommons hgvs wraps around provider + methods) is cleared before every timed pass - see clear_engine_caches(). + * Reported: median and IQR (Q1-Q3) of per-pass throughput. + * Public remote UTA is ~0.1 HGVS/s, so a full-size pass is impractical; it is + measured on the first UTA_REMOTE_N pairs of the same set (documented in Methods). Writes output/facts/benchmark.csv. Usage: python paper/scripts/compute_benchmark.py \ --refseq-grch38 cdot.refseq.grch38.json.gz \ - [--uta-uri postgresql://uta_admin@localhost/uta/uta_20210129] + [--uta-uri postgresql://postgres@127.0.0.1:5433/uta/uta_20241220] \ + [--skip-uta-remote] -Requires: hgvs, cdot +Requires: hgvs, cdot, a local SeqRepo (HGVS_SEQREPO_DIR). """ import argparse +import os +import statistics +import sys import time from pathlib import Path -import pandas as pd - - -N_WARMUP = 10 -N_BENCHMARK = 500 +# The sequence layer must be identical (and fast) for every configuration: enable +# SeqRepo file-descriptor caching before biocommons.seqrepo is imported. Without it +# every fetch re-opens a bgzf file (~1ms) and the sequence layer, not the transcript +# layer, dominates every row. +os.environ.setdefault("SEQREPO_FD_CACHE_MAXSIZE", "128") +import pandas as pd -def benchmark_cdot_local(json_gz_path: str) -> tuple[float, float]: - """Return (load_time_s, transcripts_per_second) for local JSON provider.""" - from cdot.hgvs.dataproviders import JSONDataProvider - - t0 = time.perf_counter() - hdp = JSONDataProvider([json_gz_path]) - load_time_s = time.perf_counter() - t0 - - # Get a sample of transcript accessions to query - import gzip - import json - with gzip.open(json_gz_path, "rt") as fh: - data = json.load(fh) - accessions = list(data["transcripts"].keys())[:N_BENCHMARK + N_WARMUP] - - # Warm up - for acc in accessions[:N_WARMUP]: - try: - hdp.get_tx_info(acc) - except Exception: - pass - - # Benchmark - sample = accessions[N_WARMUP:N_WARMUP + N_BENCHMARK] +REPO_ROOT = Path(__file__).resolve().parents[2] +DEFAULT_PAIRS = REPO_ROOT / "tests/test_data/clinvar_hgvs/clinvar_hgvs_500.tsv" + +N_REPEATS = 5 +UTA_REMOTE_N = 25 # public remote UTA is ~0.1 HGVS/s; a 500-pair pass is impractical + + +def clear_engine_caches(hdp): + """biocommons hgvs wraps every data-provider method (get_tx_exons, get_seq, ...) + in a per-instance 100-entry LRU at Interface construction. Clear that engine-side + memo before each timed pass, so a pass measures the provider itself, including any + internal caching the provider ships (cdot's in-memory transcript dict and exon + cache, REST's prefetched transcript cache), rather than the engine memo. Without + this, a working set smaller than 100 transcripts (the remote-UTA subset) is served + entirely from the engine memo after the first pass and the backend is never hit.""" + for name in ("get_seq", "get_tx_exons", "get_tx_info", "get_tx_identity_info", + "get_tx_mapping_options", "get_pro_ac_for_tx_ac", "get_gene_info", + "get_acs_for_protein_seq", "get_tx_for_gene", "get_tx_for_region"): + cache_clear = getattr(getattr(hdp, name, None), "cache_clear", None) + if cache_clear is not None: + cache_clear() + + +def resolve_pass(hdp, hp, pairs, build="GRCh38"): + """One timed pass: parse + project every pair; fresh AssemblyMapper per pass. + Returns (throughput_hgvs_per_s, counts).""" + from hgvs.assemblymapper import AssemblyMapper + from benchmark_resolution import classify + + am = AssemblyMapper(hdp, assembly_name=build, alt_aln_method="splign", + replace_reference=False) + counts = {} t0 = time.perf_counter() - for acc in sample: - try: - hdp.get_tx_info(acc) - except Exception: - pass + for g_hgvs, c_hgvs in pairs: + bucket, _ = classify(am, hp, g_hgvs, c_hgvs) + counts[bucket] = counts.get(bucket, 0) + 1 elapsed = time.perf_counter() - t0 - tps = len(sample) / elapsed if elapsed > 0 else 0 - - return round(load_time_s, 2), round(tps) - - -def benchmark_uta_remote(uta_uri: str) -> float: - """Return transcripts_per_second for UTA remote access.""" - import hgvs.dataproviders.uta as uta + return len(pairs) / elapsed, counts + + +def run_config(name, make_provider, hp, pairs, repeats=N_REPEATS, + fresh_provider_per_repeat=False, warmup=True): + """Run one configuration: optional untimed warm-up pass, then `repeats` timed + passes. `make_provider` is called once for persistent-provider configs, or once + per repeat when the configuration is defined by a cold provider cache (REST).""" + print(f"== {name} ==") + provider = None + if not fresh_provider_per_repeat: + provider = make_provider() + if warmup: + warm_provider = provider if provider is not None else make_provider() + tps, counts = resolve_pass(warm_provider, hp, pairs) + print(f" warm-up (untimed): {tps:7.1f} HGVS/s {counts}") + tps_values = [] + for i in range(repeats): + hdp = provider if provider is not None else make_provider() + clear_engine_caches(hdp) + tps, counts = resolve_pass(hdp, hp, pairs) + tps_values.append(tps) + print(f" repeat {i + 1}/{repeats}: {tps:7.1f} HGVS/s {counts}") + med = statistics.median(tps_values) + q1, _q2, q3 = statistics.quantiles(tps_values, n=4, method="inclusive") + print(f" median {med:.1f} IQR {q1:.1f}-{q3:.1f}") + return {"median": med, "q1": q1, "q3": q3, "repeats": tps_values} - hdp = uta.connect(db_url=uta_uri) - # Get sample accessions - sample_accessions = [ - "NM_000492.3", "NM_007294.3", "NM_000059.3", "NM_000546.5", "NM_001354689.1", - "NM_000249.3", "NM_001126112.2", "NM_005343.3", "NM_004333.5", "NM_000179.2", - ] +def main(): + parser = argparse.ArgumentParser(description="Benchmark HGVS resolution throughput by backend.") + parser.add_argument("--refseq-grch38", required=True, help="cdot GRCh38 RefSeq JSON.gz") + parser.add_argument("--uta-uri", default=None, help="local UTA PostgreSQL URI") + parser.add_argument("--pairs", default=str(DEFAULT_PAIRS), help="fixed (g.HGVS, c.HGVS) pair set") + parser.add_argument("--repeats", type=int, default=N_REPEATS) + parser.add_argument("--skip-uta-remote", action="store_true", + help="skip the (slow, ~25 min) public remote UTA rows") + parser.add_argument("--skip-rest", action="store_true", help="skip the REST configurations") + args = parser.parse_args() - # Warm up - for acc in sample_accessions[:3]: - try: - hdp.get_tx_info(acc, "GRCh38") - except Exception: - pass + sys.path.insert(0, str(Path(__file__).resolve().parent)) - t0 = time.perf_counter() - n = 0 - for acc in sample_accessions * 5: - try: - hdp.get_tx_info(acc, "GRCh38") - n += 1 - except Exception: - pass - elapsed = time.perf_counter() - t0 - return round(n / elapsed) if elapsed > 0 else 0 + import hgvs.parser + from hgvs.dataproviders.seqfetcher import SeqFetcher + from cdot.hgvs.dataproviders import JSONDataProvider, RESTDataProvider + from benchmark_resolution import load_pairs + if not os.environ.get("HGVS_SEQREPO_DIR"): + sys.exit("Set HGVS_SEQREPO_DIR to a local SeqRepo: the sequence layer must be " + "local and identical for every configuration.") -def main(): - parser = argparse.ArgumentParser(description="Benchmark cdot vs UTA throughput.") - parser.add_argument("--refseq-grch38", required=True, help="cdot GRCh38 RefSeq JSON.gz") - parser.add_argument("--uta-uri", default=None, help="UTA PostgreSQL URI for remote benchmark") - args = parser.parse_args() + pairs = load_pairs(args.pairs) + hp = hgvs.parser.Parser() + # One shared SeqRepo seqfetcher for every provider: SeqRepo's fd cache is + # per-instance, so per-provider seqfetchers would give configurations + # differently-warmed sequence layers. + shared_seqfetcher = SeqFetcher() - print("Benchmarking cdot local JSON ...") - load_time_s, cdot_local_tps = benchmark_cdot_local(args.refseq_grch38) - print(f" Load time: {load_time_s}s") - print(f" Throughput: {cdot_local_tps:,} transcripts/s") + facts = {"n_pairs": len(pairs), "n_repeats": args.repeats} - uta_remote_tps = 1 # well-established from Hart 2020 / general experience + # --- cdot local JSON ------------------------------------------------------- + t0 = time.perf_counter() + json_provider = JSONDataProvider([args.refseq_grch38], seqfetcher=shared_seqfetcher) + load_time_s = time.perf_counter() - t0 + print(f"cdot local JSON load: {load_time_s:.1f}s") + r = run_config("cdot local JSON", lambda: json_provider, hp, pairs, args.repeats) + facts.update(cdot_local_tps=round(r["median"]), cdot_local_tps_q1=round(r["q1"]), + cdot_local_tps_q3=round(r["q3"]), grch38_load_time_s=round(load_time_s, 1)) + + # --- cdot REST ------------------------------------------------------------- + if not args.skip_rest: + def make_rest(): + return RESTDataProvider(seqfetcher=shared_seqfetcher) + + r = run_config("cdot REST (cold cache, one request per transcript)", make_rest, + hp, pairs, args.repeats, fresh_provider_per_repeat=True) + facts.update(cdot_rest_tps=round(r["median"]), cdot_rest_tps_q1=round(r["q1"]), + cdot_rest_tps_q3=round(r["q3"])) + + tx_acs = sorted({c.split(":", 1)[0] for _g, c in pairs if ":" in c}) + + def make_rest_prefetched(): + hdp = RESTDataProvider(seqfetcher=shared_seqfetcher) + t0 = time.perf_counter() + n = hdp.prefetch(tx_acs) + print(f" prefetched {n}/{len(tx_acs)} transcripts in {time.perf_counter() - t0:.2f}s (untimed)") + return hdp + + r = run_config("cdot REST (after one batch prefetch())", make_rest_prefetched, + hp, pairs, args.repeats, fresh_provider_per_repeat=True) + facts.update(cdot_rest_prefetch_tps=round(r["median"]), + cdot_rest_prefetch_tps_q1=round(r["q1"]), + cdot_rest_prefetch_tps_q3=round(r["q3"])) + + # --- UTA local ------------------------------------------------------------- + import hgvs.dataproviders.uta as uta if args.uta_uri: - print("Benchmarking UTA remote ...") - uta_remote_tps = benchmark_uta_remote(args.uta_uri) - print(f" Throughput: {uta_remote_tps:,} transcripts/s") - - facts = { - "cdot_local_min_tps": [min(cdot_local_tps, 500)], # conservative bound - "cdot_local_max_tps": [max(cdot_local_tps, 1000)], # upper bound - "cdot_rest_tps": [0], # TODO: benchmark REST API separately - "uta_local_tps": [0], # TODO: benchmark UTA local if available - "uta_remote_tps": [uta_remote_tps], - "grch38_load_time_s": [load_time_s], - } + def make_uta_local(): + hdp = uta.connect(db_url=args.uta_uri) + hdp.seqfetcher = shared_seqfetcher # identical sequence layer + return hdp + + r = run_config("UTA local PostgreSQL", make_uta_local, hp, pairs, args.repeats) + facts.update(uta_local_tps=round(r["median"]), uta_local_tps_q1=round(r["q1"]), + uta_local_tps_q3=round(r["q3"])) + + # --- UTA public remote (subset: ~0.1 HGVS/s) ------------------------------- + if not args.skip_uta_remote: + remote_pairs = pairs[:UTA_REMOTE_N] + + def make_uta_remote(): + hdp = uta.connect() # default public uta.biocommons.org + hdp.seqfetcher = shared_seqfetcher + return hdp + + print(f"(public remote UTA measured on the first {len(remote_pairs)} pairs of the same set)") + r = run_config("UTA public remote", make_uta_remote, hp, remote_pairs, + args.repeats, warmup=False) + facts.update(uta_remote_tps=round(r["median"], 2), + uta_remote_tps_q1=round(r["q1"], 2), + uta_remote_tps_q3=round(r["q3"], 2), + uta_remote_n=len(remote_pairs)) out = Path("output/facts/benchmark.csv") out.parent.mkdir(parents=True, exist_ok=True) - pd.DataFrame(facts).to_csv(out, index=False) + pd.DataFrame([facts]).to_csv(out, index=False) print(f"Written: {out}") diff --git a/paper/scripts/compute_submitted_version_age.py b/paper/scripts/compute_submitted_version_age.py new file mode 100644 index 0000000..f81af12 --- /dev/null +++ b/paper/scripts/compute_submitted_version_age.py @@ -0,0 +1,203 @@ +#!/usr/bin/env python3 +"""Version-age distribution of the submitted-HGVS ClinVar corpus (paper R2). + +Of the SCV-submitted strings citing RefSeq transcripts (built by +build_clinvar_submitted_pairs.py), how many cite a transcript version that is no +longer the current version, and how many cite a version absent from the current +RefSeq annotation release? This is the number that makes the historical-depth +argument concrete: variant_summary's Name column is always at the current +version, but what laboratories actually submitted is not. + +Reference for "current": the released cdot RefSeq GRCh38 JSON. cdot merges +annotation releases newest-wins and records each transcript's winning source +``url``, so a transcript entry whose URL is the newest RefSeq annotation release +(auto-detected as the highest RS_YYYY_MM tag among ``annotation_releases`` URLs, +e.g. RS_2025_08) is in the current annotation; every other entry exists only +because cdot carries historical releases (or UTA-derived alignments). Because a +RefSeq annotation release contains exactly one version per accession: + + * version_current cited version == the current release's version + * version_not_current accession is in the current release at a different + (newer) version: the cited version is retired + * base_retired no version of the accession is in the current release + at all (transcript dropped or suppressed) + * unversioned the submitted string cites no version + +version_not_current + base_retired together are the strings whose cited version +is absent from the current annotation, i.e. resolvable only via historical +transcript data. Of those, the script also reports how many the cdot GRCh38 file +actually carries, and how many of the remainder appear in the merged all-builds +file (GRCh37-era depth). + +The cdot JSONs are scanned with a streaming regex (transcript key -> first +in-entry URL), not json.load, so peak memory stays low. + +Counts are over unique (AlleleID, string) pairs; an SCV-weighted not-current +percentage (using the corpus scv_count column) is also emitted. + +Usage: + python paper/scripts/compute_submitted_version_age.py \ + clinvar_submitted_pairs.GRCh38.tsv \ + --refseq-grch38 cdot-0.2.34.refseq.GRCh38.json.gz \ + --refseq-allbuilds cdot-0.2.34.all-builds-refseq-....json.gz \ + --out output/facts/clinvar_submitted_version_age.csv +""" +import argparse +import csv +import gzip +import re +import sys +from collections import defaultdict +from pathlib import Path + +_KEY_URL = re.compile( + r'"((?:NM|NR|XM|XR)_\d+\.\d+)": \{"biotype"|"url": "([^"]+)"' +) +_RS_TAG = re.compile(r"RS_\d{4}_\d{2}") +_CORPUS_TX = re.compile(r"^((?:NM|NR|XM|XR)_\d+)(?:\.(\d+))?:") + + +def scan_transcript_urls(json_gz, keys_only=False): + """Stream a cdot JSON.gz, returning {accession.version: first url in entry} + (url None with keys_only). Uses a chunked regex scan with overlap so the + decompressed file is never held in memory.""" + tx_url = {} + pending = None + tail = "" + with gzip.open(json_gz, "rt") as fh: + while True: + chunk = fh.read(16 << 20) + if not chunk: + break + text = tail + chunk + # Only trust matches that start before the retained tail region; + # the last 4KB is rescanned with the next chunk. + safe_end = len(text) - 4096 if len(chunk) == (16 << 20) else len(text) + for m in _KEY_URL.finditer(text): + if m.start() >= safe_end: + break + key, url = m.group(1), m.group(2) + if key: + pending = key + tx_url[key] = None + elif pending and tx_url.get(pending) is None and not keys_only: + tx_url[pending] = url + pending = None + tail = text[safe_end:] + return tx_url + + +def current_release_tag(urls): + """Highest RS_YYYY_MM tag among annotation-release URLs (the current release).""" + tags = {t for u in urls if u and "annotation_releases" in u + for t in _RS_TAG.findall(u)} + if not tags: + sys.exit("no RS_YYYY_MM annotation_releases URLs found; is this a RefSeq cdot JSON?") + return max(tags) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("pairs", help="submitted-pairs corpus TSV (build_clinvar_submitted_pairs.py)") + ap.add_argument("--refseq-grch38", required=True, help="cdot RefSeq GRCh38 JSON.gz (reference for current/historical)") + ap.add_argument("--refseq-allbuilds", help="cdot all-builds RefSeq JSON.gz (extra depth check)") + ap.add_argument("--out", default="output/facts/clinvar_submitted_version_age.csv") + args = ap.parse_args() + + print(f"scanning {args.refseq_grch38} ...", file=sys.stderr) + tx_url = scan_transcript_urls(args.refseq_grch38) + tag = current_release_tag(tx_url.values()) + print(f" {len(tx_url):,} transcript versions; current release tag: {tag}", file=sys.stderr) + + versions_of = defaultdict(dict) # base -> {version:int -> url} + for acc, url in tx_url.items(): + base, ver = acc.rsplit(".", 1) + versions_of[base][int(ver)] = url + + def is_current(base, ver): + url = versions_of[base].get(ver) + return bool(url and tag in url and "annotation_releases" in url) + + current_version_of = { + base: next((v for v, u in vers.items() if u and tag in u and "annotation_releases" in u), None) + for base, vers in versions_of.items() + } + + n_pairs = n_refseq = n_unversioned = 0 + n_current = n_not_current = n_base_retired = 0 + n_not_current_in_cdot = n_absent_cdot = 0 + scv_refseq_versioned = scv_not_current = 0 + absent_accessions = set() + with open(args.pairs) as fh: + reader = csv.DictReader(fh, delimiter="\t") + for row in reader: + n_pairs += 1 + m = _CORPUS_TX.match(row["c_hgvs"]) + if not m: + continue + n_refseq += 1 + base, ver = m.group(1), m.group(2) + scv = int(row.get("scv_count") or 1) + if ver is None: + n_unversioned += 1 + continue + ver = int(ver) + scv_refseq_versioned += scv + in_cdot = ver in versions_of.get(base, {}) + if is_current(base, ver): + n_current += 1 + else: + n_not_current += 1 + scv_not_current += scv + if current_version_of.get(base) is None: + n_base_retired += 1 + if in_cdot: + n_not_current_in_cdot += 1 + if not in_cdot: + n_absent_cdot += 1 + absent_accessions.add(f"{base}.{ver}") + + n_absent_in_allbuilds = 0 + if args.refseq_allbuilds and absent_accessions: + print(f"scanning {args.refseq_allbuilds} (keys only) ...", file=sys.stderr) + allbuilds = scan_transcript_urls(args.refseq_allbuilds, keys_only=True) + n_absent_in_allbuilds = sum(1 for a in absent_accessions if a in allbuilds) + + n_versioned = n_refseq - n_unversioned + + def pct(n, d): + return round(100 * n / d, 2) if d else 0.0 + + row = { + "n_pairs": n_pairs, + "n_refseq_pairs": n_refseq, + "n_unversioned": n_unversioned, + "n_versioned": n_versioned, + "current_release_tag": tag, + "n_version_current": n_current, + "version_current_pct": pct(n_current, n_versioned), + "n_version_not_current": n_not_current, + "version_not_current_pct": pct(n_not_current, n_versioned), + "n_base_retired": n_base_retired, + "base_retired_pct": pct(n_base_retired, n_versioned), + "n_not_current_in_cdot": n_not_current_in_cdot, + "not_current_in_cdot_pct": pct(n_not_current_in_cdot, n_not_current), + "n_absent_cdot_grch38": n_absent_cdot, + "absent_cdot_grch38_pct": pct(n_absent_cdot, n_versioned), + "n_absent_but_in_allbuilds": n_absent_in_allbuilds, + "scv_weighted_not_current_pct": pct(scv_not_current, scv_refseq_versioned), + } + 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:<32} {v}") + print(f"Written: {out}") + + +if __name__ == "__main__": + main() diff --git a/paper/scripts/compute_version_stability.py b/paper/scripts/compute_version_stability.py index 13b85fc..bc83414 100644 --- a/paper/scripts/compute_version_stability.py +++ b/paper/scripts/compute_version_stability.py @@ -16,6 +16,12 @@ reuses the packaged cdot.hgvs.version_safety helpers, so the facts and the shipped safety check are computed by the same code. +A second facts file, ``output/facts/positional_drift.csv``, bins preserved and +total coding bases by relative CDS position (deciles, 5'→3') across the same +version-bump pairs. Because drift is overwhelmingly whole-CDS (a relocation, +position independent), the positional curve is emitted both conditioned on the +partial-drift pairs (where any positional effect must live) and unconditioned. + Per CLAUDE.md, run this only against production release files on a dedicated run, never the GTF/GFF generation pipeline. Sampling (--sample) keeps it cheap and deterministic (--seed). @@ -34,6 +40,7 @@ import json import random import re +from bisect import bisect_right from collections import Counter, defaultdict from pathlib import Path @@ -103,6 +110,41 @@ def _preserved_fraction(a, b): return same / total if total else None +N_BINS = 10 + + +def _preserved_by_decile(a, b, n_bins=N_BINS): + """Per-decile (preserved, total) coding-base counts over the shared CDS. + + Same breakpoint walk as :func:`_preserved_fraction`, with the decile edges + added as extra breakpoints so no segment spans a bin. Bin 0 is the 5'-most + tenth of the shared CDS, bin ``n_bins - 1`` the 3'-most. Returns ``None`` in + exactly the cases _preserved_fraction returns ``None``; a contig/strand + relocation (whole-CDS drift) yields zero preserved bases in every bin. + """ + if a is None or b is None or a["has_gap"] or b["has_gap"]: + return None + L = min(a["cds_len"], b["cds_len"]) + if L <= 0: + return None + relocated = a["contig"] != b["contig"] or a["strand"] != b["strand"] + edges = [i * L // n_bins for i in range(1, n_bins)] + bps = sorted({0, L} | set(edges) + | {o for o, _g in a["segs"] if 0 < o < L} + | {o for o, _g in b["segs"] if 0 < o < L}) + pres = [0] * n_bins + tot = [0] * n_bins + for x0, x1 in zip(bps, bps[1:]): + length = x1 - x0 + if length <= 0: + continue + bin_ = bisect_right(edges, x0) + tot[bin_] += length + if not relocated and _genomic_at(a, x0) == _genomic_at(b, x0): + pres[bin_] += length + return pres, tot + + def drift_stats(transcripts, build, sample, seed): """Single-build drift, per-variant safety and concentration for one consortium.""" by_acc = _by_acc(transcripts) @@ -112,6 +154,10 @@ def drift_stats(transcripts, build, sample, seed): pairs = 0 preserving = full_drift = partial_drift = 0 bases_pres = bases_total = 0 + deciles = { + "all_pres": [0] * N_BINS, "all_tot": [0] * N_BINS, + "partial_pres": [0] * N_BINS, "partial_tot": [0] * N_BINS, + } acc_drifts = Counter() acc_seen = set() for acc in accs: @@ -131,6 +177,16 @@ def drift_stats(transcripts, build, sample, seed): overlap = min(a["cds_len"], b["cds_len"]) bases_total += overlap bases_pres += round(frac * overlap) + dec = _preserved_by_decile(a, b) + if dec is not None: + dec_pres, dec_tot = dec + for i in range(N_BINS): + deciles["all_pres"][i] += dec_pres[i] + deciles["all_tot"][i] += dec_tot[i] + if 0.0 < frac < 1.0: + for i in range(N_BINS): + deciles["partial_pres"][i] += dec_pres[i] + deciles["partial_tot"][i] += dec_tot[i] if frac == 1.0: preserving += 1 elif frac == 0.0: @@ -150,7 +206,7 @@ def drift_stats(transcripts, build, sample, seed): def pct(x): return round(100 * x / pairs, 1) if pairs else 0.0 - return { + stats = { "pairs": pairs, "preserving_pct": pct(preserving), "full_drift_pct": pct(full_drift), @@ -159,6 +215,36 @@ def pct(x): "accessions_drift_pct": round(100 * n_drift_acc / n_acc, 1) if n_acc else 0.0, "top5_drift_share_pct": round(top5_share, 1), } + deciles["partial_pairs"] = partial_drift + return stats, deciles + + +def positional_facts(label, deciles): + """Flatten one consortium's decile accumulators into positional-drift facts. + + ``{label}_partial_decile{i}_pct`` is the preserved fraction of coding bases in + decile i (1 = 5'-most tenth of the shared CDS, 10 = 3'-most) across the + partial-drift pairs only; ``{label}_all_decile{i}_pct`` is the unconditioned + curve over every compared pair. Halves summarise monotonicity. + """ + def curve(pres, tot, prefix): + out = {} + for i in range(N_BINS): + out[f"{prefix}_decile{i + 1}_pct"] = ( + round(100 * pres[i] / tot[i], 1) if tot[i] else 0.0) + h = N_BINS // 2 + p5, t5 = sum(pres[:h]), sum(tot[:h]) + p3, t3 = sum(pres[h:]), sum(tot[h:]) + out[f"{prefix}_5p_half_pct"] = round(100 * p5 / t5, 1) if t5 else 0.0 + out[f"{prefix}_3p_half_pct"] = round(100 * p3 / t3, 1) if t3 else 0.0 + return out + + facts = {f"{label}_partial_pairs": deciles["partial_pairs"], + f"{label}_partial_bases": sum(deciles["partial_tot"])} + facts.update(curve(deciles["partial_pres"], deciles["partial_tot"], + f"{label}_partial")) + facts.update(curve(deciles["all_pres"], deciles["all_tot"], f"{label}_all")) + return facts def crossbuild_stats(transcripts, build, other, sample, seed): @@ -220,13 +306,15 @@ def main(): args = ap.parse_args() facts = {"sample_n": args.sample, "build": args.build} + pos_facts = {"sample_n": args.sample, "build": args.build} for label, path in (("refseq", args.refseq_grch38), ("ensembl", args.ensembl_grch38)): if path and Path(path).exists(): print(f"drift: loading {path} ...", flush=True) - d = drift_stats(_load(path), args.build, args.sample, args.seed) + d, deciles = drift_stats(_load(path), args.build, args.sample, args.seed) for k, v in d.items(): facts[f"{label}_{k}"] = v + pos_facts.update(positional_facts(label, deciles)) for label, path in (("refseq", args.refseq_allbuilds), ("ensembl", args.ensembl_allbuilds)): if path and Path(path).exists(): @@ -242,6 +330,13 @@ def main(): for k, v in facts.items(): print(f" {k}: {v}") + if len(pos_facts) > 2: # got at least one single-build file + pos_out = Path("output/facts/positional_drift.csv") + pd.DataFrame([pos_facts]).to_csv(pos_out, index=False) + print(f"Written: {pos_out}") + for k, v in pos_facts.items(): + print(f" {k}: {v}") + if __name__ == "__main__": main() diff --git a/paper/scripts/inject_and_clean.py b/paper/scripts/inject_and_clean.py index a61822a..330a7f6 100644 --- a/paper/scripts/inject_and_clean.py +++ b/paper/scripts/inject_and_clean.py @@ -72,32 +72,37 @@ # --------------------------------------------------------------------------- # Injection weights — FROZEN CONSTANTS (Tier 2, cited not shipped) -# Source: cdot_private/output/cleaning_analysis_20260617.txt -# "Which clean ops do the rescuing?" (rescued rows, ops counted; N=32,752). -# Used only to weight the reproducible per-class recovery into a single -# headline number that reflects the real-world mix of errors. No corpus string -# is read or copied. +# Source: deterministic run of clean_hgvs() plus the provider-aware +# accession-prefix restoration over the production corpus, 2026-08-17 +# ("Which clean ops do the rescuing?"; rescued rows, WARNING ops counted; +# N=32,752). Used only to weight the reproducible per-class recovery into a +# single headline number that reflects the real-world mix of errors. No +# corpus string is read or copied. +# Ops with no string-level injector are excluded: RECONSTRUCTED_STRUCTURE (248 +# rescued rows; its perturbations are covered by the kind/colon injectors), +# DROPPED_EMPTY_VERSION (8; dropping the version is not round-trippable) and +# ADDED_ACCESSION_PREFIX (1; needs a data provider). # --------------------------------------------------------------------------- REAL_RESCUE_OP_COUNTS = { - "STRIPPED_WHITESPACE": 940, - "UPPERCASED_BASES": 553, + "STRIPPED_WHITESPACE": 664, + "FIXED_SEPARATOR_TYPO": 349, + "SWAPPED_GENE_TRANSCRIPT": 156, "STRIPPED_PROTEIN_SUFFIX": 130, "FIXED_GENE_WRAPPER": 100, - "SWAPPED_GENE_TRANSCRIPT": 99, "STRIPPED_SURROUNDING_PUNCTUATION": 84, - "FIXED_MULTIPLE_COLON": 60, + "FIXED_MULTIPLE_COLON": 63, "STRIPPED_UNBALANCED_BRACKETS": 58, - "FIXED_SEPARATOR_TYPO": 37, - "FIXED_MULTIPLE_DOT": 35, - "STRIPPED_LEADING_JUNK": 23, + "UPPERCASED_BASES": 42, + "FIXED_MULTIPLE_DOT": 37, + "STRIPPED_LEADING_JUNK": 32, "FIXED_MULTIPLE_KIND": 13, "UPPERCASED_TRANSCRIPT": 10, "DROPPED_DEL_DUP_COUNT": 8, - "ADDED_TRANSCRIPT_UNDERSCORE": 6, "FIXED_PREFIX_COLON": 3, "LOWERCASED_MUTATION_TYPE": 3, "ADDED_N_PREFIX": 2, + "ADDED_TRANSCRIPT_UNDERSCORE": 2, } # --------------------------------------------------------------------------- diff --git a/paper/scripts/make_positional_figure.py b/paper/scripts/make_positional_figure.py new file mode 100644 index 0000000..3a56231 --- /dev/null +++ b/paper/scripts/make_positional_figure.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python3 +"""Render Figure S1 (positional drift along the CDS) from positional_drift.csv. + +Draws the per-decile preserved-coding-base curves produced by +compute_version_stability.py as a two-panel SVG line chart: panel A conditioned +on the partial-drift version bumps (where any positional effect must live, since +whole-CDS relocation is position independent), panel B the unconditioned curve +over every compared bump. Shared 0-100% y scale so the contrast between the two +panels is visible at a glance. + +Dependency-free (writes SVG directly, like the hand-built figure1.svg; the repo +deliberately has no plotting library). Colors are categorical slots 1 and 2 of +the validated reference dataviz palette (blue #2a78d6 / green #008300, a +documented CVD-safe adjacent pair on a light surface); series identity is also +carried by marker shape (circle vs square), direct end labels, and the legend, +never by color alone. Static print figure: light surface only. + +Usage: + python paper/scripts/make_positional_figure.py \ + [--facts paper/empirical_results/positional_drift.csv] \ + [--out paper/figures/figure_s1_positional_drift.svg] +""" +import argparse +import csv +from pathlib import Path + +# Chart tokens (reference dataviz palette, light mode, print surface) +INK = "#0b0b0b" +INK2 = "#52514e" +MUTED = "#898781" +GRID = "#e1e0d9" +AXIS = "#c3c2b7" +SURFACE = "#ffffff" +SERIES = {"refseq": ("#2a78d6", "RefSeq", "circle"), + "ensembl": ("#008300", "Ensembl", "square")} +FONT = "font-family=\"system-ui, -apple-system, 'Segoe UI', sans-serif\"" + +N_BINS = 10 +# Panel geometry +PLOT_W, PLOT_H = 280, 200 +MARGIN_L, MARGIN_T = 64, 64 +PANEL_GAP = 96 +WIDTH = MARGIN_L + PLOT_W + PANEL_GAP + PLOT_W + 48 +HEIGHT = MARGIN_T + PLOT_H + 58 + + +def read_curves(path): + with open(path, newline="") as fh: + row = next(csv.DictReader(fh)) + curves = {} + for label in SERIES: + for cond in ("partial", "all"): + curves[(label, cond)] = [ + float(row[f"{label}_{cond}_decile{i}_pct"]) for i in range(1, N_BINS + 1)] + n_pairs = {label: int(row[f"{label}_partial_pairs"]) for label in SERIES} + return curves, n_pairs + + +def xpos(x0, i): + return x0 + (i + 0.5) * PLOT_W / N_BINS + + +def ypos(pct): + return MARGIN_T + PLOT_H * (1 - pct / 100.0) + + +def marker(shape, cx, cy, color): + ring = f'stroke="{SURFACE}" stroke-width="2"' + if shape == "square": + return (f'') + return f'' + + +def panel(x0, title, subtitle, curves, cond, end_labels): + parts = [f'{title}', + f'{subtitle}'] + # y gridlines + ticks + for pct in (0, 25, 50, 75, 100): + y = ypos(pct) + parts.append(f'') + parts.append(f'{pct}') + # x axis (baseline) + tick labels at each decile + y_base = ypos(0) + parts.append(f'') + for i in range(N_BINS): + parts.append(f'{i + 1}') + parts.append(f'" + "CDS position decile (5′ → 3′)") + # series lines, markers, optional end labels + for label, (color, name, shape) in SERIES.items(): + vals = curves[(label, cond)] + pts = [(xpos(x0, i), ypos(v)) for i, v in enumerate(vals)] + d = " ".join(f"{'M' if i == 0 else 'L'}{px:.1f},{py:.1f}" + for i, (px, py) in enumerate(pts)) + parts.append(f'') + for px, py in pts: + parts.append(marker(shape, px, py, color)) + if end_labels: + px, py = pts[-1] + parts.append(f'{name} {vals[-1]:.0f}%') + return "\n".join(parts) + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--facts", default="paper/empirical_results/positional_drift.csv") + ap.add_argument("--out", default="paper/figures/figure_s1_positional_drift.svg") + args = ap.parse_args() + + curves, n_pairs = read_curves(args.facts) + x_a = MARGIN_L + x_b = MARGIN_L + PLOT_W + PANEL_GAP + + # Legend inside panel A's empty lower-left corner (stacked, clear of the curves) + legend = [] + lx, ly = x_a + 14, ypos(18) + for label, (color, name, shape) in SERIES.items(): + legend.append(f'') + legend.append(marker(shape, lx + 9, ly, color)) + legend.append(f'{name}') + ly += 20 + + svg = f""" + +{''.join(legend)} +Coding bases preserved (%) +{panel(x_a, "A  Partial-drift bumps only", + f"n = {n_pairs['refseq']} RefSeq / {n_pairs['ensembl']} Ensembl pairs", + curves, "partial", end_labels=True)} +{panel(x_b, "B  All version bumps", + "unconditioned; whole-CDS drift dominates", curves, "all", end_labels=False)} + +""" + out = Path(args.out) + out.parent.mkdir(parents=True, exist_ok=True) + out.write_text(svg) + print(f"Written: {out}") + + +if __name__ == "__main__": + main() diff --git a/paper/scripts/resolve_clinvar_pass.py b/paper/scripts/resolve_clinvar_pass.py index e03a764..a9c582b 100644 --- a/paper/scripts/resolve_clinvar_pass.py +++ b/paper/scripts/resolve_clinvar_pass.py @@ -18,12 +18,16 @@ converted_g cdot's c_to_g output, or empty if it could not convert fix_codes ';'-joined HGVSFix codes fix_hgvs() would apply (empty on clean input); only populated with --with-fixes (off by default) + fixed_bucket bucket of the fix_hgvs()-repaired string, scored the same way + (equals bucket when fix_hgvs changes nothing); only populated + with --with-fixes bucket is the BASELINE resolution of the c.HGVS exactly as given (no cleaning, no version bump), scored against g_hgvs. That is what R5b (does a safe version bump preserve converted_g?) and R6 (re-normalise the incorrect bucket) both build on. -fix_codes is recorded alongside so a recovery analysis can attribute rescues -without a second pass. +fix_codes and fixed_bucket are recorded alongside so a recovery analysis can +attribute rescues (bucket != correct, fixed_bucket == correct) without a second +pass. Provider is pluggable (shared with benchmark_resolution.py): defaults to cdot REST so it runs with no local data; pass --json a cdot release for the offline run. @@ -79,21 +83,22 @@ def parse_tx_version(c_hgvs): return m.group(1), (int(m.group(2)) if m.group(2) else None) -def fix_codes_for(provider, build, c_hgvs, version_fallback): - """Return the ';'-joined WARNING-level HGVSFix codes fix_hgvs() applies to - c_hgvs (empty string when nothing fires, e.g. already-clean input).""" +def run_fix(provider, build, c_hgvs, version_fallback): + """Return (fixed_string, ';'-joined WARNING-level HGVSFix codes) from + fix_hgvs() (codes empty when nothing fires, e.g. already-clean input).""" try: - _fixed, fixes = fix_hgvs(c_hgvs, provider, build, version_fallback=version_fallback) + fixed, fixes = fix_hgvs(c_hgvs, provider, build, version_fallback=version_fallback) except Exception as e: # noqa: BLE001 - keep the pass going logging.debug("fix_hgvs error on %s: %s", c_hgvs, e) - return "" - return ";".join(f.code.name for f in fixes if f.severity == HGVSFixSeverity.WARNING) + return c_hgvs, "" + return fixed, ";".join(f.code.name for f in fixes if f.severity == HGVSFixSeverity.WARNING) -FIELDNAMES = ["g_hgvs", "c_hgvs", "tx", "version", "bucket", "converted_g", "fix_codes"] +FIELDNAMES = ["g_hgvs", "c_hgvs", "tx", "version", "bucket", "converted_g", + "fix_codes", "fixed_bucket"] -def _row(g_hgvs, c_hgvs, bucket, converted, provider, build, with_fixes, version_fallback): +def _row(g_hgvs, c_hgvs, bucket, converted, fix_codes, fixed_bucket): tx, version = parse_tx_version(c_hgvs) return { "g_hgvs": g_hgvs, @@ -102,7 +107,8 @@ def _row(g_hgvs, c_hgvs, bucket, converted, provider, build, with_fixes, version "version": version if version is not None else "", "bucket": bucket, "converted_g": converted if converted is not None else "", - "fix_codes": fix_codes_for(provider, build, c_hgvs, version_fallback) if with_fixes else "", + "fix_codes": fix_codes, + "fixed_bucket": fixed_bucket, } @@ -110,7 +116,11 @@ def gen_rows(am, hp, provider, build, pairs, with_fixes, version_fallback): """Yield one result row per (g.HGVS, c.HGVS) pair (g.HGVS-string scoring).""" for g_hgvs, c_hgvs in pairs: bucket, converted = classify(am, hp, g_hgvs, c_hgvs) - yield _row(g_hgvs, c_hgvs, bucket, converted, provider, build, with_fixes, version_fallback) + fix_codes = fixed_bucket = "" + if with_fixes: + fixed, fix_codes = run_fix(provider, build, c_hgvs, version_fallback) + fixed_bucket = bucket if fixed == c_hgvs else classify(am, hp, g_hgvs, fixed)[0] + yield _row(g_hgvs, c_hgvs, bucket, converted, fix_codes, fixed_bucket) def gen_rows_vcf(am, hp, bf, provider, build, vcf_pairs, with_fixes, version_fallback): @@ -119,7 +129,11 @@ def gen_rows_vcf(am, hp, bf, provider, build, vcf_pairs, with_fixes, version_fal g_hgvs holds the reference CLNHGVS, converted_g holds cdot's VCF call (chrom-pos-ref-alt).""" for gt, g_hgvs, c_hgvs in vcf_pairs: bucket, converted = classify_vcf(am, hp, bf, gt, c_hgvs) - yield _row(g_hgvs, c_hgvs, bucket, converted, provider, build, with_fixes, version_fallback) + fix_codes = fixed_bucket = "" + if with_fixes: + fixed, fix_codes = run_fix(provider, build, c_hgvs, version_fallback) + fixed_bucket = bucket if fixed == c_hgvs else classify_vcf(am, hp, bf, gt, fixed)[0] + yield _row(g_hgvs, c_hgvs, bucket, converted, fix_codes, fixed_bucket) def stream_to_csv(rows, out_path): diff --git a/paper/supplementary.md b/paper/supplementary.md index dedeca0..768f0fe 100644 --- a/paper/supplementary.md +++ b/paper/supplementary.md @@ -80,6 +80,46 @@ VCF normaliser; {{ clinvar_vcf_residual.ambiguity_code_allele }} IUPAC ambiguity alleles (position matches, only the degenerate base differs); and {{ clinvar_vcf_residual.identity_or_symbolic }} identity/symbolic alleles. +Binning these coordinate differences by relative CDS position (the decile scheme of +Figure S1, cited position over CDS length; a re-run of the pass, which reproduced the +committed totals within six variants) shows no 3'-end concentration: correct projections +are nearly uniform across deciles +({{ clinvar_residual_positions.correct_decile10_pct | dp(1) }}% in the 3'-most decile) +and the {{ clinvar_residual_positions.incorrect_binned | commas }} binned differences +track them with a mild mid-CDS excess and a depleted 3'-most decile +({{ clinvar_residual_positions.incorrect_decile10_pct | dp(1) }}%), consistent with +their representation and multi-mapping origins rather than positional alignment drift +(`bin_residual_positions.py`). + +### Figure S1: Positional drift along the CDS across version bumps + +![](paper/figures/figure_s1_positional_drift.svg) + +**Figure S1. Coordinate preservation by relative CDS position across consecutive +transcript version bumps** (GRCh38, seeded +{{ version_stability.sample_n | commas }}-accession sample per consortium; deciles of the +shared CDS, 1 = 5'-most tenth). **(A)** Conditioned on the partial-drift bumps +({{ positional_drift.refseq_partial_pairs }} RefSeq and +{{ positional_drift.ensembl_partial_pairs }} Ensembl pairs, the +{{ version_stability.refseq_partial_drift_pct | dp(1) }}% and +{{ version_stability.ensembl_partial_drift_pct | dp(1) }}% of bumps where some coding +bases move and others do not): preservation declines monotonically toward the 3' end, +from {{ positional_drift.refseq_partial_decile1_pct | dp(1) }}% to +{{ positional_drift.refseq_partial_decile10_pct | dp(1) }}% of coding bases for RefSeq +and {{ positional_drift.ensembl_partial_decile1_pct | dp(1) }}% to +{{ positional_drift.ensembl_partial_decile10_pct | dp(1) }}% for Ensembl, because a +partial drift keeps a 5' prefix intact up to its first alignment change and shifts the +bases downstream of it. **(B)** Unconditioned over every compared bump the curve is +essentially flat (RefSeq {{ positional_drift.refseq_all_decile1_pct | dp(1) }}% to +{{ positional_drift.refseq_all_decile10_pct | dp(1) }}%; Ensembl +{{ positional_drift.ensembl_all_decile1_pct | dp(1) }}% to +{{ positional_drift.ensembl_all_decile10_pct | dp(1) }}%): most drift is a whole-CDS +relocation, position independent and reliably flagged by the intrinsic-structure check +(Results R5). The positional effect is confined to the rare partial-drift tail, which is +what makes a 5' coding variant safer to substitute than a 3' one when a version must be +swapped. Produced by `compute_version_stability.py` (facts) and +`make_positional_figure.py` (rendering). + ### Table S7: `clean_hgvs()` operation catalogue The cleaning pipeline (Methods) applies these operation groups in a fixed canonical order. @@ -95,9 +135,12 @@ Each operation inspects the string, makes at most one class of change, and recor (`NM_000059..4` → `NM_000059.4`); repairing a misplaced colon in the accession prefix; normalising a gene symbol wedged between extra colons or stray parentheses so that `NM_000059.4:(BRCA2):c.…` and `BRCA1(NM_000059.4)c.…` become the canonical - `transcript(GENE):c.…`; collapsing a doubled kind token (`c.c.` → `c.`); and fixing a - comma or colon used in place of the kind dot, or a period used in place of a - substitution `>`. + `transcript(GENE):c.…`; collapsing a doubled kind token (`c.c.` → `c.`); fixing a + comma or colon used in place of the kind dot, a period used in place of a + substitution `>`, a semicolon, underscore, or space used in place of the + reference:allele colon (`NM_000059.4;c.68del` and `BRCA2 c.68del` → + `NM_000059.4:c.68del` / `BRCA2:c.68del`); and dropping the dangling dot of an + empty transcript version (`NM_000059.:c.68del` → `NM_000059:c.68del`). - **Casing and prefixes**: uppercasing nucleotides in substitutions and del/ins/dup edits (`c.123delg` → `c.123delG`), lowercasing an uppercased mutation type while protecting gene symbols that contain those letters (so `NM_000059.4(INSR):c.…` is @@ -110,3 +153,10 @@ Each operation inspects the string, makes at most one class of change, and recor and a final step that detects and repairs the common clinical mistake of swapping the gene symbol and transcript accession (`BRCA2(NM_000059.4):c.…` → `NM_000059.4(BRCA2):c.…`). + +One further repair sits outside `clean_hgvs()` because it needs transcript data: +`resolve_missing_accession_prefix()` (applied by `fix_hgvs()` when a data provider is +supplied) restores a fully dropped RefSeq prefix (`000059.4:c.68del` → +`NM_000059.4:c.68del`) by generating the candidates the kind letter allows (`c.` → +`NM_`/`XM_`, `n.` → `NR_`/`XR_`) and applying the fix only when exactly one candidate +accession exists in the loaded data (Methods). diff --git a/tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv b/tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv new file mode 100644 index 0000000..0fd1304 --- /dev/null +++ b/tests/test_data/clinvar_hgvs/clinvar_submitted_500.tsv @@ -0,0 +1,501 @@ +chrom pos ref alt g_hgvs c_hgvs scv_count +12 49052426 T C NC_000012.12:g.49052426T>C NM_003482.3:c.1259-2A>G 1 +6 18121581 A G NC_000006.12:g.18121581A>G NM_198586.2:c.1026T>C 1 +17 1727203 G T NC_000017.11:g.1727203G>T NM_001163809.1:c.2244G>T 1 +1 169617125 G A NC_000001.11:g.169617125G>A NM_003005.3:c.384C>T 1 +11 47342666 C T NC_000011.10:g.47342666C>T NM_000256.3:c.1536G>A 1 +16 88428910 C A NC_000016.10:g.88428910C>A NM_001127464.1:c.1440C>A 1 +6 24278050 A G NC_000006.12:g.24278050A>G NM_016356.5:c.921T>C 1 +1 65608905 C A NC_000001.11:g.65608905C>A NM_002303.5:c.1752+4C>A 1 +17 48728327 G A NC_000017.11:g.48728327G>A NM_006361.5:c.267C>T 1 +2 151514449 A G NC_000002.12:g.151514449A>G NM_001271208.1:c.23122-21T>C 1 +5 37224383 T C NC_000005.10:g.37224383T>C NM_023073.3:c.2501-50A>G 1 +17 9857438 G A NC_000017.11:g.9857438G>A NM_004246.2:c.627G>A 1 +4 95156044 A G NC_000004.12:g.95156044A>G NM_001203.2:c.*1371A>G 1 +1 25563755 G T NC_000001.11:g.25563755G>T NM_015627.2:c.711G>T 1 +2 158658252 T C NC_000002.12:g.158658252T>C NM_003628.3:c.2031T>C 1 +18 58537781 A G NC_000018.10:g.58537781A>G NM_052947.3:c.2406T>C 1 +1 231421655 C T NC_000001.11:g.231421655C>T NM_022051.2:c.234G>A 1 +19 15192300 T A NC_000019.10:g.15192300T>A NM_000435.2:c.341-2A>T 1 +21 46154234 G A NC_000021.9:g.46154234G>A NM_001320412.1:c.153C>T 1 +11 7043322 C A NC_000011.10:g.7043322C>A NM_176822.3:c.1296C>A 1 +1 32205186 G C NC_000001.11:g.32205186G>C NM_024296.4:c.549-8G>C 1 +19 47362702 T C NC_000019.10:g.47362702T>C NM_014681.5:c.1593+9T>C 1 +2 26462177 C T NC_000002.12:g.26462177C>T NM_194248.3:c.5197G>A 1 +1 154590374 A G NC_000001.11:g.154590374A>G NM_001111.4:c.2306T>C 1 +15 28272315 C T NC_000015.10:g.28272315C>T NM_004667.5:c.983G>A 1 +16 9764088 C G NC_000016.10:g.9764088C>G NM_000833.3:c.3456G>C 1 +14 23389617 G A NC_000014.9:g.23389617G>A NM_002471.3:c.3835C>T 1 +8 22532282 T C NC_000008.11:g.22532282T>C NM_001243975.1:c.1199T>C 1 +3 32158937 C T NC_000003.12:g.32158937C>T NM_015141.3:c.680C>T 1 +17 48728396 GC G NC_000017.11:g.48728397del NM_006361.5:c.197delG 1 +1 10948070 C T NC_000001.11:g.10948070C>T NM_001170754.1:c.2065G>A 1 +17 8014650 T C NC_000017.11:g.8014650T>C NM_000180.3:c.2462T>C 1 +20 36209711 G T NC_000020.11:g.36209711G>T NM_012156.2:c.1892G>T 1 +1 19227410 G T NC_000001.11:g.19227410G>T NM_015047.1:c.2105C>A 1 +2 98396851 G GGACGC NC_000002.12:g.98396852_98396853insACGCG NM_001298.3:c.1682_1683insACGCG 1 +22 39662398 G GCGGGGAT NC_000022.11:g.39662401_39662407dup NM_021096.3:c.3338_3344dup 1 +16 2317288 C T NC_000016.10:g.2317288C>T NM_001089.2:c.1106G>A 1 +9 130703086 C T NC_000009.12:g.130703086C>T NM_014285.5:c.706C>T 1 +16 89742839 G C NC_000016.10:g.89742839G>C NM_000135.2:c.3726C>G 1 +3 173604896 C T NC_000003.12:g.173604896C>T NM_014932.5:c.298C>T 1 +11 103177704 T C NC_000011.10:g.103177704T>C NM_001377.2:c.6023T>C 1 +X 41217325 G A NC_000023.11:g.41217325G>A NM_001039590.2:c.6191G>A 1 +10 62092842 CAA C NC_000010.11:g.62092844_62092845del NM_032199.3:c.3381_3382del 1 +3 121541414 G A NC_000003.12:g.121541414G>A NM_199420.3:c.409C>T 1 +3 142508104 T G NC_000003.12:g.142508104T>G NM_001184.3:c.4858A>C 1 +11 108193979 C G NC_000011.10:g.108193979C>G NM_002519.2:c.195G>C 1 +14 75047023 G T NC_000014.9:g.75047023G>T NM_001040108.1:c.2633C>A 1 +3 47408963 G T NC_000003.12:g.47408963G>T NM_015466.2:c.1518G>T 1 +17 57105864 G T NC_000017.11:g.57105864G>T NM_001242903.1:c.400G>T 1 +19 56088435 A C NC_000019.10:g.56088435A>C NM_001002836.2:c.737T>G 1 +1 155056992 G T NC_000001.11:g.155056992G>T NM_207197.1:c.1039G>T 1 +11 74457389 T C NC_000011.10:g.74457389T>C NM_005472.4:c.175A>G 1 +16 11269054 C G NC_000016.10:g.11269054C>G NM_005425.4:c.209G>C 1 +14 74798640 C T NC_000014.9:g.74798640C>T NM_019589.2:c.3343C>T 1 +17 27597347 C T NC_000017.11:g.27597347C>T NM_014238.1:c.968C>T 1 +13 102744944 C T NC_000013.11:g.102744944C>T NM_001146197.1:c.5753G>A 1 +16 2762302 C T NC_000016.10:g.2762302C>T NM_016333.3:c.1774C>T 1 +11 66315062 C T NC_000011.10:g.66315062C>T NM_020404.2:c.1966G>A 1 +19 2291716 G A NC_000019.10:g.2291716G>A NM_001101391.1:c.61C>T 1 +1 206900318 T G NC_000001.11:g.206900318T>G NM_001185156.1:c.267T>G 1 +10 86942851 C T NC_000010.11:g.86942851C>T NM_024756.2:c.1933G>A 1 +16 20423917 C T NC_000016.10:g.20423917C>T NM_017888.2:c.769C>T 1 +1 83883117 A C NC_000001.11:g.83883117A>C NM_024686.4:c.2389T>G 1 +7 157010265 G A NC_000007.14:g.157010265G>A NM_005515.3:c.86C>T 1 +11 62522244 A G NC_000011.10:g.62522244A>G NM_001620.1:c.12173T>C 1 +3 45968445 C T NC_000003.12:g.45968445C>T NM_024513.2:c.889G>A 1 +14 94497827 T C NC_000014.9:g.94497827T>C NM_173850.2:c.571A>G 1 +8 144413853 C A NC_000008.11:g.144413853C>A NM_130849.2:c.1316G>T 1 +11 58955624 T A NC_000011.10:g.58955624T>A NM_080661.3:c.599T>A 1 +X 108734121 C T NC_000023.11:g.108734121C>T NM_003604.2:c.2224G>A 1 +9 106928728 G A NC_000009.12:g.106928728G>A NM_021224.4:c.4816G>A 1 +3 47121790 C A NC_000003.12:g.47121790C>A NM_014159.6:c.2846G>T 1 +7 100124833 G A NC_000007.14:g.100124833G>A NM_152755.1:c.692G>A 1 +7 102934254 A G NC_000007.14:g.102934254A>G NM_001031692.2:c.341A>G 1 +11 124747591 G T NC_000011.10:g.124747591G>T NM_014312.3:c.928C>A 1 +19 49422544 C A NC_000019.10:g.49422544C>A NM_178449.3:c.227G>T 1 +19 14827590 A G NC_000019.10:g.14827590A>G NM_017506.1:c.652T>C 1 +1 85200857 C T NC_000001.11:g.85200857C>T NM_032184.1:c.140G>A 1 +15 81332794 C G NC_000015.10:g.81332794C>G NM_001080532.1:c.2928G>C 1 +2 33398468 G T NC_000002.12:g.33398468G>T NM_206943.2:c.5089G>T 1 +13 102739802 G C NC_000013.11:g.102739802G>C NM_001146197.1:c.10895C>G 1 +8 120502525 T C NC_000008.11:g.120502525T>C NM_022045.3:c.1643T>C 1 +11 66437358 C T NC_000011.10:g.66437358C>T NM_016050.3:c.305G>A 1 +1 224153329 C G NC_000001.11:g.224153329C>G NM_015176.2:c.704C>G 1 +19 52115461 C T NC_000019.10:g.52115461C>T NM_178523.3:c.1703G>A 1 +1 54609332 C T NC_000001.11:g.54609332C>T NM_176782.2:c.1694G>A 1 +1 74205651 G A NC_000001.11:g.74205651G>A NM_003838.3:c.1604G>A 1 +14 67650873 C T NC_000014.9:g.67650873C>T NM_001172.3:c.1018C>T 1 +6 159239546 G A NC_000006.12:g.159239546G>A NM_032532.2:c.4210G>A 1 +16 3115457 G A NC_000016.10:g.3115457G>A NM_003456.2:c.160G>A 1 +22 30372204 G C NC_000022.11:g.30372204G>C NM_001017437.2:c.1253G>C 1 +14 58329568 C G NC_000014.9:g.58329568C>G NM_002892.3:c.703C>G 1 +9 41986297 C T NC_000009.12:g.41986297C>T NM_001201380.1:c.1348G>A 1 +12 51269266 G C NC_000012.12:g.51269266G>C NM_001031628.1:c.13C>G 1 +11 93723151 C T NC_000011.10:g.93723151C>T NM_033395.1:c.6058C>T 1 +14 64487627 A C NC_000014.9:g.64487627A>C NM_006977.2:c.604T>G 1 +11 8640552 G A NC_000011.10:g.8640552G>A NM_014818.1:c.1388C>T 1 +22 46691726 G C NC_000022.11:g.46691726G>C NM_022766.5:c.1178C>G 1 +12 32701427 A C NC_000012.12:g.32701427A>C NM_012062.4:c.115A>C 1 +6 123366158 A G NC_000006.12:g.123366158A>G NM_006073.2:c.1298T>C 1 +1 155066896 T G NC_000001.11:g.155066896T>G NM_182689.1:c.280T>G 1 +1 43985351 C T NC_000001.11:g.43985351C>T NM_030587.2:c.901C>T 1 +7 90310089 A G NC_000007.14:g.90310089A>G NM_001039706.2:c.2677A>G 1 +10 46550037 G T NC_000010.11:g.46550037G>T NM_014696.3:c.700C>A 1 +19 35449760 A G NC_000019.10:g.35449760A>G NM_005306.2:c.46A>G 1 +3 52436007 C T NC_000003.12:g.52436007C>T NM_020163.1:c.1945G>A 1 +17 44210163 C T NC_000017.11:g.44210163C>T NM_014233.3:c.1587G>A 1 +17 17796240 G T NC_000017.11:g.17796240G>T NM_030665.3:c.3292G>T 1 +12 3633630 A T NC_000012.12:g.3633630A>T NM_001144958.1:c.1709T>A 1 +17 67552303 G A NC_000017.11:g.67552303G>A NM_012417.2:c.244G>A 1 +1 47439336 G A NC_000001.11:g.47439336G>A NM_004474.3:c.1201G>A 1 +12 15483987 T C NC_000012.12:g.15483987T>C NM_030667.2:c.89T>C 1 +19 38358295 C G NC_000019.10:g.38358295C>G NM_021185.4:c.1333C>G 1 +4 121859100 A C NC_000004.12:g.121859100A>C NM_176824.3:c.420T>G 2 +19 14515965 A G NC_000019.10:g.14515965A>G NM_006145.1:c.998T>C 1 +16 22266836 C T NC_000016.10:g.22266836C>T NM_013302.3:c.1724C>T 1 +10 122512013 G A NC_000010.11:g.122512013G>A NM_002775.4:c.1222G>A 1 +6 656466 G A NC_000006.12:g.656466G>A NM_148959.3:c.479C>T 1 +16 87644861 C T NC_000016.10:g.87644861C>T NM_020655.2:c.986C>T 1 +2 169509801 C T NC_000002.12:g.169509801C>T NM_006063.2:c.23C>T 1 +18 23842691 A G NC_000018.10:g.23842691A>G NM_198129.1:c.3544A>G 1 +1 225404470 C G NC_000001.11:g.225404470C>G NM_002296.3:c.1621G>C 1 +11 100979039 T C NC_000011.10:g.100979039T>C NM_152432.2:c.2446T>C 1 +2 70985034 T C NC_000002.12:g.70985034T>C NM_001115116.1:c.1327T>C 1 +7 117422279 G A NC_000007.14:g.117422279G>A NM_130768.2:c.286C>T 1 +3 49653356 T C NC_000003.12:g.49653356T>C NM_003458.3:c.3800T>C 1 +16 88877104 C T NC_000016.10:g.88877104C>T NM_005187.5:c.1834G>A 1 +16 724720 C G NC_000016.10:g.724720C>G NM_001031737.2:c.726G>C 1 +16 11178583 C T NC_000016.10:g.11178583C>T NM_015226.2:c.3055C>T 1 +15 43698765 A G NC_000015.10:g.43698765A>G NM_001015001.1:c.1136A>G 1 +20 38517796 C T NC_000020.11:g.38517796C>T NM_020336.2:c.1213C>T 1 +11 562703 G C NC_000011.10:g.562703G>C NM_003475.3:c.749G>C 1 +16 58041800 A G NC_000016.10:g.58041800A>G NM_002428.2:c.1094A>G 1 +3 52440441 C T NC_000003.12:g.52440441C>T NM_020163.2:c.1079G>A 1 +15 48136868 G A NC_000015.10:g.48136868G>A NM_205850.2:c.776G>A 1 +15 66781531 G T NC_000015.10:g.66781531G>T NM_005585.4:c.1487G>T 1 +19 37892995 C G NC_000019.10:g.37892995C>G NM_031951.3:c.2591G>C 1 +4 146640671 G T NC_000004.12:g.146640671G>T NM_004575.2:c.1093G>T 1 +18 58579590 C T NC_000018.10:g.58579590C>T NM_052947.3:c.1186G>A 1 +7 105565568 C A NC_000007.14:g.105565568C>A NM_021930.4:c.2106C>A 1 +6 99435841 T A NC_000006.12:g.99435841T>A NM_001080481.3:c.2320A>T 1 +X 153648119 A G NC_000023.11:g.153648119A>G NM_001395.2:c.166A>G 1 +4 3444105 C T NC_000004.12:g.3444105C>T NM_001528.2:c.542C>T 1 +9 116187805 G A NC_000009.12:g.116187805G>A NM_002581.3:c.1067G>A 1 +19 54095887 C G NC_000019.10:g.54095887C>G NM_206818.1:c.652G>C 1 +5 75677939 T A NC_000005.10:g.75677939T>A NM_001099271.1:c.1419A>T 1 +14 32091784 C T NC_000014.9:g.32091784C>T NM_001030055.1:c.1115C>T 1 +22 32406709 T C NC_000022.11:g.32406709T>C NM_014306.4:c.293A>G 1 +10 7205844 T C NC_000010.11:g.7205844T>C NM_001029880.2:c.1415A>G 1 +14 99175485 C A NC_000014.9:g.99175485C>A NM_022898.2:c.1138G>T 1 +2 25161144 CT C NC_000002.12:g.25161145del NM_000939.3:c.740delA 1 +16 2766181 C T NC_000016.10:g.2766181C>T NM_016333.4:c.5653C>T 1 +1 61088410 G T NC_000001.11:g.61088410G>T NM_005595.4:c.289G>T 1 +17 7512269 C G NC_000017.11:g.7512269C>G NM_000937.4:c.4417C>G 1 +1 158845720 T G NC_000001.11:g.158845720T>G NM_002432.1:c.704T>G 1 +19 11441295 G A NC_000019.10:g.11441295G>A NM_002743.2:c.406G>A 1 +8 68069843 A T NC_000008.11:g.68069843A>T NM_024870.2:c.1452A>T 1 +11 32954042 G T NC_000011.10:g.32954042G>T NM_001076786.1:c.3976G>T 1 +5 87268504 C T NC_000005.10:g.87268504C>T NM_002890.3:c.53C>T 1 +2 166405964 G C NC_000002.12:g.166405964G>C NM_002976.3:c.4665C>G 1 +X 135022298 G T NC_000023.11:g.135022298G>T NM_001078173.1:c.162C>A 1 +13 23320750 G C NC_000013.11:g.23320750G>C NM_000231.2:c.692G>C 1 +20 19685238 G A NC_000020.11:g.19685238G>A NM_020689.3:c.1201G>A 1 +17 35264578 C A NC_000017.11:g.35264578C>A NM_144975.3:c.1534C>A 1 +4 5446686 C G NC_000004.12:g.5446686C>G NM_018401.1:c.576C>G 1 +6 89333148 A G NC_000006.12:g.89333148A>G NM_016021.2:c.616T>C 1 +X 56563953 C T NC_000023.11:g.56563953C>T NM_013444.3:c.80C>T 1 +7 101164416 G C NC_000007.14:g.101164416G>C NM_003378.3:c.428C>G 1 +3 39184027 C T NC_000003.12:g.39184027C>T NM_194293.2:c.5419G>A 1 +3 44447582 C T NC_000003.12:g.44447582C>T NM_181489.5:c.2089G>A 1 +6 143884306 G A NC_000006.12:g.143884306G>A NM_001013623.2:c.31G>A 1 +1 159193882 G A NC_000001.11:g.159193882G>A NM_021189.3:c.635G>A 1 +11 85916201 T A NC_000011.10:g.85916201T>A NM_173556.3:c.1141T>A 1 +20 3671681 T A NC_000020.11:g.3671681T>A NM_025220.2:c.1805A>T 1 +16 70482153 G C NC_000016.10:g.70482153G>C NM_015386.2:c.1943C>G 1 +19 45364277 C A NC_000019.10:g.45364277C>A NM_000400.3:c.773G>T 1 +1 109750749 G A NC_000001.11:g.109750749G>A NM_139053.2:c.1684C>T 1 +2 108471109 G A NC_000002.12:g.108471109G>A NM_181453.3:c.1780G>A 1 +9 65283486 A G NC_000009.12:g.65283486A>G NM_001126334.1:c.892T>C 1 +7 24717287 T C NC_000007.14:g.24717287T>C NM_004403.2:c.664A>G 1 +8 36918835 G A NC_000008.11:g.36918835G>A NM_001031836.2:c.2534G>A 1 +10 24473864 C T NC_000010.11:g.24473864C>T NM_019590.3:c.1483C>T 1 +15 63643547 C G NC_000015.10:g.63643547C>G NM_003922.3:c.11188G>C 1 +10 89462407 C T NC_000010.11:g.89462407C>T NM_213606.4:c.172G>A 1 +16 31117821 C A NC_000016.10:g.31117821C>A NM_032188.3:c.140C>A 1 +19 5131933 G A NC_000019.10:g.5131933G>A NM_015015.3:c.1832G>A 1 +7 584538 C G NC_000007.14:g.584538C>G NM_001164761.1:c.739G>C 1 +3 12618646 T A NC_000003.12:g.12618646T>A NM_002880.3:c.76A>T 1 +9 89386480 G A NC_000009.12:g.89386480G>A NM_006378.3:c.1333C>T 1 +8 38820436 G C NC_000008.11:g.38820436G>C NM_006283.2:c.1192G>C 1 +2 137242537 G T NC_000002.12:g.137242537G>T NM_001080427.1:c.2138G>T 1 +9 35370345 C G NC_000009.12:g.35370345C>G NM_006377.3:c.1242C>G 1 +19 30009123 C T NC_000019.10:g.30009123C>T NM_003796.2:c.805C>T 1 +10 97750443 C A NC_000010.11:g.97750443C>A NM_001002261.3:c.777C>A 1 +12 1645956 C T NC_000012.12:g.1645956C>T NM_032642.2:c.784C>T 1 +15 62068476 G A NC_000015.10:g.62068476G>A NM_207322.2:c.863G>A 1 +2 206543278 G A NC_000002.12:g.206543278G>A NM_003812.2:c.682G>A 1 +5 129735107 T G NC_000005.10:g.129735107T>G NM_133638.3:c.3470T>G 1 +9 120536454 G A NC_000009.12:g.120536454G>A NM_018249.4:c.580C>T 1 +3 129257816 A G NC_000003.12:g.129257816A>G NM_016128.3:c.827A>G 1 +1 160239921 C A NC_000001.11:g.160239921C>A NM_015726.3:c.499G>T 1 +12 13214648 G C NC_000012.12:g.13214648G>C NM_001423.2:c.431G>C 1 +7 127615062 G C NC_000007.14:g.127615062G>C NM_006193.2:c.154C>G 1 +10 94246152 T C NC_000010.11:g.94246152T>C NM_016341.3:c.2627T>C 1 +19 18163296 G A NC_000019.10:g.18163296G>A NM_005027.2:c.1324G>A 1 +18 2694582 T C NC_000018.10:g.2694582T>C NM_015295.2:c.929T>C 1 +2 71148038 A G NC_000002.12:g.71148038A>G NM_005791.2:c.1597A>G 1 +17 29096879 G A NC_000017.11:g.29096879G>A NM_078471.3:c.4267C>T 1 +7 158652328 G A NC_000007.14:g.158652328G>A NM_017760.5:c.2899C>T 1 +9 97668905 A G NC_000009.12:g.97668905A>G NM_002486.4:c.2076A>G 1 +5 36975895 A G NC_000005.10:g.36975895A>G NM_133433.3:c.988A>G 1 +1 24458929 T C NC_000001.11:g.24458929T>C NM_020448.4:c.815T>C 1 +1 5927718 G A NC_000001.11:g.5927718G>A NM_015102.3:c.1372C>T 1 +7 150857504 G T NC_000007.14:g.150857504G>T NM_001091.2:c.1034G>T 1 +16 4693780 G A NC_000016.10:g.4693780G>A NM_032349.3:c.54G>A 1 +20 1443801 G T NC_000020.11:g.1443801G>T NM_001206736.1:c.1067C>A 1 +15 44675670 G C NC_000015.10:g.44675670G>C NM_001145112.1:c.38C>G 1 +10 68284224 T C NC_000010.11:g.68284224T>C NM_022129.3:c.820A>G 1 +5 141489477 G C NC_000005.10:g.141489477G>C NM_018929.2:c.237G>C 1 +2 158642499 G A NC_000002.12:g.158642499G>A NM_003628.3:c.1709G>A 1 +9 6012720 A G NC_000009.12:g.6012720A>G NM_012416.3:c.2888T>C 1 +1 173909616 C G NC_000001.11:g.173909616C>G NM_000488.3:c.1089G>C 1 +1 234427692 C G NC_000001.11:g.234427692C>G NM_005646.3:c.3135G>C 1 +11 1295721 T G NC_000011.10:g.1295721T>G NM_019009.3:c.107A>C 1 +15 82810121 C T NC_000015.10:g.82810121C>T NM_001080435.1:c.395C>T 1 +X 70238310 C G NC_000023.11:g.70238310C>G NM_001013579.2:c.559C>G 1 +5 77453892 A G NC_000005.10:g.77453892A>G NM_018268.2:c.448T>C 1 +6 31644192 G C NC_000006.12:g.31644192G>C NM_004639.3:c.1468C>G 1 +16 48211010 G A NC_000016.10:g.48211010G>A NM_032583.3:c.1546C>T 1 +2 112309978 A G NC_000002.12:g.112309978A>G NM_198581.2:c.430A>G 1 +2 112324571 A G NC_000002.12:g.112324571A>G NM_198581.2:c.1760A>G 1 +8 94170967 C T NC_000008.11:g.94170967C>T NM_001144663.1:c.802G>A 1 +5 79044430 C G NC_000005.10:g.79044430C>G NM_013391.2:c.868G>C 1 +11 6566370 T C NC_000011.10:g.6566370T>C NM_144666.2:c.11183T>C 1 +5 55108314 A C NC_000005.10:g.55108314A>C NM_006144.3:c.547A>C 1 +4 113355986 G C NC_000004.12:g.113355986G>C NM_001148.4:c.7368G>C 1 +11 64800767 C T NC_000011.10:g.64800767C>T NM_004579.3:c.722G>A 1 +6 142770009 T A NC_000006.12:g.142770009T>A NM_006734.3:c.4730A>T 1 +5 141209982 A G NC_000005.10:g.141209982A>G NM_018932.3:c.1075A>G 1 +3 15644541 C T NC_000003.12:g.15644541C>T NM_001370658.1:c.625C>T 1 +17 31258057 G A NC_000017.11:g.31258057G>A NM_000267.3:c.4111-287G>A 1 +2 151655888 A C NC_000002.12:g.151655888A>C NM_001271208.2:c.6631T>G 1 +20 18525912 A G NC_000020.11:g.18525912A>G NM_006363.6:c.814A>G 1 +9 136197447 T G NC_000009.12:g.136197447T>G NM_014564.3:c.1087A>C 1 +2 73230168 C T NC_000002.12:g.73230168C>T NM_032319.1:c.113G>A 1 +13 42888375 T C NC_000013.11:g.42888375T>C NM_001002264.1:c.1108A>G 1 +4 56469663 C G NC_000004.12:g.56469663C>G NM_006947.3:c.120C>G 1 +17 49798345 G A NC_000017.11:g.49798345G>A NM_007067.4:c.367G>A 1 +12 47077942 C T NC_000012.12:g.47077942C>T NM_181847.4:c.1061G>A 1 +21 44638070 C G NC_000021.9:g.44638070C>G NM_181688.1:c.653C>G 1 +20 10038483 T G NC_000020.11:g.10038483T>G NM_022096.4:c.182T>G 1 +17 58267479 C A NC_000017.11:g.58267479C>A NM_006151.2:c.1824C>A 1 +17 45397952 C G NC_000017.11:g.45397952C>G NM_199282.2:c.816G>C 1 +20 35842685 C T NC_000020.11:g.35842685C>T NM_016436.4:c.196C>T 1 +16 269324 G T NC_000016.10:g.269324G>T NM_183337.1:c.1349C>A 1 +1 210241979 C T NC_000001.11:g.210241979C>T NM_019605.3:c.713C>T 1 +17 7042574 A T NC_000017.11:g.7042574A>T NM_153357.1:c.608T>A 1 +9 15423160 T A NC_000009.12:g.15423160T>A NM_001039697.1:c.281T>A 1 +15 74202255 G C NC_000015.10:g.74202255G>C NM_022369.3:c.13C>G 1 +17 49168861 G T NC_000017.11:g.49168861G>T NM_153446.2:c.1456G>T 1 +9 129868240 C T NC_000009.12:g.129868240C>T NM_006676.6:c.926C>T 1 +19 7622370 C A NC_000019.10:g.7622370C>A NM_020196.2:c.1578G>T 1 +1 35392256 C G NC_000001.11:g.35392256C>G NM_005095.2:c.2632C>G 1 +19 58067830 A G NC_000019.10:g.58067830A>G NM_007134.1:c.1418A>G 1 +16 547527 C T NC_000016.10:g.547527C>T NM_005632.2:c.689C>T 1 +7 21906428 G A NC_000007.14:g.21906428G>A NM_018719.4:c.782C>T 1 +12 77044778 T C NC_000012.12:g.77044778T>C NM_203394.2:c.847A>G 1 +1 8865353 T C NC_000001.11:g.8865353T>C NM_001428.3:c.797A>G 1 +12 132570214 C T NC_000012.12:g.132570214C>T NM_001142641.1:c.980C>T 1 +1 108820578 G T NC_000001.11:g.108820578G>T NM_152763.3:c.2216C>A 1 +12 47789529 G A NC_000012.12:g.47789529G>A NM_015401.3:c.2141C>T 1 +9 21350385 C G NC_000009.12:g.21350385C>G NM_021002.2:c.503G>C 1 +X 118610335 A G NC_000023.11:g.118610335A>G NM_144658.3:c.3013A>G 1 +6 170289712 GC G NC_000006.12:g.170289716del NM_005618.3:c.150del 1 +15 48301386 A G NC_000015.10:g.48301386A>G NM_000338.3:c.3164+4A>G 1 +3 52348969 C T NC_000003.12:g.52348969C>T NM_015512.4:c.2188C>T 1 +20 10639899 C T NC_000020.11:g.10639899C>T NM_000214.2:c.3256G>A 1 +2 9543222 G A NC_000002.12:g.9543222G>A NM_003183.4:c.161C>T 1 +19 19110746 G A NC_000019.10:g.19110746G>A NM_178526.5:c.827G>A 1 +1 48471962 C G NC_000001.11:g.48471962C>G NM_019073.2:c.47G>C 1 +5 128114212 G A NC_000005.10:g.128114212G>A NM_001046.3:c.877G>A 1 +7 44112696 G T NC_000007.14:g.44112696G>T NM_001129.5:c.2356G>T 1 +11 65046348 G C NC_000011.10:g.65046348G>C NM_005468.2:c.1696C>G 1 +17 42021615 C T NC_000017.11:g.42021615C>T NM_001144927.1:c.38C>T 1 +12 123485720 A G NC_000012.12:g.123485720A>G NM_178314.3:c.887T>C 1 +7 77778822 A C NC_000007.14:g.77778822A>C NM_198467.2:c.2195A>C 1 +10 100499579 G A NC_000010.11:g.100499579G>A NM_015490.3:c.1430C>T 1 +19 16865566 C T NC_000019.10:g.16865566C>T NM_015260.2:c.1636C>T 1 +11 63096048 A G NC_000011.10:g.63096048A>G NM_001136506.2:c.1013T>C 1 +18 46086420 G A NC_000018.10:g.46086420G>A NM_001001937.1:c.1251C>T 1 +2 171055109 T G NC_000002.12:g.171055109T>G NM_012290.4:c.613A>C 1 +7 135187033 G A NC_000007.14:g.135187033G>A NM_014149.3:c.2018C>T 1 +17 10695795 T C NC_000017.11:g.10695795T>C NM_004589.2:c.310A>G 1 +2 121367788 C G NC_000002.12:g.121367788C>G NM_015282.2:c.3686G>C 1 +16 66609986 C T NC_000016.10:g.66609986C>T NM_144601.3:c.503C>T 1 +3 194686884 G A NC_000003.12:g.194686884G>A NM_153690.4:c.58G>A 1 +1 171324169 A G NC_000001.11:g.171324169A>G NM_002022.1:c.353A>G 1 +19 55359088 T C NC_000019.10:g.55359088T>C NM_001145402.1:c.1780A>G 1 +12 57471004 C T NC_000012.12:g.57471004C>T NM_005269.2:c.2264C>T 1 +1 18365593 T C NC_000001.11:g.18365593T>C NM_032880.4:c.911T>C 1 +6 76005346 A G NC_000006.12:g.76005346A>G NM_001563.2:c.1076T>C 1 +8 91084834 G A NC_000008.11:g.91084834G>A NM_016023.5:c.848G>A 1 +10 71446400 G A NC_000010.11:g.71446400G>A NM_022124.5:c.145+5G>A 1 +11 2588801 C A NC_000011.10:g.2588801C>A NM_000218.2:c.1340C>A 1 +12 2567572 C T NC_000012.12:g.2567572C>T NM_000719.6:c.1673C>T 2 +9 128222476 C T NC_000009.12:g.128222476C>T NM_004408.2:c.1008C>T 1 +X 111724993 G T NC_000023.11:g.111724993G>T NM_001099922.2:c.1661G>T 1 +13 23331229 T C NC_000013.11:g.23331229T>C NM_014363.5:c.12647A>G 1 +6 64590780 T G NC_000006.12:g.64590780T>G NM_001142800.2:c.5087A>C 1 +21 37625289 A G NC_000021.9:g.37625289A>G NM_002240.5:c.1142T>C 1 +16 89779955 G C NC_000016.10:g.89779955G>C NM_000135.2:c.1629C>G 1 +19 35249008 G A NC_000019.10:g.35249008G>A NM_205834.2:c.130G>A 1 +12 122490539 C T NC_000012.12:g.122490539C>T NM_017612.3:c.346G>A 1 +11 64296965 G A NC_000011.10:g.64296965G>A NM_033310.2:c.277G>A 2 +15 40932351 G A NC_000015.10:g.40932351G>A NM_019074.3:c.754G>A 1 +21 46272582 A G NC_000021.9:g.46272582A>G NM_003906.3:c.2444T>C 1 +9 77228192 C T NC_000009.12:g.77228192C>T NM_033305.2:c.1523C>T 1 +3 139348194 G A NC_000003.12:g.139348194G>A NM_020191.2:c.374G>A 1 +16 58009966 G A NC_000016.10:g.58009966G>A NM_024598.3:c.303G>A 1 +14 102027516 G A NC_000014.9:g.102027516G>A NM_001376.4:c.9020G>A 1 +X 41344350 C T NC_000023.11:g.41344350C>T NM_001356.3:c.976C>T 1 +7 25124041 G A NC_000007.14:g.25124041G>A NM_018947.6:c.79C>T 3 +10 102396726 G C NC_000010.11:g.102396726G>C NM_001077494.3:c.146G>C 2 +1 172553276 C T NC_000001.11:g.172553276C>T NM_014283.3:c.194C>T 1 +20 32435519 C T NC_000020.11:g.32435519C>T NM_015338.5:c.2807C>T 1 +11 47349904 GAGA G NC_000011.10:g.47349907_47349909del NM_000256.3:c.521_523delTCT 2 +9 127661161 ACT A NC_000009.12:g.127661162CT[1] NM_003165.3:c.388_389delCT 1 +15 51576176 A G NC_000015.10:g.51576176A>G NM_001174116.1:c.93T>C 1 +3 129437838 G C NC_000003.12:g.129437838G>C NM_001276270.2:c.217C>G 1 +22 50220869 C T NC_000022.11:g.50220869C>T NM_020461.3:c.3490G>A 2 +16 88434270 G A NC_000016.10:g.88434270G>A NM_001127464.1:c.6716G>A 2 +20 62824440 G A NC_000020.11:g.62824440G>A NM_001853.4:c.520-5G>A 1 +14 23389417 T C NC_000014.9:g.23389417T>C NM_002471.3:c.3954A>G 1 +17 58357806 C G NC_000017.11:g.58357806C>G NM_017763.6:c.1970G>C 1 +20 23049022 C A NC_000020.11:g.23049022C>A NM_000361.2:c.483G>T 1 +17 80368140 A G NC_000017.11:g.80368140A>G NM_001256071.1:c.12152A>G 1 +12 109788550 C CA NC_000012.12:g.109788551dup NM_021625.4:c.2057dupT 1 +8 67161813 CAT C NC_000008.11:g.67161814_67161815del NM_024790.6:c.2527_2528delAT 2 +6 79925034 T C NC_000006.12:g.79925034T>C NM_022726.3:c.289-2A>G 1 +16 10907778 C A NC_000016.10:g.10907778C>A NM_000246.3:c.2286C>A 2 +21 43059272 G A NC_000021.9:g.43059272G>A NM_000071.2:c.1177C>T 1 +18 57551312 T C NC_000018.10:g.57551312T>C NM_000140.4:c.1137+3A>G 1 +16 2072337 C T NC_000016.10:g.2072337C>T NM_000548.3:c.2194C>T 3 +7 117590379 A G NC_000007.14:g.117590379A>G NM_000492.3:c.1706A>G 2 +20 63350612 G A NC_000020.11:g.63350612G>A NM_000744.5:c.799C>T 1 +9 108926548 C T NC_000009.12:g.108926548C>T NM_003640.3:c.441G>A 3 +8 102212897 G T NC_000008.11:g.102212897G>T NM_015713.5:c.790-8C>A 2 +2 178651699 C T NC_000002.12:g.178651699C>T NM_133378.4:c.32128G>A 1 +3 37008894 A G NC_000003.12:g.37008894A>G NM_000249.3:c.534A>G 4 +17 61686079 G A NC_000017.11:g.61686079G>A NM_032043.2:c.2662C>T 8 +7 117535257 T C NC_000007.14:g.117535257T>C NM_000492.4:c.589T>C 1 +19 50416622 C T NC_000019.10:g.50416622C>T NM_002691.2:c.2966C>T 1 +12 132638096 C T NC_000012.12:g.132638096C>T NM_006231.4:c.5596G>A 1 +19 11111585 C T NC_000019.10:g.11111585C>T NM_000527.5:c.1132C>T 2 +4 78470102 T C NC_000004.12:g.78470102T>C NM_025074.7:c.7371+11T>C 1 +10 71778352 C G NC_000010.11:g.71778352C>G NM_022124.6:c.5187+44C>G 2 +12 21604811 A C NC_000012.12:g.21604811A>C NM_021957.3:c.-219T>G 1 +4 169393800 T C NC_000004.12:g.169393800T>C NM_012224.2:c.*710A>G 1 +6 152344164 A C NC_000006.12:g.152344164A>C NM_033071.3:c.11929T>G 3 +17 7221969 TTCTG T NC_000017.11:g.7221973_7221976del NM_000018.4:c.644_647del 1 +7 128845086 C T NC_000007.14:g.128845086C>T NM_001458.4:c.3621C>T 5 +X 32365169 C T NC_000023.11:g.32365169C>T NM_004006.2:c.4876G>A 4 +2 29193286 C T NC_000002.12:g.29193286C>T NM_004304.3:c.4801G>A 1 +8 89953246 A G NC_000008.11:g.89953246A>G NM_002485.5:c.1843T>C 2 +14 95115790 T C NC_000014.9:g.95115790T>C NM_177438.3:c.1784A>G 1 +5 240381 C A NC_000005.10:g.240381C>A NM_004168.4:c.1456C>A 1 +13 32319185 C G NC_000013.11:g.32319185C>G NM_000059.3:c.176C>G 2 +2 178601535 CTT C NC_000002.12:g.178601537_178601538del NM_001256850.1:c.50537_50538delAA 1 +2 214781054 T C NC_000002.12:g.214781054T>C NM_000465.2:c.820A>G 1 +13 32376727 C G NC_000013.11:g.32376727C>G NM_000059.3:c.8690C>G 3 +2 47445563 T TA NC_000002.12:g.47445564dup NM_000251.1:c.1293dupA 1 +2 73519760 C A NC_000002.12:g.73519760C>A NM_015120.4:c.9543-15C>A 1 +2 73519994 T TAC NC_000002.12:g.73519995_73519996dup NM_015120.4:c.9763_9764dupAC 1 +5 78885674 G A NC_000005.10:g.78885674G>A NM_000046.4:c.1052C>T 1 +7 117540155 G A NC_000007.14:g.117540155G>A NM_000492.3:c.925G>A 2 +1 198706748 A G NC_000001.11:g.198706748A>G NM_002838.3:c.694A>G 1 +9 131515444 CCT C NC_000009.12:g.131515445_131515446del NM_001077365.2:c.1195_1196del 1 +11 108247004 CT C NC_000011.10:g.108247006del NM_000051.3:c.944del 1 +12 32878474 C A NC_000012.12:g.32878474C>A NM_004572.4:c.406G>T 1 +19 50401832 T G NC_000019.10:g.50401832T>G NM_002691.2:c.371T>G 1 +8 144512557 GA G NC_000008.11:g.144512558del NM_004260.4:c.2889del 1 +2 47799318 T G NC_000002.12:g.47799318T>G NM_000179.2:c.1335T>G 1 +18 51067092 C T NC_000018.10:g.51067092C>T NM_005359.5:c.1213C>T 2 +16 2054324 T A NC_000016.10:g.2054324T>A NM_000548.5:c.365T>A 1 +18 46545308 C T NC_000018.10:g.46545308C>T NM_144612.6:c.3619+9G>A 3 +3 10073363 G A NC_000003.12:g.10073363G>A NM_033084.4:c.2715+1G>A 2 +1 40091984 G A NC_000001.11:g.40091984G>A NM_000310.3:c.362+61C>T 2 +19 48304065 C T NC_000019.10:g.48304065C>T NM_144577.3:c.630G>A 1 +9 95482182 C T NC_000009.12:g.95482182C>T NM_000264.3:c.606G>A 1 +12 39341575 C T NC_000012.12:g.39341575C>T NM_017641.3:c.1812G>A 2 +4 657483 C T NC_000004.12:g.657483C>T NM_000283.3:c.1390C>T 2 +19 50403490 C A NC_000019.10:g.50403490C>A NM_002691.2:c.1138-3C>A 1 +5 132588817 A G NC_000005.10:g.132588817A>G NM_005732.3:c.1182A>G 1 +6 56619308 C T NC_000006.12:g.56619308C>T NM_001723.5:c.4726G>A 1 +16 88842720 G C NC_000016.10:g.88842720G>C NM_000512.4:c.230C>G 2 +1 20645675 G A NC_000001.11:g.20645675G>A NM_032409.2:c.1075G>A 2 +11 31794664 G A NC_000011.10:g.31794664G>A NM_000280.4:c.648C>T 4 +12 123697081 T C NC_000012.12:g.123697081T>C NM_024809.4:c.1394-6T>C 2 +14 76500014 T C NC_000014.9:g.76500014T>C NM_004452.3:c.1448T>C 2 +16 50732832 A G NC_000016.10:g.50732832A>G NM_022162.1:c.*1013A>G 2 +16 53601225 G T NC_000016.10:g.53601225G>T NM_015272.2:c.*851C>A 3 +17 75763963 G A NC_000017.11:g.75763963G>A NM_000154.1:c.289C>T 1 +20 33412667 G A NC_000020.11:g.33412667G>A NM_003098.2:c.817C>T 2 +4 79903476 C T NC_000004.12:g.79903476C>T NM_058172.5:c.*3953G>A 1 +6 7584732 T C NC_000006.12:g.7584732T>C NM_004415.2:c.7470T>C 2 +1 237648602 T G NC_000001.11:g.237648602T>G NM_001035.3:c.7501T>G 2 +13 32316451 G A NC_000013.11:g.32316451G>A NM_000059.3:c.-10G>A 1 +6 7569191 C T NC_000006.12:g.7569191C>T NM_004415.2:c.1425C>T 2 +1 45332192 C T NC_000001.11:g.45332192C>T NM_001128425.1:c.907G>A 1 +7 116759381 G T NC_000007.14:g.116759381G>T NM_001127500.1:c.2309G>T 2 +2 214769228 C T NC_000002.12:g.214769228C>T NM_000465.4:c.1395+4G>A 1 +X 108695336 C T NC_000023.11:g.108695336C>T NM_033380.3:c.4891C>T 1 +7 87443759 G A NC_000007.14:g.87443759G>A NM_000443.3:c.1134C>T 1 +7 55205297 G A NC_000007.14:g.55205297G>A NM_005228.3:c.3313G>A 1 +2 188984776 T C NC_000002.12:g.188984776T>C NM_000090.3:c.96T>C 1 +X 154357504 G A NC_000023.11:g.154357504G>A NM_001110556.1:c.4875C>T 1 +9 12704526 G C NC_000009.12:g.12704526G>C NM_000550.3:c.1082G>C 1 +X 38301338 C T NC_000023.11:g.38301338C>T NM_000328.2:c.968G>A 2 +10 98422394 G C NC_000010.11:g.98422394G>C NM_000195.5:c.1718C>G 2 +17 75522027 C T NC_000017.11:g.75522027C>T NM_207346.3:c.946C>T 1 +14 23433665 C T NC_000014.9:g.23433665C>T NM_000257.4:c.68G>A 2 +2 233760930 A G NC_000002.12:g.233760930A>G NM_000463.3:c.643A>G 2 +18 51058145 G A NC_000018.10:g.51058145G>A NM_005359.6:c.688G>A 1 +22 49907841 G A NC_000022.11:g.49907841G>A NM_024105.3:c.872C>T 1 +7 116699246 T C NC_000007.14:g.116699246T>C NM_001127500.1:c.162T>C 1 +20 4699856 G A NC_000020.11:g.4699856G>A NM_000311.5:c.636G>A 1 +11 68416313 C T NC_000011.10:g.68416313C>T NM_002335.4:c.2828-15C>T 1 +11 78478493 AAAG A NC_000011.10:g.78478496_78478498del NM_024678.6:c.922-21_922-19delCTT 1 +11 67490106 G T NC_000011.10:g.67490106G>T NM_003977.2:c.537G>T 1 +13 32340798 C CTA NC_000013.11:g.32340800_32340801dup NM_000059.3:c.6445_6446dupAT 1 +17 65557855 T C NC_000017.11:g.65557855T>C NM_004655.3:c.766A>G 1 +22 20991769 G A NC_000022.11:g.20991769G>A NM_006767.3:c.933G>A 1 +17 43094545 T C NC_000017.11:g.43094545T>C NM_007294.4:c.986A>G 1 +17 43093249 T C NC_000017.11:g.43093249T>C NM_007294.3:c.2282A>G 1 +4 127921938 C G NC_000004.12:g.127921938C>G NM_152778.4:c.1024G>C 1 +16 2054305 T G NC_000016.10:g.2054305T>G NM_000548.3:c.346T>G 1 +5 132385477 G A NC_000005.10:g.132385477G>A NM_003060.3:c.802G>A 1 +16 89281453 G A NC_000016.10:g.89281453G>A NM_013275.6:c.5089C>T 1 +X 17728193 A G NC_000023.11:g.17728193A>G NM_198270.2:c.4024A>G 1 +1 193142038 G A NC_000001.11:g.193142038G>A NM_024529.5:c.701G>A 1 +5 74720499 A C NC_000005.10:g.74720499A>C NM_000521.4:c.1489A>C 1 +16 2086821 T G NC_000016.10:g.2086821T>G NM_000548.3:c.4939T>G 1 +15 89274245 GA G NC_000015.10:g.89274247del NM_001113378.2:c.1055del 1 +22 23803342 C T NC_000022.11:g.23803342C>T NM_003073.3:c.548C>T 1 +12 110281629 G A NC_000012.12:g.110281629G>A NM_170665.4:c.-161G>A 1 +16 67435979 C T NC_000016.10:g.67435979C>T NM_000196.4:c.501C>T 2 +6 52024709 C T NC_000006.12:g.52024709C>T NM_138694.4:c.5101G>A 1 +19 38721574 C T NC_000019.10:g.38721574C>T NM_004924.4:c.1328C>T 1 +19 50409172 A G NC_000019.10:g.50409172A>G NM_002691.2:c.1943A>G 1 +11 71444151 A G NC_000011.10:g.71444151A>G NM_001360.2:c.163T>C 1 +9 108929758 CT C NC_000009.12:g.108929760del NM_003640.4:c.303+10del 1 +6 144187411 CAGGTGCG C NC_000006.12:g.144187412_144187418del NM_003764.4:c.785_791del 1 +17 31259095 G C NC_000017.11:g.31259095G>C NM_000267.3:c.4333G>C 1 +X 154030810 T TA NC_000023.11:g.154030811dup NM_001110792.2:c.1053dup 1 +13 51974987 T C NC_000013.11:g.51974987T>C NM_000053.4:c.233A>G 1 +15 89776908 AGGGGCAGG A NC_000015.10:g.89776911_89776918del NM_001039958.1:c.554_561delGGCAGGGG 1 +7 117504365 T C NC_000007.14:g.117504365T>C NM_000492.4:c.164+2T>C 3 +16 53649101 T C NC_000016.10:g.53649101T>C NM_015272.5:c.2167A>G 1 +3 186785958 CGAAAT C NC_000003.12:g.186785960AAATG[1] NM_001967.4:c.431_435del 1 +5 90629418 G T NC_000005.10:g.90629418G>T NM_032119.3:c.1718G>T 4 +7 5977588 C T NC_000007.14:g.5977588C>T NM_000535.7:c.2445G>A 2 +13 110166257 G A NC_000013.11:g.110166257G>A NM_001845.6:c.3996C>T 1 +19 48297312 G A NC_000019.10:g.48297312G>A NM_144577.3:c.1677C>T 2 +16 28900625 G A NC_000016.10:g.28900625G>A NM_173201.4:c.1809G>A 1 +12 51913244 C T NC_000012.12:g.51913244C>T NM_000020.3:c.207C>T 2 +14 95107692 A G NC_000014.9:g.95107692A>G NM_177438.3:c.2720T>C 2 +2 178609740 G A NC_000002.12:g.178609740G>A NM_003319.4:c.24488C>T 1 +19 11041430 G A NC_000019.10:g.11041430G>A NM_001128849.1:c.4390G>A 2 +7 128848826 C T NC_000007.14:g.128848826C>T NM_001458.4:c.4771C>T 4 +2 1943100 C T NC_000002.12:g.1943100C>T NM_015025.2:c.387G>A 1 +2 165095539 A C NC_000002.12:g.165095539A>C NM_006922.3:c.4403T>G 1 +1 7984986 A G NC_000001.11:g.7984986A>G NM_007262.4:c.502A>G 1 +9 99149272 G A NC_000009.12:g.99149272G>A NM_004612.3:c.1479G>A 1 +1 181715381 C A NC_000001.11:g.181715381C>A NM_000721.3:c.1215C>A 2 +3 49530951 A G NC_000003.12:g.49530951A>G NM_004393.6:c.440A>G 1 +11 17395689 C T NC_000011.10:g.17395689C>T NM_000352.6:c.4228G>A 1 +15 48412654 G A NC_000015.10:g.48412654G>A NM_000138.4:c.8141C>T 1 +16 2076528 C A NC_000016.10:g.2076528C>A NM_000548.5:c.2780C>A 3 +9 16419194 T C NC_000009.12:g.16419194T>C NM_017637.5:c.3095A>G 1 +8 11708487 G A NC_000008.11:g.11708487G>A NM_002052.3:c.175G>A 1 +15 20534995 C A NC_000015.10:g.20534995C>A NM_001145004.1:c.1517G>T 1 +7 116795952 T C NC_000007.14:g.116795952T>C NM_001127500.1:c.4055T>C 1 +16 10682083 G A NC_000016.10:g.10682083G>A NM_144674.1:c.773C>T 1 +17 42103632 T C NC_000017.11:g.42103632T>C NM_024119.2:c.1730A>G 1 +2 238427110 G A NC_000002.12:g.238427110G>A NM_001040445.1:c.40G>A 1 +20 32435682 A G NC_000020.11:g.32435682A>G NM_015338.5:c.2970A>G 1 +19 46755906 C G NC_000019.10:g.46755906C>G NM_024301.4:c.456C>G 6 +4 990302 C G NC_000004.12:g.990302C>G NM_213613.2:c.637G>C 1 +3 39130990 A T NC_000003.12:g.39130990A>T NM_001366900.1:c.2459-2A>T 1 +11 112094890 T G NC_000011.10:g.112094890T>G NM_003002.3:c.400T>G 1 +10 132785714 G GC NC_000010.11:g.132785721dup NM_177400.3:c.234dup 2 +12 6592031 C T NC_000012.12:g.6592031C>T NM_001273.2:c.2975G>A 1 +17 43074362 C T NC_000017.11:g.43074362C>T NM_007294.4:c.4644G>A 1 +17 43124077 C A NC_000017.11:g.43124077C>A NM_007294.3:c.20G>T 1 +2 178547242 C T NC_000002.12:g.178547242C>T NM_001267550.2:c.94283G>A 2 +17 43091755 T G NC_000017.11:g.43091755T>G NM_007294.3:c.3776A>C 6 +X 154380001 TTACTC T NC_000023.11:g.154380005_154380009del NM_000117.2:c.251_255del5 1 +11 19182687 CTTG C NC_000011.10:g.19182689TGT[1] NM_003476.4:c.565_567delCAA 1 +19 55401545 C T NC_000019.10:g.55401545C>T NM_014501.2:c.560G>A 1 +17 34156387 A G NC_000017.11:g.34156387A>G NM_001094.4:c.146T>C 1 From 20469bd712ea3cc442c0ed3380a65bf0195277e5 Mon Sep 17 00:00:00 2001 From: Dave Lawrence Date: Mon, 17 Aug 2026 21:52:49 +0000 Subject: [PATCH 3/3] Paper: LOVD syntax checker head-to-head on the injection corpus - #112 Run the exact inject_and_clean.py corpus (3,419 cases) through the LOVD HGVS syntax checker (v1.2.2, local PHP CLI) and clean_hgvs(), scored identically (gene-annotation-insensitive exact match to the canonical target). Weighted by the production error mix LOVD's top-ranked correction recovers 71.7% vs 100% (by construction) for clean_hgvs(); neither tool alters a valid input. New lovd_comparison Snakefile rule with frozen-constants fallback, facts CSV, and Methods/Discussion prose. --- claude/20260817_lovd_head_to_head_plan.md | 71 +++++ paper/Snakefile | 26 +- paper/discussion.md | 21 ++ paper/empirical_results/lovd_comparison.csv | 2 + paper/methods.md | 14 +- paper/scripts/lovd_head_to_head.py | 309 ++++++++++++++++++++ 6 files changed, 441 insertions(+), 2 deletions(-) create mode 100644 claude/20260817_lovd_head_to_head_plan.md create mode 100644 paper/empirical_results/lovd_comparison.csv create mode 100644 paper/scripts/lovd_head_to_head.py diff --git a/claude/20260817_lovd_head_to_head_plan.md b/claude/20260817_lovd_head_to_head_plan.md new file mode 100644 index 0000000..8047f69 --- /dev/null +++ b/claude/20260817_lovd_head_to_head_plan.md @@ -0,0 +1,71 @@ +# LOVD syntax checker head-to-head plan + +Follow-up to item 3 of the (now deleted) 20260808 paper feedback plan. The related-work +paragraph in `paper/discussion.md` now cites the LOVD HGVS syntax checker +(github.com/LOVDnl/HGVS-syntax-checker) as a repair tool. A Bioinformatics reviewer is +likely to ask how `clean_hgvs()` compares, so run the comparison pre-emptively: score both +tools on the same corrupted-HGVS corpus and report recovery head-to-head. + +## Feasibility facts (assessed 2026-08-17) + +- The checker is a self-contained PHP library, PHP >= 7.4 with json and zlib extensions, + runnable fully offline from the CLI: `php -f HGVS.php "NM_004006.3:c.157C>T"` emits JSON + with `valid` and `corrected_values` (corrections ranked by likelihood, each with a + confidence score). +- Installable via composer as `lovd/hgvs-syntax-checker`, or plain clone of the repo. +- Public REST fallback: api.lovd.nl checkHGVS v2 (single or batch), rate-limited to about + 5 variants/s or 1 batch/s, no registration. Local CLI is minutes for the whole corpus, + REST would be ~12 minutes; use the local CLI. +- Estimated effort: about half a day including the PHP install and a subprocess wrapper. + +## Corpus + +`paper/scripts/inject_and_clean.py` generates ~3,400 corrupted strings by injecting known +error classes into committed ClinVar c.HGVS (all public data, safe to use anywhere). Each +corrupted string has a known canonical target, so scoring is exact-match recovery. Reuse +its injection machinery and categories; do not invent a new corpus. + +## Method + +1. Install PHP + the checker locally (pin the checker version and record it in Methods for + reproducibility; a composer.lock or pinned git SHA in the script docstring is enough). +2. New script `paper/scripts/lovd_head_to_head.py`: subprocess wrapper around the PHP CLI, + feeding each corrupted string, parsing the JSON, and scoring. Batch via the CLI's own + batching if it has one, otherwise one call per string is fine at local speed. +3. Scoring must be decided and documented up front, then applied identically to both tools: + - Primary metric: top-ranked LOVD correction equals the canonical target, versus + `clean_hgvs()` output equals the canonical target. + - Secondary (report, do not headline): target appears anywhere in LOVD's ranked + `corrected_values` list. +4. Fairness stratification, the key design point: LOVD's checker is deliberately + sequence-agnostic, so injection categories that require transcript data (anything a + data-provider-aware fix handles, eg accession-prefix restoration, version work) are + structurally out of its scope. Split results into two strata: string-repairable + (both tools eligible) and data-requiring (cdot only). Report the head-to-head on the + first stratum and the second as cdot-only capability. Do not present a blended number + that makes LOVD look artificially weak. +5. Also record LOVD's behaviour on the uncorrupted originals (false-correction rate), since + the paper claims `clean_hgvs()` guarantees no regressions; the comparison is only + meaningful with the same check applied to LOVD. + +## Deliverables + +- `paper/scripts/lovd_head_to_head.py` (plus whatever tiny install notes it needs in its + docstring; the PHP install itself stays out of the repo). +- Facts CSV under `paper/empirical_results/` (eg `lovd_comparison.csv`) following the + existing pattern, wired into `paper/Snakefile` as a frozen-constants rule with a + provenance docstring, matching how the other measured facts are handled. +- One short paragraph in `paper/discussion.md` extending the existing related-work + paragraph with the measured comparison (or a sentence there plus a supplementary table + if the strata need a table). Methods gets the scoring and stratification rules. +- No CHANGELOG entry (paper/analysis only). No em-dashes in prose. Nothing from + `../cdot_private` (the injection corpus is public ClinVar, so this is naturally safe). + +## Caveats + +- If the checker cannot be installed offline on this machine (no PHP, composer blocked), + fall back to the REST API with the rate limit respected, and note the service date in + Methods instead of a version pin. +- LOVD may legitimately return multiple plausible corrections where the injection is + ambiguous; that is what the secondary metric captures. Do not count ambiguity as failure + in prose without saying the top-1 rule caused it. diff --git a/paper/Snakefile b/paper/Snakefile index 99f5fdd..8eca596 100644 --- a/paper/Snakefile +++ b/paper/Snakefile @@ -66,7 +66,7 @@ EMPIRICAL = "paper/empirical_results" GEN_FACTS = "output/facts" FACT_FILES = ["literature.csv", "coverage.csv", "benchmark.csv", "clinvar.csv", "clinvar_submitted.csv", "clinvar_submitted_residual.csv", - "cleaning.csv", "sources.csv", "historical.csv", + "cleaning.csv", "lovd_comparison.csv", "sources.csv", "historical.csv", "version_stability.csv", "positional_drift.csv"] PAPER_SOURCES = [ @@ -429,6 +429,30 @@ rule cleaning: "{PYTHON} paper/scripts/inject_and_clean.py" +rule lovd_comparison: + """Tier-1 head-to-head: clean_hgvs() vs the LOVD HGVS syntax checker (issue #112). + + Runs the exact injection corpus of the `cleaning` rule through the LOVD + checker (github.com/LOVDnl/HGVS-syntax-checker, pinned v1.2.2, local PHP + CLI) and through clean_hgvs(), scored identically; see + paper/scripts/lovd_head_to_head.py for the scoring rules. Needs a local + clone of the checker: pass --config lovd_checker=/path/to/HGVS.php (PHP + >= 7.4 CLI required). Without it, the frozen measured constants from + paper/empirical_results/ are copied (measured 2026-08-17, checker v1.2.2, + commit 4cd074ba9cbd, PHP 8.5). + """ + output: "output/facts/lovd_comparison.csv" + run: + import os, shutil + os.makedirs("output/facts", exist_ok=True) + checker = config.get("lovd_checker", "") + if checker and Path(checker).exists(): + shell(f"{PYTHON} paper/scripts/lovd_head_to_head.py --checker {checker}") + else: + print("No lovd_checker - copying the frozen measured LOVD comparison facts") + shutil.copy(f"{EMPIRICAL}/lovd_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 971d352..98a53a2 100644 --- a/paper/discussion.md +++ b/paper/discussion.md @@ -49,6 +49,27 @@ anything that consumes the string downstream benefits; every change is returned `HGVSFix` the caller can audit; and it guarantees no regressions, never breaking a description that already parsed. +The LOVD checker is the closest comparator, being the one tool in this group that runs +offline without a reference sequence, so we ran it head-to-head with `clean_hgvs()` on +the reproducible injection corpus ({{ lovd_comparison.n_cases | commas }} corrupted +strings, checker {{ lovd_comparison.lovd_version }}; scoring in Methods). The corpus +injects the error classes `clean_hgvs()` targets, so `clean_hgvs()` recovers every case +by construction and the comparison measures how much of that territory a general syntax +checker also covers, not the overall quality of either tool. Weighted by the production +error mix, LOVD's top-ranked correction restored +{{ lovd_comparison.lovd_top1_weighted_pct | dp(0) }}% of cases +({{ lovd_comparison.lovd_top1_pct | dp(0) }}% unweighted across categories), matching +`clean_hgvs()` on the common single-token errors: whitespace, letter case, separator +typos and trailing protein annotations. The rest splits into inputs that must be treated +as free text rather than a malformed variant description (surrounding quotes, unbalanced +brackets, leading assembly text), which sit outside the checker's intended input, and +repairs applied to only one accession family (it restores a swapped gene and transcript +for RefSeq but not Ensembl accessions, and re-cases a lowercase Ensembl accession but +not a RefSeq one). Neither tool altered any valid input, though LOVD flags intronic +descriptions on a transcript reference +({{ lovd_comparison.lovd_flagged_invalid_pct | dp(0) }}% of the valid originals) as +requiring a genomic reference, a deliberate design position rather than a defect. + Beyond HGVS resolution, the JSON representation is useful in its own right. It parses far faster than the GTF/GFF files it is built from and loads trivially over HTTP, so cdot doubles as a lightweight, queryable gene/transcript reference. We publish the per-release diff --git a/paper/empirical_results/lovd_comparison.csv b/paper/empirical_results/lovd_comparison.csv new file mode 100644 index 0000000..4f0e2a4 --- /dev/null +++ b/paper/empirical_results/lovd_comparison.csv @@ -0,0 +1,2 @@ +lovd_version,n_cases,cdot_pct,lovd_top1_pct,lovd_anyrank_pct,cdot_weighted_pct,lovd_top1_weighted_pct,lovd_anyrank_weighted_pct,originals_n,lovd_false_corrections,lovd_flagged_invalid,lovd_flagged_invalid_pct,cdot_false_corrections +v1.2.2,3419,100.0,33.0,33.0,100.0,71.7,71.7,998,0,266,26.7,0 diff --git a/paper/methods.md b/paper/methods.md index bce594c..61b8f19 100644 --- a/paper/methods.md +++ b/paper/methods.md @@ -224,7 +224,19 @@ a pluggable provider (local JSON, REST, or UTA) and reports resolution rate, rec cleaning and version fallback, and speed; the ClinVar pair set is built by `build_clinvar_pairs.py`. Cleaning is evaluated on a production query corpus and, as a reproducible control, with `inject_and_clean.py`, which injects each fix category into -clean ClinVar strings. Version-fallback safety is measured by `compute_version_stability.py` +clean ClinVar strings. `lovd_head_to_head.py` runs the same injected cases (same seed and +per-category caps) through both `clean_hgvs()` and the LOVD HGVS syntax checker +[@LovdHgvsChecker] ({{ lovd_comparison.lovd_version }}, run locally as a PHP CLI), scored +with one rule: a case is recovered when the tool's output (for LOVD, its top-ranked +suggested correction) exactly matches the known canonical target. The two ecosystems +canonicalise a parenthesised gene symbol in opposite directions (biocommons keeps +`NM_x.y(GENE):c.`, LOVD removes the symbol), so the comparison ignores that annotation on +both sides. A secondary metric accepting the target anywhere in LOVD's ranked correction +list gave identical results, so top-1 ranking cost LOVD nothing. Every injected category +is string-repairable by construction (`inject_and_clean.py` does not inject errors whose +repair needs transcript data, such as a missing accession prefix), so both tools are +eligible on every case; the same false-correction check (does the tool alter a valid +input) is applied to both over the uncorrupted originals. Version-fallback safety is measured by `compute_version_stability.py` on GRCh38, using a seeded {{ version_stability.sample_n | commas }}-accession sample drawn from accessions cdot holds at two or more versions (the only accessions where a version bump can be assessed). The same run bins preserved coding bases by relative CDS position diff --git a/paper/scripts/lovd_head_to_head.py b/paper/scripts/lovd_head_to_head.py new file mode 100644 index 0000000..170d377 --- /dev/null +++ b/paper/scripts/lovd_head_to_head.py @@ -0,0 +1,309 @@ +#!/usr/bin/env python3 +""" +Head-to-head: cdot ``clean_hgvs()`` vs the LOVD HGVS syntax checker (issue #112). + +Runs the exact injection corpus from ``inject_and_clean.py`` (same seed, same +per-class caps, same injectors, so the case list is identical) through both +tools and scores them with the same rule. + +LOVD HGVS syntax checker +------------------------ +github.com/LOVDnl/HGVS-syntax-checker: a deliberately sequence-agnostic PHP +library that validates HGVS syntax and suggests corrections for invalid +descriptions, ranked by a confidence score. Run locally from the CLI:: + + git clone https://github.com/LOVDnl/HGVS-syntax-checker + (cd HGVS-syntax-checker && git checkout v1.2.2) + php -f HGVS-syntax-checker/HGVS.php "NM_004006.3:c.157C>T" # emits JSON + +Version pinned for the paper: **v1.2.2** (commit 4cd074ba9cbd, 2026-06-15), +PHP 8.5 CLI, fully offline. Pass the path to ``HGVS.php`` via ``--checker``. + +Scoring (decided up front, applied identically to both tools) +------------------------------------------------------------- +Each injected case has a known canonical target (``expected``). A tool +recovers a case iff its output string equals the target, compared +*gene-annotation-insensitively*: a parenthesised gene symbol between the +accession and the colon is stripped from both sides before comparing, because +the two ecosystems canonicalise that form in opposite directions (biocommons +keeps ``NM_x.y(GENE):c.``, LOVD corrects it to ``NM_x.y:c.``) and both denote +the same variant. Concretely: + +- **cdot**: ``clean_hgvs(perturbed)`` == target. +- **LOVD top-1** (primary): the checker's highest-confidence entry in + ``corrected_values`` == target. +- **LOVD any-rank** (secondary): target appears anywhere in the ranked + ``corrected_values`` list. + +False-correction check on the uncorrupted originals: a tool "falsely corrects" +a valid input if its output differs from the input (same insensitive +comparison). For LOVD this also counts inputs it flags invalid. + +Stratification note: the fairness plan called for splitting string-repairable +vs data-requiring categories, but the injection corpus contains *only* +string-repairable categories by construction (``inject_and_clean.py`` already +excludes ops that need a data provider, eg accession-prefix restoration), so +both tools are eligible on every case and no split is needed. Per-category +results are still emitted so scope disputes stay visible. + +Usage:: + + python paper/scripts/lovd_head_to_head.py --checker /path/to/HGVS.php +""" + +import argparse +import csv +import json +import random +import re +import subprocess +import sys +import time +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 (injectors, seed, sample, weights) + +FACTS_DIR = iac.FACTS_DIR +LOVD_VERSION = "v1.2.2" + +# Strip a parenthesised gene symbol between accession and colon: +# "NM_000059.4(BRCA2):c.316G>A" -> "NM_000059.4:c.316G>A" +_GENE_ANNOT = re.compile(r"^([A-Za-z]{2,4}_?\d+\.\d+)\(([^()]+)\):") + + +def norm(s): + return _GENE_ANNOT.sub(r"\1:", s) + + +# --------------------------------------------------------------------------- +# Case generation: mirrors inject_and_clean.run() exactly (same seeds, same +# skip rules, same caps) so the attempted case list is identical. +# --------------------------------------------------------------------------- + +def build_cases(parser, pool): + cases = [] # (category, orig, perturbed, expected) + for idx, (name, injector) in enumerate(iac.INJECTORS): + order = list(pool) + rng_c = random.Random(iac.SEED + idx) + rng_c.shuffle(order) + attempted = 0 + gi = 0 + for orig in order: + res = injector(orig, iac.GENE_POOL[gi % len(iac.GENE_POOL)]) + if res is None: + continue + if attempted >= iac.TARGET_PER_CLASS: + continue + gi += 1 + perturbed, expected = res + if not iac.parses(parser, expected): + continue + if perturbed == expected: + continue + attempted += 1 + cases.append((name, orig, perturbed, expected)) + return cases + + +# --------------------------------------------------------------------------- +# LOVD CLI driver +# --------------------------------------------------------------------------- + +def run_lovd(checker, strings, chunk=500): + """Run the PHP CLI over strings (batched argv), return {input: result}.""" + results = {} + todo = [s for s in dict.fromkeys(strings)] # dedupe, keep order + for i in range(0, len(todo), chunk): + batch = todo[i:i + chunk] + proc = subprocess.run( + ["php", "-f", str(checker), *batch], + capture_output=True, text=True, check=True) + for item in json.loads(proc.stdout): + results[item["input"]] = item + missing = [s for s in todo if s not in results] + if missing: + raise RuntimeError(f"LOVD checker returned no result for {len(missing)} " + f"inputs, eg {missing[0]!r}") + return results + + +def lovd_ranked(result): + """Ranked correction list (best first) from a checker result.""" + cv = result.get("corrected_values") or {} + if not isinstance(cv, dict): # empty PHP array serialises as [] + return [] + return sorted(cv, key=lambda k: -cv[k]) + + +# --------------------------------------------------------------------------- +# Scoring +# --------------------------------------------------------------------------- + +def score(cases, lovd_results): + per_class = {} + for name, _ in iac.INJECTORS: + per_class[name] = { + "weight": iac.REAL_RESCUE_OP_COUNTS[name], + "n_attempted": 0, "cdot_recovered": 0, + "lovd_top1": 0, "lovd_anyrank": 0, + } + for name, orig, perturbed, expected in cases: + c = per_class[name] + c["n_attempted"] += 1 + target = norm(expected) + cleaned, _ = clean_hgvs(perturbed) + if norm(cleaned) == target: + c["cdot_recovered"] += 1 + ranked = [norm(v) for v in lovd_ranked(lovd_results[perturbed])] + if ranked and ranked[0] == target: + c["lovd_top1"] += 1 + if target in ranked: + c["lovd_anyrank"] += 1 + + for c in per_class.values(): + n = c["n_attempted"] + for k in ("cdot_recovered", "lovd_top1", "lovd_anyrank"): + c[k + "_pct"] = round(100.0 * c[k] / n, 1) if n else None + return per_class + + +def weighted(per_class, key): + wsum = wnum = 0.0 + for c in per_class.values(): + if c["n_attempted"]: + wsum += c["weight"] + wnum += c["weight"] * c[key + "_pct"] + return round(wnum / wsum, 1) if wsum else None + + +def false_corrections(pool, lovd_results): + lovd_fc = cdot_fc = lovd_invalid = 0 + examples = [] + for orig in pool: + target = norm(orig) + cleaned, _ = clean_hgvs(orig) + if norm(cleaned) != target: + cdot_fc += 1 + r = lovd_results[orig] + ranked = [norm(v) for v in lovd_ranked(r)] + changed = not ranked or ranked[0] != target + if not r.get("valid"): + lovd_invalid += 1 + if changed: + lovd_fc += 1 + if len(examples) < 10: + examples.append({"input": orig, "valid": r.get("valid"), + "top": ranked[0] if ranked else None}) + return lovd_fc, lovd_invalid, cdot_fc, examples + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("--checker", required=True, type=Path, + help="path to the LOVD checker's HGVS.php (pinned v1.2.2)") + ap.add_argument("--limit", type=int, default=None, + help="cap cases per category (smoke test only; not for facts)") + args = ap.parse_args() + if not args.checker.exists(): + sys.exit(f"checker not found: {args.checker}") + + 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] + + t0 = time.time() + lovd_results = run_lovd(args.checker, [c[2] for c in cases] + list(pool)) + lovd_seconds = round(time.time() - t0, 1) + + per_class = score(cases, lovd_results) + n = sum(c["n_attempted"] for c in per_class.values()) + totals = {k: sum(c[k] for c in per_class.values()) + for k in ("cdot_recovered", "lovd_top1", "lovd_anyrank")} + lovd_fc, lovd_invalid, cdot_fc, fc_examples = false_corrections(pool, lovd_results) + + facts = { + "issue": "SACGF/cdot#112", + "tier": 1, + "description": ( + "Head-to-head on the reproducible injection corpus: cdot clean_hgvs() " + "vs the LOVD HGVS syntax checker (local PHP CLI, pinned " + f"{LOVD_VERSION}). Same cases, same gene-annotation-insensitive " + "exact-match scoring; see paper/scripts/lovd_head_to_head.py." + ), + "lovd_version": LOVD_VERSION, + "n_cases": n, + "sample_size": len(pool), + "seed": iac.SEED, + "cdot_recovered": totals["cdot_recovered"], + "cdot_pct": round(100.0 * totals["cdot_recovered"] / n, 1), + "lovd_top1": totals["lovd_top1"], + "lovd_top1_pct": round(100.0 * totals["lovd_top1"] / n, 1), + "lovd_anyrank": totals["lovd_anyrank"], + "lovd_anyrank_pct": round(100.0 * totals["lovd_anyrank"] / n, 1), + "cdot_weighted_pct": weighted(per_class, "cdot_recovered"), + "lovd_top1_weighted_pct": weighted(per_class, "lovd_top1"), + "lovd_anyrank_weighted_pct": weighted(per_class, "lovd_anyrank"), + "originals_n": len(pool), + "lovd_false_corrections": lovd_fc, + "lovd_flagged_invalid": lovd_invalid, + "lovd_flagged_invalid_pct": round(100.0 * lovd_invalid / len(pool), 1), + "cdot_false_corrections": cdot_fc, + "lovd_seconds": lovd_seconds, + "per_class": per_class, + "false_correction_examples": fc_examples, + } + + FACTS_DIR.mkdir(parents=True, exist_ok=True) + json_path = FACTS_DIR / "lovd_comparison.json" + json_path.write_text(json.dumps(facts, indent=2) + "\n") + csv_path = FACTS_DIR / "lovd_comparison.csv" + row = {k: facts[k] for k in ( + "lovd_version", "n_cases", "cdot_pct", "lovd_top1_pct", + "lovd_anyrank_pct", "cdot_weighted_pct", "lovd_top1_weighted_pct", + "lovd_anyrank_weighted_pct", "originals_n", "lovd_false_corrections", + "lovd_flagged_invalid", "lovd_flagged_invalid_pct", + "cdot_false_corrections")} + with open(csv_path, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=list(row)) + w.writeheader() + w.writerow(row) + + print(f"\nLOVD head-to-head ({LOVD_VERSION}, {n} injected cases, " + f"LOVD wall time {lovd_seconds}s)") + print(f" {'category':<34} {'n':>4} {'cdot%':>7} {'lovd@1%':>8} {'lovd@any%':>9}") + for name, c in per_class.items(): + if not c["n_attempted"]: + continue + print(f" {name:<34} {c['n_attempted']:>4} {c['cdot_recovered_pct']:>7} " + f"{c['lovd_top1_pct']:>8} {c['lovd_anyrank_pct']:>9}") + print(f" {'-'*66}") + print(f" overall : cdot {facts['cdot_pct']}% lovd top-1 {facts['lovd_top1_pct']}% " + f"lovd any-rank {facts['lovd_anyrank_pct']}%") + print(f" weighted: cdot {facts['cdot_weighted_pct']}% " + f"lovd top-1 {facts['lovd_top1_weighted_pct']}% " + f"lovd any-rank {facts['lovd_anyrank_weighted_pct']}%") + print(f" originals ({len(pool)}): lovd false-corrections {lovd_fc} " + f"(flagged invalid {lovd_invalid}), cdot false-corrections {cdot_fc}") + if fc_examples: + print(" lovd false-correction examples:") + for e in fc_examples: + print(f" - {e['input']!r} -> {e['top']!r} (valid={e['valid']})") + print(f"\nWrote: {json_path}") + print(f"Wrote: {csv_path}") + + +if __name__ == "__main__": + main()