Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 19 additions & 0 deletions docs/extras/vcf_annotator.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
30 changes: 30 additions & 0 deletions src/ga4gh/vrs/config.py
Original file line number Diff line number Diff line change
@@ -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()
6 changes: 1 addition & 5 deletions src/ga4gh/vrs/extras/annotator/vcf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
3 changes: 2 additions & 1 deletion src/ga4gh/vrs/extras/translator.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion src/ga4gh/vrs/normalize.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__)
Expand Down Expand Up @@ -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
Expand Down
91 changes: 91 additions & 0 deletions tests/extras/test_rle_config.py

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm wondering if we could avoid using subprocess with a full script (it's a little hard to read and I imagine it could be hard to maintain). Could we instead look into mocking the environment variables in our existing tests (such as

def test_rle_seq_limit(tlr):
and https://github.com/ga4gh/vrs-python/blob/main/tests/extras/test_annotate_vcf.py)?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yeah generally the recommended way to injectj/manage env vars in tests is with the pytest monkeypatch module

https://docs.pytest.org/en/stable/how-to/monkeypatch.html#monkeypatching-environment-variables

Original file line number Diff line number Diff line change
@@ -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()