Repository: Mykrobe-tools/mykrobe

Issue:
mykrobe calls the wrong amino acid change for first drug resistance SNP on a bottom strand gene
I ran Mykrobe on the FASTQ data for BioProject PRJEB5280 Run ERR751367 from NCBI. Mykrobe identifies a SNP conferring isoniazid resistance, katG_S315G-GCT2155167GGT. The resulting amino acid, in this case the G in S315G, is incorrect. The last part of the SNP, GCT2155167GGT, identifies the codon that was changed and the position within the genome. These nucleotides are always the top strand nucleotides. Therefore, for genes on the bottom strand, we must take the reverse complement. GCT becomes AGC, which is serine, the S in S315G. So far so good. However, the reverse complement of GGT is ACC, which is threonine, but the output S315G says the resulting amino acid is glycine. I believe the code is neglecting to use the reverse complement, so it is incorrectly translating GGT to glycine.

I have also tried running Mykrobe on .bam files that I generated by running snippy (version 3.2) on the same FASTQ data downloaded from NCBI. I see the same error, and looking manually at the bam file in IGV confirms that S315G is incorrect and should be S315T.

I suspect that this bug only affects the first drug resistant SNP listed in the json output that is in a bottom strand gene. katG_S315G-GCT2155167GGT just happens to be an isoniazid SNP, so it is near the top of the file and thus often the SNP affected.

Identical cases for katG_S315G-GCT2155167GGT occur with the following NCBI BioProject and Run pairs:
PRJEB5280, ERR779909
PRJEB10385, ERR1034759
PRJEB11653, ERR1213902

