Skip to content

Enable use of RefGetStore for sequence and alias data - #626

Draft
theferrit32 wants to merge 40 commits into
ga4gh:v3from
theferrit32:kf/dataproxy-backend-refactor
Draft

theferrit32 wants to merge 40 commits into
ga4gh:v3from
theferrit32:kf/dataproxy-backend-refactor

Conversation

@theferrit32

@theferrit32 theferrit32 commented Apr 16, 2026 •

Copy link
Copy Markdown
Contributor

NOTE: This draft PR contains work-in-progress and analysis scripts and output files, and should not be merged as-is. Files from the fork branch will be cherry picked into a later clean PR.

The purpose is to provide a view into proposed source changes, and challenges being worked through.

theferrit32 and others added 13 commits February 12, 2026 15:26
…o backend

Extract SeqRepo-specific implementations (SeqRepoDataProxy,
SeqRepoRESTDataProxy) from dataproxy.py into dataproxy_seqrepo.py.
Replace monolithic _DataProxy ABC with composable SequenceDataProxy
and AliasDataProxy interfaces, combined as DataProxy. Replace
get_metadata() with get_sequence_length() and get_aliases(). Add
provider registry for extensible backend support.
Add CompositeDataProxy that composes a SequenceDataProxy and
AliasDataProxy into a full DataProxy, with automatic identifier
resolution fallback via alias lookup. Add RefgetSequenceDataProxy
backed by gtars RefgetStore for high-performance sequence retrieval
by GA4GH digest. Add gtars as optional "refget" dependency group.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
…e archive

Build and verify a gtars RefgetStore archive with the alias namespaces
vrs-python needs (GRCh38, GRCh38.p14, GRCh37, GRCh37.p13, refseq, insdc),
as an eventual replacement for seqrepo alias mapping. Independent of
bioutils/seqrepo — reads NCBI assembly reports directly.

- build_store.py: ingests NCBI _genomic.fna.gz files via RefgetStore and
  writes per-sequence and per-collection aliases from the matching
  assembly_report.txt. Rows with RefSeq-Accn=na are tracked as
  rows_skipped_non_refseq (GenBank-only contigs omitted from RefSeq's
  _genomic.fna.gz have no sequence data to alias to).
- assemblies.toml: declarative list of assembly FASTA/report URLs.
- verify_store.py: standalone checks that every RefSeq row in every
  assembly report resolves under refseq:/<assembly>:/insdc:, asserts
  ground-truth digests match seqrepo 2024-12-20, confirms cross-patch
  and cross-build invariants (chr1 stable within builds, chrM stable
  across builds), and asserts no legacy namespaces leaked in.
- README.md: how-to and config reference.
- .gitignore: exclude misc/refgetstore/{store,downloads}/.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Add RefgetAliasDataProxy, backed by a gtars RefgetStore's sequence alias
index, so the refget backend can answer the full AliasDataProxy interface
(not just sequence retrieval by digest). Wire up `refget+file://` so
create_dataproxy returns a CompositeDataProxy of the existing
RefgetSequenceDataProxy and the new alias proxy over a single shared
store instance.

- dataproxy_refget.py: new RefgetAliasDataProxy class with
  _resolve_to_digest handling every input form vrs-python callers
  actually construct (ga4gh:SQ.<d>, sha512t24u:<d>, SQ.<d>,
  namespace:alias, and bare accessions via bioutils coerce_namespace).
  _extract_digest gains a sha512t24u: branch and is reused by both the
  sequence and alias proxies. _create_refget_dataproxy now returns a
  CompositeDataProxy(sequence_proxy, alias_proxy).
- dataproxy.py: create_dataproxy recognizes the `refget+file://` URI
  scheme and delegates to the refget factory. Also updates the docstring
  with the new accepted scheme. No changes to the abstract interfaces or
  CompositeDataProxy beyond the existing get_aliases forwarding.
