Skip to content
Merged
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
1,022 changes: 510 additions & 512 deletions notebooks/Extras.ipynb

Large diffs are not rendered by default.

933 changes: 471 additions & 462 deletions notebooks/HGVS Translation.ipynb

Large diffs are not rendered by default.

81 changes: 61 additions & 20 deletions src/ga4gh/vrs/extras/translator.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
Output formats: VRS (serialized), hgvs, spdi, gnomad (vcf)

"""

from typing import Tuple
from collections.abc import Mapping
import logging
import re
Expand All @@ -25,6 +25,12 @@
_logger = logging.getLogger(__name__)


class ValidationError(Exception):
"""Class for validation errors during translation"""

pass


class Translator:
"""Translates various variation formats to and from GA4GH VRS models

Expand All @@ -37,7 +43,7 @@ class Translator:
"""

beacon_re = re.compile(r"(?P<chr>[^-]+)\s*:\s*(?P<pos>\d+)\s*(?P<ref>\w+)\s*>\s*(?P<alt>\w+)")
gnomad_re = re.compile(r"(?P<chr>[^-]+)-(?P<pos>\d+)-(?P<ref>\w+)-(?P<alt>\w+)")
gnomad_re = re.compile(r"(?P<chr>[^-]+)-(?P<pos>\d+)-(?P<ref>[ACGTN]+)-(?P<alt>[ACGTN]+|\*|\.)", re.IGNORECASE)
hgvs_re = re.compile(r"[^:]+:[cgnpr]\.")
spdi_re = re.compile(r"(?P<ac>[^:]+):(?P<pos>\d+):(?P<del_len_or_seq>\w*):(?P<ins_seq>\w*)")

Expand All @@ -56,24 +62,28 @@ def __init__(self,
self.hgvs_tools = None


def translate_from(self, var, fmt=None):
def translate_from(self, var, fmt=None, require_validation=True):
"""Translate variation `var` to VRS object

If `fmt` is None, guess the appropriate format and return the variant.
If `fmt` is specified, try only that format.
If `require_validation` is `True` then validation checks must pass in order to
return a VRS object. A `ValidationError` will be raised if validation checks
fail. If `False` then VRS object will be returned even if validation checks
fail.

See also notes about `from_` and `to_` methods.
"""

if fmt:
t = self.from_translators[fmt]
o = t(self, var)
o = t(self, var, require_validation=require_validation)
Comment thread
korikuzma marked this conversation as resolved.
if o is None:
raise ValueError(f"Unable to parse data as {fmt} variation")
return o

for fmt, t in self.from_translators.items():
o = t(self, var)
o = t(self, var, require_validation=require_validation)
if o:
return o

Expand Down Expand Up @@ -106,18 +116,18 @@ def ensure_allele_is_latest_model(self, allele):
############################################################################
## INTERNAL

def _from_beacon(self, beacon_expr, assembly_name=None):
def _from_beacon(self, beacon_expr, assembly_name=None, require_validation=True):
"""Parse beacon expression into VRS Allele

#>>> a = tlr.from_beacon("13 : 32936732 G > C")
#>>> a = tlr.from_beacon("19 : 44908822 C > T")
#>>> a.as_dict()
{'location': {'interval': {
'end': {'value': 32936732, 'type': Number},
'start': {'value': 32936731, 'type': Number},
'end': {'value': 44908822, 'type': Number},
'start': {'value': 44908821, 'type': Number},
'type': 'SequenceInterval'},
'sequence_id': 'GRCh38:13 ',
'sequence_id': 'GRCh38:19',
'type': 'SequenceLocation'},
'state': {'sequence': 'C', 'type': 'LiteralSequenceExpression'},
'state': {'sequence': 'T', 'type': 'LiteralSequenceExpression'},
'type': 'Allele'}

"""
Expand Down Expand Up @@ -148,9 +158,14 @@ def _from_beacon(self, beacon_expr, assembly_name=None):
return allele


def _from_gnomad(self, gnomad_expr, assembly_name=None):
def _from_gnomad(self, gnomad_expr, assembly_name=None, require_validation=True):
"""Parse gnomAD-style VCF expression into VRS Allele

:param str gnomad_expr: chr-pos-ref-alt
:param str assembly_name: `gnomad_expr`'s assembly
:param bool require_validation: `True` if validation checks must pass in order
for a VRS Allele to be returned. `False` otherwise.

