#591: Use mapper boundary rather than seq length to determine when to…

… infer ref from genomic sequence

hgvs/sequencevariant.py

@@ -59,14 +59,14 @@ def __repr__(self):
     def fill_ref(self, hdp):
         # TODO: Refactor. SVs should not operate on themselves when
         # external resources are required
-        hm = hgvs.variantmapper.VariantMapper(hdp)
+        vm = hgvs.variantmapper.VariantMapper(hdp)
         type = None
         if isinstance(self.posedit, hgvs.posedit.PosEdit) and isinstance(
                 self.posedit.edit, hgvs.edit.Edit):
             type = self.posedit.edit.type
         if type in ["del", "delins", "identity", "dup", "repeat"
                     ] and self.posedit.edit.ref_s is None:
-            hm._replace_reference(self)
+            vm._replace_reference(self)
         if type == "identity" and isinstance(self.posedit.edit, hgvs.edit.NARefAlt):
             self.posedit.edit.alt = self.posedit.edit.ref
         return self

hgvs/variantmapper.py

@@ -175,7 +175,9 @@ def g_to_n(self, var_g, tx_ac, alt_aln_method=hgvs.global_config.mapping.alt_aln
         pos_n.uncertain = var_g.posedit.pos.uncertain
         var_n = hgvs.sequencevariant.SequenceVariant(
             ac=tx_ac, type="n", posedit=hgvs.posedit.PosEdit(pos_n, edit_n))
-        if self.replace_reference:
+        if (self.replace_reference
+            and var_n.posedit.pos.start.base >= 0
+            and var_n.posedit.pos.end.base < mapper.tgt_len):
             self._replace_reference(var_n)
         if self.add_gene_symbol:
             self._update_gene_symbol(var_n, var_g.gene)

tests/issues/test_437.py

@@ -118,24 +118,6 @@ def test_enforce_strict_bounds(self):
             e2_g = self.am37.n_to_g(self.e2_n)
 
 
-
-def test_invitae_examples(parser, am37):
-    tests = [
-        # ("NC_000009.11:g.35657741A>G", "NR_003051.3:n.*7T>C"), # Grammar doesn't support n.*
-        ("NC_000009.11:g.35658020C>T", "NR_003051.3:n.-5G>A"),
-        ]
-
-    hgvs.global_config.mapping.strict_bounds = False
-
-    for hgvs_g, hgvs_t in tests:
-        var_g = parser.parse(hgvs_g)
-        var_t = parser.parse(hgvs_t)
-        assert hgvs_t == str(am37.g_to_t(var_g, var_t.ac))
-
-    hgvs.global_config.mapping.strict_bounds = True
-
-
-
 def test_oob_dup(parser, am37):
     """Intentionally preserve dup, derived from genomic sequence, when
     projecting to out-of-bounds transcript coordinates
@@ -156,6 +138,28 @@ def test_oob_dup(parser, am37):
 
 
 
+def test_invitae_examples(parser, am37):
+    tests = [
+        ("NC_000009.11:g.35658020C>T", "NR_003051.3:n.-5G>A"),
+
+        # ("NC_000009.11:g.35657741A>G", "NR_003051.3:n.*7T>C"),
+        # See #592: hgvs doesn't support n.* coordinates, so
+        # rewrite of previous variant with SEQ_START coordinate
+        ("NC_000009.11:g.35657741A>G", "NR_003051.3:n.275T>C"),
+        ]
+
+    hgvs.global_config.mapping.strict_bounds = False
+
+    for hgvs_g, hgvs_t in tests:
+        var_g = parser.parse(hgvs_g)
+        var_t = parser.parse(hgvs_t)
+        assert hgvs_t == str(am37.g_to_t(var_g, var_t.ac))
+
+    hgvs.global_config.mapping.strict_bounds = True
+
+
+
+
 if __name__ == "__main__":
     from hgvs.easy import *
     test_437(parser=parser, am=am37)