I have ~200 more identical cases. I can provide them if they would be helpful for debugging.

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 86b3c90..a2cc13b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,10 +32,15 @@ this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 
 - Bug where the `--ont` flag resulted in mykrobe predict crashing if the panel
   contained any presence/absence genes [[#123][123]].
+- Bug reporting variants that are "amino acid to any other amino acid" changes,
+  eg `katG_S315X`. If gene is on reverse strand, then the wrong new amino
+  acid change was being reported [[#127][127]].
+
 
 [113]: https://github.com/Mykrobe-tools/mykrobe/issues/113
 [114]: https://github.com/Mykrobe-tools/mykrobe/issues/114
 [123]: https://github.com/Mykrobe-tools/mykrobe/issues/123
+[127]: https://github.com/Mykrobe-tools/mykrobe/issues/127
 [Unreleased]: https://github.com/Mykrobe-tools/mykrobe/compare/v0.9.0...HEAD
 [mkdtemp]: https://docs.python.org/3.6/library/tempfile.html#tempfile.mkdtemp
 
diff --git a/src/mykrobe/utils.py b/src/mykrobe/utils.py
index 889d8c8..047b071 100644
--- a/src/mykrobe/utils.py
+++ b/src/mykrobe/utils.py
@@ -92,25 +92,53 @@ def _x_mutation_fixed_var_name(var_name):
     """Takes mutation name from results base_json. If it is an "X" amino acid
     mutation, returns new name, where X is fixed with correct amino acid name.
     Otherwise returns None"""
-    #                              (prefix-) (--middle-)
+    #                             prefix pos1     pos2
+    #                              ||||| |||     |||||||
     # Example var_name with an X: "katG_S315X-GCT2155167GGT"
-    #                                       |           |||
-    #                                       aa          codon
+    #                                   |   | |||       |||
+    #                                 aa1 aa2 codon1   codon2
     match = re.match(
-        r"""(?P<prefix>.*_[A-Z][0-9]+)(?P<aa>[A-Z])(?P<middle>-[ACGT]{3}[0-9]+)(?P<codon>[ACGT]{3})""",
+        r"""(?P<prefix>.*_)(?P<aa1>[A-Z])(?P<pos1>[0-9]+)(?P<aa2>[A-Z])-(?P<codon1>[ACGT]{3})(?P<pos2>[0-9]+)(?P<codon2>[ACGT]{3})""",
         var_name,
     )
-    if match is None or match.group("aa") != "X":
+    if match is None or match.group("aa2") != "X":
         return None
+
     try:
-        amino_acid = str(Seq(match.group("codon")).translate())
+        codon1_translate = str(Seq(match.group("codon1")).translate())
+        codon1_rev_translate = str(
+            Seq(match.group("codon1")).reverse_complement().translate()
+        )
     except:
         return None
+
+    # Check which strand the gene is on by checking if first codon or its
+    # reverse complement match the amino acid. This is a safe way of
+    # determining the strand because there is no codon where its translation
+    # and the translation of its reverse complement results in the same
+    # amino acid (this is checked in the tests of this repo)
+    if codon1_translate == match.group("aa1"):
+        try:
+            new_aa = str(Seq(match.group("codon2")).translate())
+        except:
+            return None
+    elif codon1_rev_translate == match.group("aa1"):
+        try:
+            new_aa = str(Seq(match.group("codon2")).reverse_complement().translate())
+        except:
+            return None
+    else:
+        return None
+
     return (
         match.group("prefix")
-        + amino_acid
-        + match.group("middle")
-        + match.group("codon")
+        + match.group("aa1")
+        + match.group("pos1")
+        + new_aa
+        + "-"
+        + match.group("codon1")
+        + match.group("pos2")
+        + match.group("codon2")
     )
 
 


Test patch:
diff --git a/tests/utils_test.py b/tests/utils_test.py
index 1a9ac1e..fbc730c 100644
--- a/tests/utils_test.py
+++ b/tests/utils_test.py
@@ -1,26 +1,58 @@
 import pytest
 
+from Bio.Seq import Seq
+
 from mykrobe import utils
 
+# When we have an amino acid to any other amino mutation (eg I42X), we want to
+# replace the X with the actual amino acid. This means getting the codon from
+# the variant name and translating it. We have the variants in the form eg:
+#   katG_S315X-GCT2155167GGT
+# Problem is, if the gene is on the reverse strand then we have to reverse
+# complement the codon. Now, we don't actually know the strand that the gene is
+# on, but we do know the first codon and the amino acid - in that example they
+# are GCT and S. This test is not testing the mykrobe code, but is testing
+# that we can use the first codon to check if a gene is on the reverse strand.
+# For this to work, we can't have a codon and its reverse complement translating
+# to the same amino acid.
+# Eg for GCT -> S, the reverse complement is AGC, which translates to A. This
+# test runs through all codons and checks that we never get the same amino
+# acid from a codon and the codon's reverse complement.
+def test_x_mutation_no_amino_acid_revcomp_clash():
+    nucs = ["A", "C", "G", "T"]
+    codons = [f"{x}{y}{z}" for x in nucs for y in nucs for z in nucs]
+
+    for codon in codons:
+        amino_acid = str(Seq(codon).translate())
+        amino_acid_rev = str(Seq(codon).reverse_complement().translate())
+        assert amino_acid != amino_acid_rev
+
 
 def test_x_mutation_fixed_var_name():
-    assert "katG_S315G-GCT2155167GGT" == utils._x_mutation_fixed_var_name(
+    # This is on the reverse strand, so the GGT should get reverse complemented,
+    # resulting in amino acid T (=ACC), not G (=GGT).
+    assert "katG_S315T-GCT2155167GGT" == utils._x_mutation_fixed_var_name(
         "katG_S315X-GCT2155167GGT"
     )
     assert utils._x_mutation_fixed_var_name("katG_S315.-GCT2155167GGT") is None
     assert utils._x_mutation_fixed_var_name("katG_S315X-GCT2155167NNN") is None
 
+    # This is on forward strand, so should get TAT -> Y
+    assert "embB_D328Y-GAT4247495TAT" == utils._x_mutation_fixed_var_name(
+        "embB_D328X-GAT4247495TAT"
+    )
+
 
 def test_fix_amino_acid_X_variants_keys():
     test_dict = {
         "foo": "bar",
         "katG_S315X-GCT2155167GGT": "baz",
-        "katG_S315C-GCT2155167CTT": "baz",
+        "katG_S315C-GCT2155167TTA": "baz",
     }
 
     utils.fix_amino_acid_X_variants_keys(test_dict)
     assert test_dict == {
         "foo": "bar",
-        "katG_S315G-GCT2155167GGT": "baz",
-        "katG_S315C-GCT2155167CTT": "baz",
+        "katG_S315T-GCT2155167GGT": "baz",
+        "katG_S315C-GCT2155167TTA": "baz",
     }