- misc/refgetstore/compare_against_seqrepo.py: three-phase harness
  verifying the refget proxy against the seqrepo proxy over a curated
  set of canonical identifiers. Phase 1 asserts functional equivalence
  (derive_refget_accession, translate ns=ga4gh/refseq, get_sequence_length)
  on every identifier both backends natively support. Phase 2 covers
  refget-only extensions (GRCh38.p14, insdc:, sha512t24u:). Phase 3
  checks that both proxies emit refseq:<accn> and ga4gh:SQ.<digest>
  entries from get_aliases, the only aliases vrs-python callers
  actually consume. 74/0 EQUIVALENT against seqrepo 2024-12-20.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Keeps the spec-documented DataProxy.get_metadata method as a concrete
default that composes get_sequence_length + get_aliases. Omits the
alphabet key — SeqRepo's per-sequence character set doesn't have a
faithful equivalent in the refget backend and no code in vrs-python
reads it. SeqRepoRESTDataProxy overrides to preserve its single
round-trip behavior.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
RefgetStore.open_local() reads sequences.rgsi and populates the in-memory
store with a Stub record for every sequence, so metadata lookups and alias
resolution work immediately. But byte-level reads (get_substring) need
the per-digest payload to be fetched via load_sequence(digest) first, which
was never being called on the dataproxy path — so get_sequence() raised
on every on-disk store.

Add a small _ensure_sequence_loaded() helper that calls load_sequence(digest)
exactly once per digest, memoized in a set. Nothing is loaded until a caller
actually asks for a sequence, and the working set is bounded by what the
calling application actually touches — not by the total store size.
Extend build_store.py with a [[seqset]] config kind for flat FASTAs where
the header name is the accession (NCBI RefSeq mRNA_Prot shards). Each shard
is ingested as its own collection; its header names become refseq:<accn>
aliases on the resulting digests. Adds ~470k aliases covering NM_/NR_/XM_/
XR_ (16 RNA shards, 273k entries) and NP_/XP_ (16 protein shards, 198k
entries), on top of the existing 4 assembly collections.

Critical perf fix: gtars.RefgetStore.add_sequence_alias rewrites the full
namespace TSV on every call (O(n) per insert → O(n^2) to populate), which
made per-alias inserts take ~7.5 ms each and blew up to hours on a 21k-seq
shard. Rewire process_seqset to accumulate (alias, digest) pairs across
shards into a pending dict, write them to a temp TSV, and hand the whole
file to store.load_sequence_aliases() which merges and persists the
namespace in one pass — ~5 µs per alias, 1500x speedup. The full seqset
build now finishes in well under a minute after the FASTAs are cached.

Other changes:
  - assemblies.toml: add refseq_human_rna + refseq_human_protein entries
    with shard_range=[1,16] (the max NCBI actually publishes; 17–32 404
    and are reported as warnings rather than failures).
  - download.sh: pre-fetcher that curls FASTAs and assembly reports into
    the loader's cache dir with the filenames build_store.py expects.
    Supports resume via curl -C -.
  - verify_store.py: relax the "n_collections == 4" inventory check to
    ">= 4" so the assembly-only and assembly+seqsets builds both pass.
  - README.md: document the [[seqset]] kind, new CLI flags, and how to
    build or skip each category.

Validated end-to-end with 41/41 verify_store checks and 74/74
compare_against_seqrepo checks passing post-build.
Drives ga4gh.vrs.extras.annotator.VcfAnnotator over the same input VCF
with two dataproxy configurations — SeqRepoDataProxy (unbounded fd cache)
and CompositeDataProxy(RefgetSequenceDataProxy, RefgetAliasDataProxy) —
and compares the resulting VRS_Allele_IDs / VRS_Starts / VRS_Ends /
VRS_States / VRS_Lengths / VRS_Error INFO fields for exact agreement.

Both backends run single-threaded (pysam.VariantFile defaults to
threads=1, neither refget nor seqrepo uses worker pools). Seqrepo is
always constructed with fd_cache_size=None so the baseline reflects
steady-state performance, not per-record file-open cost. Suggests
exporting a larger SEQREPO_LRU_CACHE_MAXSIZE if the workload has more
than ~1M distinct sequence fetches.

