Repository: Mykrobe-tools/mykrobe

Issue:
X mutations in 201901 panel
Dear mykrobe team,

In the 201901 probe set, many (most?) putative mutations are annotated as "X". For example, all resistance mutations at position 445 in the rpoB gene is now annotated as H445X.

Our clinicians are very interested in the actual amino acid change, so at the moment I am doing this with a fairly flimsy in-house script. What I'm wondering is if there are any plans to change the probe annotations back to their specific amino acid changes, or if there are any other tricks to getting the actual changes annotated in the output?

Best regards,
Ola Brynildsrud
Norwegian Institute of Public Health

Discussion:

Patch:
diff --git a/src/mykrobe/cmds/amr.py b/src/mykrobe/cmds/amr.py
index 8901564..fce6c0b 100755
--- a/src/mykrobe/cmds/amr.py
+++ b/src/mykrobe/cmds/amr.py
@@ -25,6 +25,7 @@ from mykrobe.predict import MykrobePredictorSusceptibilityResult
 from mykrobe.metagenomics import AMRSpeciesPredictor
 from mykrobe.species_data import DataDir
 from mykrobe.utils import load_json
+from mykrobe.utils import fix_amino_acid_X_variants_keys
 from mykrobe.version import __version__ as predictor_version
 from mykrobe.version import __version__ as atlas_version
 
@@ -145,7 +146,9 @@ class ConfThresholder:
 def ref_data_from_args(args):
     if args.species == "custom":
         if args.custom_probe_set_path is None:
-            raise ValueError("Must use --custom_probe_set_path option if the species is 'custom'")
+            raise ValueError(
+                "Must use --custom_probe_set_path option if the species is 'custom'"
+            )
         ref_data = {
             "fasta_files": [args.custom_probe_set_path],
             "var_to_res_json": args.custom_variant_to_resistance_json,
@@ -184,7 +187,9 @@ def detect_species_and_get_depths(cov_parser, hierarchy_json, wanted_phylo_group
         return {}, depths
 
     species_predictor = AMRSpeciesPredictor(
-        phylo_group_covgs=cov_parser.covgs.get("complex", cov_parser.covgs.get("phylo_group", {})),
+        phylo_group_covgs=cov_parser.covgs.get(
+            "complex", cov_parser.covgs.get("phylo_group", {})
+        ),
         sub_complex_covgs=cov_parser.covgs.get("sub-complex", {}),
         species_covgs=cov_parser.covgs["species"],
         lineage_covgs=cov_parser.covgs.get("sub-species", {}),
@@ -193,7 +198,11 @@ def detect_species_and_get_depths(cov_parser, hierarchy_json, wanted_phylo_group
     phylogenetics = species_predictor.run()
 
     if wanted_phylo_group in species_predictor.out_json["phylogenetics"]["phylo_group"]:
-        depths = [species_predictor.out_json["phylogenetics"]["phylo_group"][wanted_phylo_group]["median_depth"]]
+        depths = [
+            species_predictor.out_json["phylogenetics"]["phylo_group"][
+                wanted_phylo_group
+            ]["median_depth"]
+        ]
     return phylogenetics, depths
 
 
@@ -231,15 +240,33 @@ def write_outputs(args, base_json):
             print(output)
 
 
+def fix_X_amino_acid_variants(sample_json):
+    if "susceptibility" in sample_json:
+        for drug_dict in sample_json["susceptibility"].values():
+            if "called_by" in drug_dict:
+                fix_amino_acid_X_variants_keys(drug_dict["called_by"])
+
+    if "variant_calls" in sample_json:
+        fix_amino_acid_X_variants_keys(sample_json["variant_calls"])
+
+
 def run(parser, args):
     logger.info(f"Start runnning mykrobe predict. Command line: {' '.join(sys.argv)}")
     base_json = {args.sample: {}}
     args = parser.parse_args()
     ref_data = ref_data_from_args(args)
-    if args.species == "custom" and ref_data["var_to_res_json"] is None and ref_data["lineage_json"] is None:
-        logger.info("Forcing --report_all_calls because species is 'custom' and options --custom_variant_to_resistance_json,--custom_lineage_json were not used")
+    if (
+        args.species == "custom"
+        and ref_data["var_to_res_json"] is None
+        and ref_data["lineage_json"] is None
+    ):
+        logger.info(
+            "Forcing --report_all_calls because species is 'custom' and options --custom_variant_to_resistance_json,--custom_lineage_json were not used"
+        )
         args.report_all_calls = True
-    logger.info(f"Running mykrobe predict using species {args.species}, and panel version {ref_data['version']}")
+    logger.info(
+        f"Running mykrobe predict using species {args.species}, and panel version {ref_data['version']}"
+    )
 
     # Run Cortex
     cp = CoverageParser(
@@ -260,7 +287,9 @@ def run(parser, args):
         phylogenetics = {}
         depths = [cp.estimate_depth()]
     else:
-        phylogenetics, depths = detect_species_and_get_depths(cp, ref_data["hierarchy_json"], ref_data["species_phylo_group"])
+        phylogenetics, depths = detect_species_and_get_depths(
+            cp, ref_data["hierarchy_json"], ref_data["species_phylo_group"]
+        )
 
     # Genotype
     variant_calls_dict = {}
@@ -326,7 +355,10 @@ def run(parser, args):
         if args.conf_percent_cutoff < 100:
             logger.debug("Expected depth: " + str(depths[0]))
             conf_thresholder = ConfThresholder(
-                kmer_count_error_rate, depths[0], ref_data["kmer"], incorrect_kmer_to_pc_cov
+                kmer_count_error_rate,
+                depths[0],
+                ref_data["kmer"],
+                incorrect_kmer_to_pc_cov,
             )
             time_start = time.time()
             conf_threshold = conf_thresholder.get_conf_threshold(
@@ -369,7 +401,11 @@ def run(parser, args):
         depths = [cp.estimate_depth()]
 
     mykrobe_predictor_susceptibility_result = MykrobePredictorSusceptibilityResult()
-    if gt is not None and (max(depths) > args.min_depth or args.force) and ref_data["var_to_res_json"] is not None:
+    if (
+        gt is not None
+        and (max(depths) > args.min_depth or args.force)
+        and ref_data["var_to_res_json"] is not None
+    ):
         predictor = BasePredictor(
             variant_calls=gt.variant_calls,
             called_genes=gt.sequence_calls_dict,
@@ -383,15 +419,22 @@ def run(parser, args):
         logger.info("Progress: finished making AMR predictions")
 
     base_json[args.sample] = {
-        "susceptibility": list(mykrobe_predictor_susceptibility_result.to_dict().values())[0],
-        "phylogenetics": {} if phylogenetics == {} else list(phylogenetics.to_dict().values())[0],
+        "susceptibility": list(
+            mykrobe_predictor_susceptibility_result.to_dict().values()
+        )[0],
+        "phylogenetics": {}
+        if phylogenetics == {}
+        else list(phylogenetics.to_dict().values())[0],
         "variant_calls": variant_calls_dict,
         "sequence_calls": sequence_calls_dict,
         "lineage_calls": lineage_calls_dict,
         "kmer": ref_data["kmer"],
         "probe_sets": ref_data["fasta_files"],
         "files": args.seq,
-        "version": {"mykrobe-predictor": predictor_version, "mykrobe-atlas": atlas_version},
+        "version": {
+            "mykrobe-predictor": predictor_version,
+            "mykrobe-atlas": atlas_version,
+        },
         "genotype_model": args.model,
     }
     if len(lineage_predict_dict) > 0:
@@ -401,5 +444,6 @@ def run(parser, args):
         cp.remove_temporary_files()
 
     logger.info("Progress: writing output")
+    fix_X_amino_acid_variants(base_json[args.sample])
     write_outputs(args, base_json)
     logger.info("Progress: finished")
diff --git a/src/mykrobe/utils.py b/src/mykrobe/utils.py
index 2a7f58c..889d8c8 100644
--- a/src/mykrobe/utils.py
+++ b/src/mykrobe/utils.py
@@ -3,13 +3,16 @@ import hashlib
 import re
 import json
 
+from Bio.Seq import Seq
+
 
 def check_args(args):
     if args.db_name is None:
         args.db_name = os.environ.get("DB_NAME")
     if args.db_name is None:
         raise ValueError(
-            "db_name needs to be set. Either run with --db_name :db_name or export DB_NAME=:db_name")
+            "db_name needs to be set. Either run with --db_name :db_name or export DB_NAME=:db_name"
+        )
     return args
 
 
@@ -18,8 +21,7 @@ def make_hash(s):
 
 
 def make_var_hash(ref, pos, alts):
-    var = "".join(
-        [ref, str(pos), "/".join(alts)])
+    var = "".join([ref, str(pos), "/".join(alts)])
     return make_hash(var)
 
 
@@ -44,7 +46,7 @@ def get_params(url):
     except IndexError:
         return params
     p_str = p_str.split(" ")[0]
-    p_str = p_str.split('&')
+    p_str = p_str.split("&")
     for p in p_str:
         k, v = p.split("=")
         params[k] = v
@@ -58,27 +60,73 @@ def median(lst):
     lstLen = len(lst)
     index = (lstLen - 1) // 2
 
-    if (lstLen % 2):
+    if lstLen % 2:
         return sortedLst[index]
     else:
         return (sortedLst[index] + sortedLst[index + 1]) / 2.0
 
 
 def load_json(f):
-    with open(f, 'r') as infile:
+    with open(f, "r") as infile:
         return json.load(infile)
 
 
 def lazyprop(fn):
-    attr_name = '_' + fn.__name__
+    attr_name = "_" + fn.__name__
 
     @property
     def _lazyprop(self):
         if not hasattr(self, attr_name):
             setattr(self, attr_name, fn(self))
         return getattr(self, attr_name)
+
     return _lazyprop
 
+
 def seq_to_kmers(seq, kmer_size):
-    for i in range(len(seq)-kmer_size+1):
-        yield seq[i:i+kmer_size]    
+    for i in range(len(seq) - kmer_size + 1):
+        yield seq[i : i + kmer_size]
+
+
+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-)
+    # Example var_name with an X: "katG_S315X-GCT2155167GGT"
+    #                                       |           |||
+    #                                       aa          codon
+    match = re.match(
+        r"""(?P<prefix>.*_[A-Z][0-9]+)(?P<aa>[A-Z])(?P<middle>-[ACGT]{3}[0-9]+)(?P<codon>[ACGT]{3})""",
+        var_name,
+    )
+    if match is None or match.group("aa") != "X":
+        return None
+    try:
+        amino_acid = str(Seq(match.group("codon")).translate())
+    except:
+        return None
+    return (
+        match.group("prefix")
+        + amino_acid
+        + match.group("middle")
+        + match.group("codon")
+    )
+
+
+def fix_amino_acid_X_variants_keys(dict_to_fix):
+    """The way panels and variants work mean that the 'any amino acid' variants
+    look like eg H445X, where X is any amino acid other than H. Users want to
+    know the actual change. This function changes all keys in dict_to_fix that
+    have those variants, changing the X to the actual amino acid."""
+    keys_to_replace = {}
+    for key in dict_to_fix:
+        new_key = _x_mutation_fixed_var_name(key)
+        if new_key is not None:
+            assert new_key not in keys_to_replace
+            assert new_key not in dict_to_fix
+            keys_to_replace[key] = new_key
+
+    for key, new_key in keys_to_replace.items():
+        dict_to_fix[new_key] = dict_to_fix[key]
+        del dict_to_fix[key]


Test patch:
diff --git a/tests/utils_test.py b/tests/utils_test.py
new file mode 100644
index 0000000..1a9ac1e
--- /dev/null
+++ b/tests/utils_test.py
@@ -0,0 +1,26 @@
+import pytest
+
+from mykrobe import utils
+
+
+def test_x_mutation_fixed_var_name():
+    assert "katG_S315G-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
+
+
+def test_fix_amino_acid_X_variants_keys():
+    test_dict = {
+        "foo": "bar",
+        "katG_S315X-GCT2155167GGT": "baz",
+        "katG_S315C-GCT2155167CTT": "baz",
+    }
+
+    utils.fix_amino_acid_X_variants_keys(test_dict)
+    assert test_dict == {
+        "foo": "bar",
+        "katG_S315G-GCT2155167GGT": "baz",
+        "katG_S315C-GCT2155167CTT": "baz",
+    }