Update genotype API for latest version of mykatlas

mykrobe/cmds/amr.py

@@ -186,7 +186,7 @@ def run(parser, args):
         depths = cp.estimate_depth()
     args.quiet = q
     mykrobe_predictor_susceptibility_result = MykrobePredictorSusceptibilityResult()
-    if Predictor is not None:
+    if Predictor is not None and max(depths) > args.min_depth:
         predictor = Predictor(variant_calls=gt.variant_calls,
                               called_genes=gt.gene_presence_covgs,
                               base_json=base_json[args.sample],

mykrobe/mykrobe_predictor.py

@@ -78,26 +78,21 @@ def main():
         metavar='min_depth',
         type=int,
         help='min_depth',
-        default=3)
+        default=1)
     parser_amr.set_defaults(func=run_subtool)
 
     # ##########
     # # Genotype
     # ##########
-    # ##########
-    # # Genotype
-    # ##########
+
     parser_geno = subparsers.add_parser(
         'genotype',
         parents=[
             sequence_or_binary_parser_mixin,
             probe_set_mixin,
-            force_mixin],
+            force_mixin,
+            genotyping_mixin],
         help='genotype a sample using a probe set')
-    parser_geno.add_argument(
-        '--ignore_filtered',
-        help="don't include filtered genotypes",
-        default=False)
     parser_geno.set_defaults(func=run_subtool)
 
     args = parser.parse_args()

mykrobe/predict/amr.py

@@ -10,6 +10,7 @@
 from ga4ghmongo.schema import SequenceCall
 
 from pprint import pprint
+import logging
 
 
 DEFAULT_MIN_GENE_CN = 0.03
@@ -27,10 +28,10 @@ def copy_number(call):
 
     return round(float(alternate_depth) / (alternate_depth + wt_depth), 2)
 
-def depth_on_allele(call):
+def depth_on_alternate(call):
     coverage = call.info.get("coverage")
     try:
-        alternate_depth = coverage.get("alternate").get("median_depth") + coverage.get("reference").get("median_depth")
+        alternate_depth = coverage.get("alternate").get("median_depth")
     except:
         alternate_depth = coverage.get("median_depth")
     return alternate_depth
@@ -62,7 +63,7 @@ def __init__(self, variant_calls, called_genes, base_json={}, depth_threshold =
 
     def _create_initial_resistance_prediction(self):
         self.result =  MykrobePredictorSusceptibilityResult(dict(
-            (k, {"predict": "I"}) for k in self.drugs))
+            (k, {"predict": "N"}) for k in self.drugs))
         self.resistance_predictions = self.result.susceptibility
 
     def _get_drug_list_from_variant_to_resistance_drug(self):
@@ -80,14 +81,15 @@ def _update_resistance_prediction(self, allele_name, variant_or_gene):
             drugs = self._get_drugs(name)
             resistance_prediction = self._resistance_prediction(
                 variant_or_gene, variant_or_gene_names)
+
             for drug in drugs:
                 current_resistance_prediction = self.resistance_predictions[
                     drug]["predict"]
                 assert resistance_prediction is not None
-                if current_resistance_prediction in ["I", "N"]:
+                if current_resistance_prediction == "N":
                     self.resistance_predictions[drug][
-                        "predict"] = resistance_prediction
-                elif current_resistance_prediction == "S":
+                            "predict"] = resistance_prediction
+                elif current_resistance_prediction in ["I", "S"]:
                     if resistance_prediction in ["r", "R"]:
                         self.resistance_predictions[drug][
                             "predict"] = resistance_prediction
@@ -146,26 +148,26 @@ def _get_drugs(self, name, lower=False):
     def _resistance_prediction(self, variant_or_gene, names):
         __resistance_prediction = None
         if sum(variant_or_gene.genotype) == 2:
-            if (variant_or_gene.is_filtered() and self.ignore_filtered) or depth_on_allele(variant_or_gene) < self.depth_threshold:
-                __resistance_prediction =  "I"            
+            if (variant_or_gene.is_filtered() and self.ignore_filtered) or depth_on_alternate(variant_or_gene) < self.depth_threshold:
+                __resistance_prediction =  "N"            
             elif self._coverage_greater_than_threshold(variant_or_gene, names):
                 __resistance_prediction =  "R"
             else:
                 __resistance_prediction =  "S"
         elif sum(variant_or_gene.genotype) == 1:
-            if (variant_or_gene.is_filtered() and self.ignore_filtered) or depth_on_allele(variant_or_gene) < self.depth_threshold:
-                __resistance_prediction =  "I"             
+            if (variant_or_gene.is_filtered() and self.ignore_filtered) or depth_on_alternate(variant_or_gene) < self.depth_threshold:
+                __resistance_prediction =  "N"             
             elif self._coverage_greater_than_threshold(variant_or_gene, names):
                 __resistance_prediction =  "r"
             else:
                 __resistance_prediction =  "S"
         elif sum(variant_or_gene.genotype) == 0:
-            if depth_on_allele(variant_or_gene) < self.depth_threshold:
-                __resistance_prediction = "I"
-            else:
-                __resistance_prediction =  "S"
+            # if depth_on_allele(variant_or_gene) < self.depth_threshold:
+                # __resistance_prediction = "I"
+            # else:
+            __resistance_prediction =  "S"
         else:
-            __resistance_prediction =  "I"
+            __resistance_prediction =  "N"
         return __resistance_prediction
 
     def _coverage_greater_than_threshold(self, variant_or_gene, names):

scripts/compare.py

@@ -125,6 +125,11 @@ def IS(self):
         """Pheno is S predict is inconclusive """
         return self.count_comparision.get('IS', 0)
 
+    @property
+    def I(self):
+        """Inconcusinve or Null calls"""
+        return self.IR + self.IS
+
     @property
     def FP(self):
         return self.count_comparision.get('FP', 0)
@@ -153,6 +158,10 @@ def VME_str(self):
     def ME_str(self):
         return "%s%% (%s%%-%s%%)" % (self.ME, self.ME_LB, self.ME_UB)
 
+    @property
+    def MinorE_str(self):
+        return "%s%% (%s%%-%s%%)" % (self.MinorE, self.MinorE_LB, self.MinorE_UB)
+
     @property
     def sensitivity_str(self):
         return "%s%% (%s%%-%s%%)" % (self.sensitivity,
@@ -165,11 +174,11 @@ def specificity_str(self):
 
     @property
     def FN_str(self):
-        return "%s (%s)" % (self.FN, self.P)
+        return "%s:%s (%s)" % (self.FN, self.IR, self.P)
 
     @property
     def FP_str(self):
-        return "%s (%s)" % (self.FP, self.N)
+        return "%s:%s (%s)" % (self.FP, self.IS, self.N)
 
     def percentage(self, num, denom):
         try:
@@ -185,6 +194,10 @@ def VME(self):
     def ME(self):
         return self.percentage(self.FP, self.N)
 
+    @property
+    def MinorE(self):
+        return self.percentage(self.I, self.total)
+
     @property
     def sensitivity(self):
         return self.percentage(self.TP, self.P)
@@ -209,6 +222,14 @@ def ME_UB(self):
     def ME_LB(self):
         return round(100 * (self.ME_conf[0]), 1)
 
+    @property
+    def MinorE_UB(self):
+        return round(100 * (self.MinorE_conf[1]), 1)
+
+    @property
+    def MinorE_LB(self):
+        return round(100 * (self.MinorE_conf[0]), 1)        
+
     @property
     def sensitivity_UB(self):
         return round(100 * (self.sensitivity_conf[1]), 1)
@@ -273,6 +294,10 @@ def VME_conf(self):
     def ME_conf(self):
         return self.binom_interval(success=self.FP, total=self.N)
 
+    @property
+    def MinorE_conf(self):
+        return self.binom_interval(success=self.I, total=self.N)
+
     @property
     def sensitivity_conf(self):
         return self.binom_interval(success=self.TP, total=self.P)
@@ -291,7 +316,7 @@ def binom_interval(self, success, total, confint=0.95):
 
     @property
     def row_long_header(self):
-        header = ["Total", "TP", "FP", "P", "TN", "FN", "N", "VME", "ME", "Sens", "Spec",
+        header = ["Total", "TP", "FP", "P", "TN", "FN", "N", "VME", "ME",  "Sens", "Spec",
                   "VME_LB", "VME_UB", "ME_LB", "ME_UB", "Sens_LB", "Sens_UB", "Spec_LB", "Spec_UB",
                     "PPV", "PPV_LB", "PPV_UB",
                     "NPV", "NPV_LB", "NPV_UB"
@@ -311,17 +336,19 @@ def row_long(self):
     @property
     def row_short(self):
         return [self.num_samples, self.FN_str, self.FP_str,
-                    self.VME_str, self.ME_str, self.sensitivity_str, self.specificity_str,
+                    self.VME_str, self.ME_str, self.MinorE_str, 
+                    self.sensitivity_str, self.specificity_str,
                     self.PPV_str, self.NPV_str]
 
     @property
     def row_short_header(self):
         header = [
             "Total",
-            "FN(R)",
-            "FP(S)",
+            "FN:IR (R)",
+            "FP:IS (S)",
             "VME",
             "ME",
+            "MinorE",
             "sensitivity",
             "specificity",
             "PPV", "NPV"]
@@ -414,15 +441,15 @@ def compare_analysis_to_truth(sample_ids, truth, ana, ana_name = ""):
                 elif ana_drug_predict not in ["R", "NA", "S", "I"]:
                     sys.stderr.write("failed predict check %s %s - %s - %s \n" % (ana_name, sample_id,drug, ana_drug_predict))
                 else:
-                    if truth_drug_predict == "NA" or ana_drug_predict == "NA":
+                    if truth_drug_predict == "NA":
                         compare = "UNKNOWN"
-                    elif ana_drug_predict == "I":
+                    elif ana_drug_predict == "I" or ana_drug_predict == "NA" :
                         if truth_drug_predict == "R":
                             compare = "IR"
                         elif truth_drug_predict == "S":
                             compare = "IS" 
                         else:
-                            compare = "NA"                       
+                            compare = "UNKNOWN"                       
                     elif truth_drug_predict == ana_drug_predict:
                         if truth_drug_predict == "R":
                             compare = "TP"