Exits non-zero on any INFO-field mismatch. On a ~1M-record gnomAD chr1
VCF (GRCh38): 1_000_000/1_000_000 VRS_Allele_IDs match, 1_000_000/1_000_000
full-INFO match, refget_native 120.43s vs seqrepo 137.10s (1.14x).
Complements vcf_equivalence.py by driving the translator path directly
rather than the VCF annotator. Iterates every record of a large VCF,
constructs gnomad-style chrom-pos-ref-alt expressions, and feeds each
through ga4gh.vrs.extras.translator.AlleleTranslator under two backends:
SeqRepoDataProxy (unbounded fd cache) and CompositeDataProxy(refget).

Exercises the same dataproxy call sites HGVS/SPDI import paths hit, not
just the annotator's gnomad-to-allele shortcut. Reports match rate, first
N mismatches with coordinates, and end-to-end wall time per backend.
Writes an optional JSON report via --report. Single-threaded throughout.

Exits non-zero on any Allele ID mismatch. On ~1M-record gnomAD chr1 VCF
(GRCh38): 1_000_000/1_000_000 match, refget_native 53.97s vs seqrepo
62.01s (1.15x).
…parison

Add three Ensembl release-113 seqsets (cdna, ncrna, pep) to
assemblies.toml, bringing ~534k ensembl-namespace aliases into the
store. Update verify_store.py with ground-truth checks for VHL,
BRCA1, XIST, and TP53 across all three biotypes. Update README to
document Ensembl coverage and the seqset schema.

Include a detailed comparison report against seqrepo 2024-12-20:
- RefSeq: 464,830 common aliases, zero digest mismatches
- Ensembl: 322,286 common (83%), 116 ENSP digest mismatches (0.1%)
- ~212k refget-only Ensembl aliases, mostly ncRNA coverage

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Root-caused the 116 Ensembl digest mismatches between seqrepo and
refgetstore. 115/116 are caused by bioutils.normalize_sequence()
stripping '*' (stop codon) characters before computing sha512t24u,
while gtars hashes raw FASTA bytes including '*'. The underlying
sequence bytes are identical in both stores. The remaining 1/116
(ENST00000668831.1) is a genuine sequence divergence across Ensembl
releases.

All 115 affected ENSPs are Ensembl-only annotations (no RefSeq
counterpart) from pseudogenes, IG D-segments, KIR receptors, and
olfactory receptors. NCBI RefSeq protein FASTAs contain zero '*'
characters, which is why this was invisible until Ensembl protein
data was loaded.

New files:
- diagnose_ensembl_mismatches.py: 4-phase diagnostic script
- diagnosis_report.json: full JSON report for all 116 mismatches
- ensembl_known_divergent.txt: TSV of affected accessions
- ensembl_translator_equivalence.py: SPDI-based Ensembl equivalence
  test (25/25 pass)
- Two investigation reports documenting findings and spec analysis

Updated files:
- verify_store.py: add [I] known-divergent regression section
- README.md: version-drift policy and * divergence documentation
- comparison report: update section 4 with confirmed root cause

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Add `refget+http://` and `refget+https://` URI schemes, opening a gtars
RefgetStore via `open_remote()` against a static HTTP(S) directory tree
(S3/R2 bucket endpoint, CDN origin, etc.).

Remote stores cache their local working set — index, alias sidecars, and
fetched sequence payloads — under a deterministic per-remote subdirectory.
The cache root resolves from `GA4GH_VRS_REFGET_CACHE_DIR`, falling back to
`$XDG_CACHE_HOME/gtars/refget` and then `~/.cache/gtars/refget`.

Document the new schemes in the `create_dataproxy` docstring, the
`vrs-annotate vcf` CLI help, the README, and the VCF annotator docs.
…izes

Resolve refget accessions without reading sequence payloads or scanning
reverse aliases:

* RefgetSequenceDataProxy and RefgetAliasDataProxy now use the store's
  get_sequence_metadata()/get_sequence_metadata_by_alias() so length and
  digest lookups no longer materialize sequence bytes.