#>>> a = tlr.from_gnomad("1-55516888-G-GA")
#>>> a.as_dict()
{'location': {'interval': {
Expand All @@ -163,7 +178,6 @@ def _from_gnomad(self, gnomad_expr, assembly_name=None):
'type': 'Allele'}

"""

if not isinstance(gnomad_expr, str):
return None
m = self.gnomad_re.match(gnomad_expr)
Expand All @@ -175,11 +189,16 @@ def _from_gnomad(self, gnomad_expr, assembly_name=None):
assembly_name = self.default_assembly_name
sequence_id = assembly_name + ":" + g["chr"]
start = int(g["pos"]) - 1
ref = g["ref"]
alt = g["alt"]
ref = g["ref"].upper()
alt = g["alt"].upper()
end = start + len(ref)
ins_seq = alt

# validation checks
valid_ref_seq, err_msg = self._is_valid_ref_seq(sequence_id, start, end, ref)
if require_validation and not valid_ref_seq:
raise ValidationError(err_msg)

interval = models.SequenceInterval(start=models.Number(value=start),
end=models.Number(value=end))
location = models.SequenceLocation(sequence_id=sequence_id, interval=interval)
Expand All @@ -189,7 +208,7 @@ def _from_gnomad(self, gnomad_expr, assembly_name=None):
return allele


def _from_hgvs(self, hgvs_expr):
def _from_hgvs(self, hgvs_expr, require_validation=True):
"""parse hgvs into a VRS object (typically an Allele)

#>>> a = tlr.from_hgvs("NM_012345.6:c.22A>T")
Expand Down Expand Up @@ -258,7 +277,7 @@ def _from_hgvs(self, hgvs_expr):
return allele


def _from_spdi(self, spdi_expr):
def _from_spdi(self, spdi_expr, require_validation=True):

"""Parse SPDI expression in to a GA4GH Allele

Expand Down Expand Up @@ -303,7 +322,7 @@ def _from_spdi(self, spdi_expr):
return allele


def _from_vrs(self, var):
def _from_vrs(self, var, require_validation=True):
"""convert from dict representation of VRS JSON to VRS object"""
if not isinstance(var, Mapping):
return None
Expand Down Expand Up @@ -451,6 +470,28 @@ def _to_spdi(self, vo, namespace="refseq"):
return spdis


def _is_valid_ref_seq(self, sequence_id: str, start_pos: int, end_pos: int,
ref: str) -> Tuple[bool, str]:
"""Return wether or not the expected reference sequence matches the actual reference sequence

:param str sequence_id: Sequence ID to use
:param int start_pos: Start pos (inter-residue) on the sequence_id
:param int end_pos: End pos (inter-residue) on the sequence_id
:param str ref: The expected reference sequence on the sequence_id given the
start_pos and end_pos
:return: Tuple containing whether or not actual reference sequence matches
the expected reference sequence and error message if mismatch
"""
actual_ref = self.data_proxy.get_sequence(sequence_id, start_pos, end_pos)
is_valid = actual_ref == ref
err_msg = ""
if not is_valid:
err_msg = f"Expected reference sequence {ref} on {sequence_id} at positions "\
f"({start_pos}, {end_pos}) but found {actual_ref}"
_logger.warning(err_msg)
return is_valid, err_msg


def is_valid_allele(self, vo):
return (vo.type == "Allele"
and vo.location.type == "SequenceLocation"
Expand Down Expand Up @@ -518,8 +559,8 @@ def _seq_id_mapper(self, ir):
expressions = [
"bogus",
"1-55516888-G-GA",
"13 : 32936732 G > C",
"NC_000013.11:g.32936732G>C",
"19 : 44908822 C > T",
"NC_000019.10:g.44908822C>T",
"NM_000551.3:21:1:T", {
"location": {
"interval": {
Expand Down
4 changes: 2 additions & 2 deletions src/ga4gh/vrs/extras/vcf_annotation.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
from python_jsonschema_objects import ValidationError

from ga4gh.vrs.dataproxy import SeqRepoDataProxy, SeqRepoRESTDataProxy
from ga4gh.vrs.extras.translator import Translator
from ga4gh.vrs.extras.translator import Translator, ValidationError as TranslatorValidationError



Expand Down Expand Up @@ -196,7 +196,7 @@ def _get_vrs_object(
"""
try:
vrs_obj = self.tlr._from_gnomad(vcf_coords, assembly_name=assembly)
except ValidationError as e:
except (ValidationError, TranslatorValidationError) as e:
_logger.error("ValidationError when translating %s from gnomad: %s", vcf_coords, str(e))
except KeyError as e:
_logger.error("KeyError when translating %s from gnomad: %s", vcf_coords, str(e))
Expand Down
Loading