diff --git a/docs/extras/vcf_annotator.md b/docs/extras/vcf_annotator.md index eced28b0..8187a526 100644 --- a/docs/extras/vcf_annotator.md +++ b/docs/extras/vcf_annotator.md @@ -35,6 +35,25 @@ Or an alternate REST path: vrs-annotate vcf --dataproxy-uri="seqrepo+http://mylabwebsite.org/seqrepo" --vcf-out=out.vcf.gz input.vcf.gz ``` +### Configuring optional RLE sequences + +Set `GA4GH_VRS_RLE_SEQ_LIMIT` before starting Python or `vrs-annotate` to change +the maximum length of optional `ReferenceLengthExpression.sequence` values. +The default is `50` bases. Use a non-negative integer, `0` to omit the optional +sequence, or `none` (case-insensitive) for no limit. Invalid values raise a +`ValueError` when VRS-Python is imported. + +The normalizer, translator, and VCF writer share this setting. Explicit +`rle_seq_limit` arguments to normalization or translation override the default; +the VCF writer still applies the configured output limit. Required +`LiteralSequenceExpression.sequence` values are unaffected. + +For example: + +```shell +GA4GH_VRS_RLE_SEQ_LIMIT=100 vrs-annotate vcf --vrs-attributes --vcf-out=out.vcf input.vcf +``` + ### Other Options `--vrs-attributes` >Will include VRS_Start, VRS_End, VRS_State fields in the INFO field. diff --git a/src/ga4gh/vrs/config.py b/src/ga4gh/vrs/config.py new file mode 100644 index 00000000..d5c865fc --- /dev/null +++ b/src/ga4gh/vrs/config.py @@ -0,0 +1,30 @@ +"""Process-level configuration for VRS-Python.""" + +import os + + +def _get_rle_seq_limit() -> int | None: + """Read the maximum optional RLE sequence length from the environment. + + Set ``GA4GH_VRS_RLE_SEQ_LIMIT`` to change the maximum length of optional + ``ReferenceLengthExpression.sequence`` values. The default is 50 bases. + Use a non-negative integer, 0 to omit the optional sequence, or ``none`` + (case-insensitive) for no limit. + + :raise ValueError: If the value is not an integer or is negative, unless + it is ``none`` (case-insensitive). + """ + value = os.environ.get("GA4GH_VRS_RLE_SEQ_LIMIT", "50") + if value.lower() == "none": + return None + message = "GA4GH_VRS_RLE_SEQ_LIMIT must be a non-negative integer or 'none'" + try: + limit = int(value) + except ValueError as exc: + raise ValueError(message) from exc + if limit < 0: + raise ValueError(message) + return limit + + +RLE_SEQ_LIMIT = _get_rle_seq_limit() diff --git a/src/ga4gh/vrs/extras/annotator/vcf.py b/src/ga4gh/vrs/extras/annotator/vcf.py index 3dd70dcf..f118d758 100644 --- a/src/ga4gh/vrs/extras/annotator/vcf.py +++ b/src/ga4gh/vrs/extras/annotator/vcf.py @@ -14,6 +14,7 @@ use_ga4gh_compute_identifier_when, ) from ga4gh.vrs import VRS_VERSION, VrsType, __version__ +from ga4gh.vrs.config import RLE_SEQ_LIMIT from ga4gh.vrs.dataproxy import _DataProxy from ga4gh.vrs.extras.translator import AlleleTranslator from ga4gh.vrs.models import Allele, Range @@ -64,11 +65,6 @@ class FieldName(str, Enum): } ) -# ReferenceLengthExpression .sequence values will be included in output VCF if -# length <= this value. This field is optional for RLE since it can be derived -# from the reference sequence. Set to None to always include the sequence. -RLE_SEQ_LIMIT = 50 - def dump_alleles_to_pkl(alleles: list[Allele], output_pkl_path: Path) -> None: """Create pkl file of dictionary mapping VRS IDs to ingested alleles. diff --git a/src/ga4gh/vrs/extras/translator.py b/src/ga4gh/vrs/extras/translator.py index 91c45532..0adcdba5 100644 --- a/src/ga4gh/vrs/extras/translator.py +++ b/src/ga4gh/vrs/extras/translator.py @@ -15,6 +15,7 @@ from ga4gh.core import ga4gh_identify from ga4gh.vrs import models, normalize +from ga4gh.vrs.config import RLE_SEQ_LIMIT from ga4gh.vrs.dataproxy import SequenceProxy, _DataProxy from ga4gh.vrs.extras.decorators import lazy_property from ga4gh.vrs.normalize import denormalize_reference_length_expression @@ -74,7 +75,7 @@ def __init__( data_proxy: _DataProxy, default_assembly_name: str = "GRCh38", identify: bool = True, - rle_seq_limit: int | None = 50, + rle_seq_limit: int | None = RLE_SEQ_LIMIT, ) -> None: self.default_assembly_name = default_assembly_name self.data_proxy = data_proxy diff --git a/src/ga4gh/vrs/normalize.py b/src/ga4gh/vrs/normalize.py index fad73d87..8e3b2511 100644 --- a/src/ga4gh/vrs/normalize.py +++ b/src/ga4gh/vrs/normalize.py @@ -15,6 +15,7 @@ from ga4gh.core import ga4gh_digest, is_pydantic_instance, pydantic_copy from ga4gh.vrs import models +from ga4gh.vrs.config import RLE_SEQ_LIMIT from ga4gh.vrs.dataproxy import SequenceProxy, _DataProxy _logger = logging.getLogger(__name__) @@ -85,7 +86,9 @@ def _get_new_allele_location_pos( def _normalize_allele( - input_allele: models.Allele, data_proxy: _DataProxy, rle_seq_limit: int = 50 + input_allele: models.Allele, + data_proxy: _DataProxy, + rle_seq_limit: int | None = RLE_SEQ_LIMIT, ): """Normalize Allele using "fully-justified" normalization adapted from NCBI's VOCA. Fully-justified normalization expands such ambiguous representation over the diff --git a/tests/extras/test_rle_config.py b/tests/extras/test_rle_config.py new file mode 100644 index 00000000..03b87702 --- /dev/null +++ b/tests/extras/test_rle_config.py @@ -0,0 +1,91 @@ +from importlib import import_module, reload +from pathlib import Path + +import pysam +import pytest +import vcr + +from ga4gh.vrs import config +from ga4gh.vrs.dataproxy import SeqRepoRESTDataProxy + + +@pytest.fixture +def configured_modules(monkeypatch, request): + """Reload import-time defaults, restoring them after each environment case.""" + modules = [ + config, + import_module("ga4gh.vrs.normalize"), + import_module("ga4gh.vrs.extras.translator"), + import_module("ga4gh.vrs.extras.annotator.vcf"), + ] + original_namespaces = [vars(module).copy() for module in modules] + try: + with monkeypatch.context() as env: + env.delenv("GA4GH_VRS_RLE_SEQ_LIMIT", raising=False) + if request.param is not None: + env.setenv("GA4GH_VRS_RLE_SEQ_LIMIT", request.param) + for module in modules: + reload(module) + yield modules[1:] + finally: + # Preserve class identities already imported by other test modules. + for module, namespace in zip(modules, original_namespaces, strict=True): + vars(module).clear() + vars(module).update(namespace) + + +@pytest.mark.parametrize( + ("configured_modules", "expected_sequence"), + [ + (None, "CTTTCTTT"), + ("0", None), + ("4", None), + ("8", "CTTTCTTT"), + ("none", "CTTTCTTT"), + ("NoNe", "CTTTCTTT"), + ], + indirect=["configured_modules"], +) +def test_rle_sequence_limit_environment( + configured_modules, expected_sequence, tmp_path +): + """Normalization, translation and VCF output share the configured default.""" + normalization, translator, vcf_module = configured_modules + output = tmp_path / "annotated.vcf" + with vcr.use_cassette( + "tests/extras/cassettes/test_annotate_vcf_rle.yaml", + record_mode="none", + allow_playback_repeats=True, + ): + proxy = SeqRepoRESTDataProxy( + base_url="http://localhost:5000/seqrepo", + disable_healthcheck=True, + ) + tlr = translator.AlleleTranslator(proxy) + variant = "1-102995989-CTTT-CTTTCTTT" + state = tlr.translate_from(variant, fmt="gnomad").state + assert ( + state.sequence.root if state.sequence is not None else None + ) == expected_sequence + literal = tlr.translate_from(variant, fmt="gnomad", do_normalize=False) + normalized = normalization.normalize(literal, proxy).state + assert ( + normalized.sequence.root if normalized.sequence is not None else None + ) == expected_sequence + explicit = tlr.translate_from(variant, fmt="gnomad", rle_seq_limit=None).state + assert explicit.sequence.root == "CTTTCTTT" + vcf_module.VcfAnnotator(proxy).annotate( + Path("tests/extras/data/test_rle.vcf"), + output, + vrs_attributes=True, + ) + with pysam.VariantFile(output) as vcf: + records = list(vcf) + assert records[1].info["VRS_States"][-1] == (expected_sequence or ".") + + +@pytest.mark.parametrize("limit", ["-1", "invalid", "1.5"]) +def test_invalid_rle_sequence_limit_environment(limit, monkeypatch): + monkeypatch.setenv("GA4GH_VRS_RLE_SEQ_LIMIT", limit) + with pytest.raises(ValueError, match="GA4GH_VRS_RLE_SEQ_LIMIT"): + config._get_rle_seq_limit()