* RefgetAliasDataProxy gains a derive_refget_accession() override that
  resolves an accession straight to a digest, and CompositeDataProxy
  routes _resolve_ga4gh_identifier through it rather than fetching and
  filtering the full alias list.
* Memoize _resolve_to_digest, and quiet the gtars store's own output.

Downgrade the alias-miss log in derive_refget_accession from exception to
debug; a missing accession is normal backend coverage, not an error, and
formatting a traceback per miss is costly.

Guard the single-accession assumption instead of relying on it:
translate_sequence_identifier filters on the "ga4gh:" namespace rather
than "SQ.", and returns a set-derived list with arbitrary order, so
silently taking [0] of a longer list would pick unpredictably. Raise
ValueError if more than one alias comes back. An empty list still means
"not found" and still returns None, which callers depend on.

Make the lru_cache sizes configurable, replacing hardcoded maxsize
literals across all three backends:

* GA4GH_VRS_DATAPROXY_CACHE_SIZE (default 8192) for lengths, digests,
  and alias lists.
* GA4GH_VRS_REST_DATAPROXY_SEQUENCE_CACHE_SIZE (default 128) for REST
  sequence payloads, kept small because an entry is bounded only by the
  caller's span.

Both accept an integer, 0 to disable, or -1/"none" for unbounded,
following the biocommons.seqrepo convention. They are read at import
time, since the decorators are applied when the classes are defined.

Add docs/dataproxy_backends.md covering the backends, URI schemes,
remote refget caching, the in-memory caches, and the stop-codon digest
divergence between refget and SeqRepo.
models is a module, so models[var["type"]] raised TypeError for every input and the "vrs" format was unusable. Resolve the class via VrsType, returning None for unknown types as intended.
SeqRepo silently truncates out-of-range fetches, so translators minted permanent VRS identifiers for locations that do not exist on their sequence, or reported them as a misleading reference mismatch against ''. Zero-width insertions past the end got no diagnostic at all, and the CNV and vrs paths never fetch sequence.

Add _DataProxy.validate_location_bounds, which raises DataProxyValidationError (unconditionally; there is no require_validation escape hatch) when any defined start/end value lies outside [0, len]. start and end are checked independently so circular start > end remains valid, pos == len is a valid insertion point, and undefined Range endpoints are skipped.

Call it from _create_allele (hgvs, spdi, beacon, gnomad), from _from_gnomad before validate_ref_seq, from CnvTranslator._from_hgvs, and from _from_vrs. The original input accession is threaded through as sequence_id and namespace-coerced with the new coerce_accession_namespace (shared with derive_refget_accession) so the length lookup is a metadata cache hit, adding no requests.
Second layer for routes the translator checks do not cover: normalize() called directly on a hand-built Allele, the annotator, and denormalize/renormalize round trips. Uses the length SequenceProxy already fetched, so it adds no I/O, and runs before the definite-range early return so those locations are checked too. Fails before bioutils can compute on truncated sequence.
Parameterized cases for validate_location_bounds (stub proxy), the normalize guard (local test SeqRepo), and every translator input path: hgvs g./n./c./p. (including p.Ter<len+1>del), spdi (including a zero-width insertion past the end), gnomad (bounds message rather than reference mismatch, with and without require_validation), beacon, vrs, and CNV hgvs. In-bounds edge cases (terminal residue, insertion at end, circular start > end, indefinite ranges) are accepted.

Each translator test uses a fresh REST dataproxy so its cassette is self-contained regardless of test order.
- Import Range at runtime instead of under TYPE_CHECKING; there is no
  import cycle. Use isinstance(pos, Range) rather than duck-typing .root.
- Fix validate_location_bounds docstring: every defined Range member is
  checked, not just the largest.
- Drop the warning log before raising; the error is always raised.
- Make the accession coercion helper private.
- Remove the passthrough _Translator._validate_location_bounds wrapper and
  call the dataproxy method directly.
- Trim tests that duplicated coverage across the helper, normalize, and
  translator levels, along with their cassettes.
