Skip to content

Enforce sequence coordinate bounds - #663

Open
theferrit32 wants to merge 10 commits into
mainfrom
kf/location-bounds
Open

theferrit32 wants to merge 10 commits into
mainfrom
kf/location-bounds

Conversation

@theferrit32

@theferrit32 theferrit32 commented Oct 1, 2026 •

Copy link
Copy Markdown
Contributor

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. 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_location_bounds(sequence_id, start, end) is itself a new public method that clients can call.

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

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.
@theferrit32 theferrit32 self-assigned this Oct 1, 2026
@theferrit32
theferrit32 marked this pull request as ready for review October 1, 2026 05:11
@theferrit32
theferrit32 requested review from a team as code owners October 1, 2026 05:11

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