Address presentation feedback: submitted-SCV corpus, benchmark rework, cleaning rules, positional drift

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)

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,
 )

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),

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):