Call bioutils coerce_namespace directly, which already leaves namespaced
identifiers unchanged, and restore derive_refget_accession to its
original form. Fold _defined_values into _check_location_bounds, its
only caller.
Look up the sequence length by the refget accession already in the allele
values instead of threading the input accession through every translator.
This is the same metadata lookup normalize makes next, so it adds no
request on the default path. Out-of-bounds errors from _create_allele now
name the ga4gh:SQ id. test_from_beacon (do_normalize=False) gains the
ga4gh metadata request in its cassette.
- Leave _from_vrs input as-is: like the rest of that path, it is not
  validated or normalized. Test the VrsType lookup fix directly, including
  the unknown and missing type cases.
- Validate CNV bounds against the refget accession, matching the allele
  paths.
- Call validate_location_bounds from normalize instead of importing the
  private helper, and fold the helper into the method.
- Drop redundant bounds test rows and their cassettes; add type hints.
Out-of-bounds errors now name the sequence as given (e.g. GRCh38:1)
followed by the refget accession it resolved to, taken from the aliases
in the same metadata lookup. Pass the input accession from every
translator path again so errors are consistent across formats, and so
the length lookup reuses the cached input-id metadata. This restores
test_from_beacon.yaml to its original form; bounds cassettes are
re-recorded.
The fix landed on main separately (ga4gh#656, follow-up ga4gh#668), so remove it and
its tests to keep this branch focused on location bounds validation.
RefgetStore stores written with gtars that includes databio/gtars#273 use
residue codes gtars <= 0.10 does not know, and older gtars reads them back as
different letters without raising (selenocysteine U as X). No gtars release
includes ga4gh#273 yet, so the refget extra requires >=0.11.0 and
[tool.uv.sources] pins theferrit32/gtars@7831e09f: the PR head plus the fix
that lets its Python bindings compile.

The git source is uv-only, and gtars 0.11.0 is not on PyPI, so this commit
is for local evaluation. Before merging, drop the [tool.uv.sources] entry and
keep only the version floor once a release includes ga4gh#273.
When a translator recognized its input format but the data proxy had no
refget accession for the sequence, it returned None, which translate_from
reads as "not my format". The input then fell through every other format
and surfaced as "Unable to parse data as beacon, gnomad, hgvs, spdi, vrs",
hiding that the expression was well-formed and the backend simply lacked
the sequence.

Translators now resolve accessions through require_refget_accession, which
raises UnresolvableAccessionError (a KeyError subclass) naming the
accession. The VCF annotator catches exactly that exception and leaves the
record's fields empty, as it did when the translator returned None, so a
contig the backend lacks still does not abort a file.

AliasDataProxy.derive_refget_accession now returns None, rather than
raising ValueError, for a bare accession whose namespace cannot be inferred,
matching RefgetAliasDataProxy.
bioutils.accessions.coerce_namespace has no pattern for Locus Reference
Genomic accessions, so a bare LRG_1, LRG_1t1 or LRG_1p1 could not be
resolved even by a backend holding lrg: aliases: the alias proxies treated
it as a namespace-inference failure and SeqRepo raised ValueError.

coerce_accession_namespace recognizes LRG genomic, transcript and protein
accessions as lrg: and defers to bioutils for everything else. The base,
refget and SeqRepo proxies all use it, so every backend treats a bare LRG
accession the same way.
Now that a recognized format with an unresolvable accession raises instead
of returning None, translate_from without fmt would let that exception end
the guess, so a spurious earlier match could hide a later format that
succeeds. Collect KeyError, ValueError and DataProxyValidationError per
format and continue. If no format succeeds, raise the first recognized
failure (noting any others on Python 3.11+), or the existing "Unable to
parse" ValueError when nothing recognized the input. Other exceptions are
programming errors and still abort.
validate_ref_seq fetches the reference it compares against, and an
out-of-bounds fetch may be silently truncated and misreported as a
reference mismatch. Checking bounds in the method itself protects every
caller instead of relying on each to check first, so _from_gnomad no
longer needs its own call. The bounds error is raised regardless of
require_validation, as it already was for gnomAD input.
vrs-python treats every sequence as linear. VRS allows start > end only
on a circular sequence, for a location spanning its origin, which
vrs-python does not implement. Such inputs were handled inconsistently:
pysam raised ValueError, seqrepo-rest-service answered HTTP 422, and the
CNV translator built a CopyNumberChange on the inverted location.

validate_start_le_end now rejects them with DataProxyValidationError.
validate_location_bounds applies it (for a Range, using the lowest
possible start and the highest possible end), and
_SeqRepoDataProxyBase.get_sequence applies it before reaching the
backend. Zero-width intervals (start == end) remain valid.
get_position_and_state fetches the reference for identity and
duplication edits before _create_allele validates the location. An
out-of-bounds fetch may be silently truncated, and a fetch whose start
is after its end fails with the backend's own error, so both are now
checked first and reported as DataProxyValidationError.
Some recent work done on clinvar validation surfaced a small set of
discrepancies which informed this change.

htslib does not enforce coordinates are in bounds for the sequence. This
seems like intentional and defined behavior. And pysam mostly delegates
to htslib for fetches (but disallows negative coords and `start>end`),
and seqrepo delegates to pysam, and vrs-python currently delegates to
seqrepo.

Entry points that now check bounds:
- AlleleTranslator.translate_from with fmt="hgvs", "spdi", "beacon" or
"gnomad", or with no fmt, which tries each format in turn. The check
runs before the Allele is built. For HGVS identity (`=`) and dup, it
also runs before the reference is fetched. On gnomAD it runs before the
reference comparison, so the error isn't reported as a reference
mismatch.
- CnvTranslator.translate_from (fmt="hgvs"). This path never fetches
sequence, so it previously had no protection at all.
- normalize() called directly on an Allele with a
LiteralSequenceExpression state and a SequenceReference. The check runs
before any fetch.
- vrs-annotate / the VCF annotator. It goes through the gnomAD path, so
out-of-bounds records now get VRS_Error instead of allele IDs.
- DataProxy.validate_ref_seq now checks bounds before its fetch, and
raises the bounds error regardless of `require_validation`.
- DataProxy.validate_location_bounds(sequence_id, start, end) is itself
a new public method that clients can call.

Locations whose start is after their end are also rejected now. VRS
allows `start > end` only on a circular sequence, for a location
spanning its origin, and vrs-python treats every sequence as linear.
ClinVar has a few such expressions, e.g.
`NC_000012.11:g.110593351_110576466dup`, which previously gave different
results depending on the backend.
- validate_location_bounds rejects them; a `Range` is compared by its
lowest possible start and highest possible end. This covers all the
entry points above.
- SeqRepo-backed `get_sequence` rejects them before reaching the
backend. Previously the local proxy raised pysam's `ValueError` and the
REST proxy raised an `HTTPError` (422); both now raise
`DataProxyValidationError`.
- normalize() now raises for an Allele whose `Range` start is after its
end, where it previously returned the Allele unchanged.

At a high level these are the layers and their bounds handling (prior to
this branch which changes vrs-python):

##### (begin ai)

## htslib (`faidx_fetch_seq`) (converted to interbase for ease of
comparison)

```
s = 'ABCD'
(1) fetch(s, 0, 4)      -> 'ABCD'
(2) fetch(s, -1, 3)     -> 'ABC'    # negative start clamped to 0
(3) fetch(s, 3, 1)      -> 'A'      # start clamped down to end
(4) fetch(s, 1, 100)    -> 'BCD'    # end clamped to length
(5) fetch(s, 100, 101)  -> ''
(6) fetch(s, 100, 100)  -> n/a      # inclusive-end API: can't express a zero-width interval
```

The only signal that anything was clamped is the length htslib reports
back for what it
actually fetched. The `samtools faidx` command prints a "Truncated
sequence" warning to
stderr but still exits 0.

## pysam (`FastaFile.fetch`), compared with htslib

```
(1) fetch(s, 0, 4)      -> 'ABCD'
(2) fetch(s, -1, 3)     -> ValueError: start out of range (-1)                  # new: rejects negatives
(3) fetch(s, 3, 1)      -> ValueError: invalid coordinates: start (3) > stop (1)  # new: rejects reversed
(4) fetch(s, 1, 100)    -> 'BCD'    # passed to htslib; the clamped length it reports is discarded
(5) fetch(s, 100, 101)  -> ''       # same result, but pysam returns it itself; htslib is never called
(6) fetch(s, 100, 100)  -> ''       # zero-width also returned early by pysam
```

## seqrepo (`SeqRepo.fetch`), compared with pysam

No differences: start/end are passed straight through.

```
(1) fetch(s, 0, 4)      -> 'ABCD'
(2) fetch(s, -1, 3)     -> ValueError: start out of range (-1)
(3) fetch(s, 3, 1)      -> ValueError: invalid coordinates: start (3) > stop (1)
(4) fetch(s, 1, 100)    -> 'BCD'
(5) fetch(s, 100, 101)  -> ''
(6) fetch(s, 100, 100)  -> ''
```

## vrs-python `main`, local proxy (`SeqRepoDataProxy.get_sequence`),
compared with seqrepo

No differences: `get_sequence` passes start/end straight through, and
nothing else on
`main` checks location bounds.

The one related check, `validate_ref_seq`, runs only on the gnomAD
(VCF-style) input
path. It compares the fetched reference with the input REF.

- A non-empty REF fails as a misleading "Reference mismatch … correct
ref is ''" rather
  than a bounds error.
- With `require_validation=False` (the `vrs-annotate` default), that
failure is only a
  logged warning, and an Allele is still returned.

HGVS, SPDI and beacon inputs never compare against the reference, and
neither do
`normalize()` or `translate_to`. An out-of-bounds location from those
paths is accepted
silently and given a permanent identifier.

That is why the ClinVar issue went unnoticed on SeqRepo.
`NP_001346993.1:p.Ter194del`
(a protein of 193 aa) parses to a deletion at `[193, 194)`,
normalization's fetch returns
`''`, and a `ga4gh:VA` Allele is produced. It only surfaced because
RefgetStore raised an
error on the same fetch. The comparison between the two backends then
showed different
object types: `ga4gh:VA` on SeqRepo, and on RefgetStore the harness's
CNV fallback,
`ga4gh:CX`.

```
(1) fetch(s, 0, 4)      -> 'ABCD'
(2) fetch(s, -1, 3)     -> ValueError: start out of range (-1)
(3) fetch(s, 3, 1)      -> ValueError: invalid coordinates: start (3) > stop (1)
(4) fetch(s, 1, 100)    -> 'BCD'
(5) fetch(s, 100, 101)  -> ''
(6) fetch(s, 100, 100)  -> ''
RefgetStore.open_remote defers loading the sequence index, and gtars'
alias lookups do not load it, so every accession on a refget+http(s)
store was unresolvable. Load the index when opening a remote store, via
a metadata lookup, which does load it.

The new tests build a RefgetStore from the committed test SeqRepo and
run against it both locally and over HTTP from a localhost server; the
remote accession lookups fail without this change.
@theferrit32

Copy link
Copy Markdown
Contributor Author

This is going to be re-targeted at v3 instead of main.

Brings ga4gh#663 (sequence coordinate bounds) into v3. Conflicts resolved:

- HgvsTools.get_position_and_state keeps v3's _hgvs_pos_to_vrs
  coordinates and restores ga4gh#663's bounds checks before the reference is
  fetched for identity and dup.
- CnvTranslator._from_hgvs keeps v3's build_cnv_location and checks the
  bounds of the location it builds.
- The beacon translator stays removed, along with ga4gh#663's beacon bounds
  test case and its cassette.
Brings in v3 (VRS 2.1.1 models, hgvs 2 with uncertain-range HGVS, beacon
removal) and the merge of main's ga4gh#663 into v3 (ga4gh#681), so
this branch can target v3. The beacon translator stays removed, along
with this branch's beacon case in test_from_unresolvable_accession.
@theferrit32
theferrit32 changed the base branch from main to v3 October 28, 2026 02:38

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant