allow for new yaml rules

jodyphelan · jodyphelan · commit 92371f4fa568 · 2025-02-14T18:59:35.000+01:00
diff --git a/db/tbdb.rules.yml b/db/tbdb.rules.yml
@@ -0,0 +1,25 @@
+mmpR5_rule:
+  type: epistasisRule
+  source:
+    gene_name: mmpL5
+    type: LoF
+  target:
+    drug: bedaquiline
+    gene_name: mmpR5
+  source_inactivation_freq_cutoff: 100
+  target_escape_freq_cutoff: 50
+
+eis_rule:
+  type: epistasisRule
+  source:
+    gene_name: eis
+    type: LoF
+  target:
+    drug:
+    - kanamycin
+    - amikacin
+    gene_name: eis
+  source_inactivation_freq_cutoff: 100
+  target_escape_freq_cutoff: 50
+  
+
diff --git a/db/tbdb.variables.json b/db/tbdb.variables.json
@@ -43,6 +43,7 @@
     "spoligotype_spacers": "tbdb.spoligotype_spacers.txt",
     "spoligotype_annotations": "tbdb.spoligotype_list.csv",
     "bedmask": "tbdb.mask.bed",
+    "rules": "tbdb.rules.yml",
     "barcode": "tbdb.barcode.bed"
   }
 }
diff --git a/tb-profiler b/tb-profiler
@@ -18,7 +18,7 @@ from tqdm import tqdm
 import logging
 from rich.logging import RichHandler
 from copy import deepcopy
-
+import yaml
 import importlib
 import pkgutil
 
@@ -143,7 +143,11 @@ def main_profile(args):
     for plugin in tbp.ProfilePlugin.__subclasses__():
         plugin().process_variants(args,variants_profile)
 
-    
+
+    rules_file = args.db_dir+"/"+args.conf['files']['rules']
+    for rule in yaml.load(open(rules_file),Loader=yaml.SafeLoader).values():
+        tbp.apply_epistasis_rule(args,variants_profile,rule)
+
     tbp.clean_up_duplicate_annotations(variants_profile)
 
     # Convert variant objects to DrVariant if they cause resistance
@@ -273,6 +277,7 @@ def main_create_db(args):
         "spoligotype_spacers": args.spoligotypes,
         "spoligotype_annotations": args.spoligotype_annotations,
         "bedmask": {"name":args.bedmask,"convert":1},
+        "rules": args.rules,
     }
     if args.barcode:
         extra_files["barcode"] = args.barcode
@@ -584,6 +589,7 @@ parser_sub.add_argument('--spoligotypes',default="spoligotype_spacers.txt",type=
 parser_sub.add_argument('--spoligotype_annotations','--spoligotype-annotations',default="spoligotype_list.csv")
 parser_sub.add_argument('--barcode',default="barcode.bed",type=str,help='A bed file containing lineage barcode SNPs')
 parser_sub.add_argument('--bedmask',default="mask.bed",type=str,help='A bed file containing a list of low-complexity regions')
+parser_sub.add_argument('--rules',default="rules.yml",type=str,help='A yaml file containing rules for the resistance library')
 parser_sub.add_argument('--amplicon_primers','--amplicon-primers',type=str,help='A file containing a list of amplicon primers')
 parser_sub.add_argument('--match_ref','--match-ref',type=str,help='Match the chromosome name to the given fasta file')
 parser_sub.add_argument('--custom',action="store_true",help='Tells the script this is a custom database, this is used to alter the generation of the version definition')
diff --git a/tbprofiler/rules.py b/tbprofiler/rules.py
@@ -40,80 +40,124 @@ def search_variant(variants: List[Variant], **kwargs) -> List[Variant]:
 
     return list(found_variants)
 
-def inactivate_drug_resistance(variant: Variant):
+def inactivate_drug_resistance(variants: List[Variant]):
     """
     Inactivate a drug resistance variant
     """
-    for ann in variant.consequences[0].annotation:
-        if ann['type']=='drug_resistance':
-            ann['type'] = 'inactivated_drug_resistance'
-    for csq in variant.consequences:
-        for ann in csq.annotation:
-            if ann['type']=='drug_resistance':
-                ann['type'] = 'inactivated_drug_resistance'
+    for var in variants:
+        for csq in var.consequences:
+            for ann in csq.annotation:
+                if ann['type']=='drug_resistance':
+                    ann['type'] = 'inactivated_drug_resistance'
             
+class Rule(ProfilePlugin):
+    pass
 
-
-class MmpR5WHORule(ProfilePlugin):
+class epistasisRule(Rule):
     """
-    Epistasis rule for mmpL5/mmpR5
+    Epistasis rule
     """
 
-    def process_variants(self,args,variants: List[Variant]):
+    def process_variants(
+            self, 
+            source:dict, 
+            target:dict, 
+            args: argparse.Namespace,
+            variants: List[Variant], 
+            source_inactivation_freq_cutoff:int=100, 
+            target_escape_freq_cutoff:int=10,
+            **kwargs
+        ):
         """Generic variant processing method"""
 
-        v1 = search_variant(variants,drug='bedaquiline',gene_name='mmpR5')
-        v2 = search_variant(variants,gene_name='mmpL5',type='LoF')
+        source_vars = search_variant(variants,**source)
+        target_vars = search_variant(variants,**target)
 
-        v1_total_freq = math.ceil(sum([x.freq*100 for x in v1]))
-        v1_total_freq = min(v1_total_freq,100)
+        source_vars_total_freq = math.ceil(sum([x.freq*100 for x in source_vars]))
+        source_vars_total_freq = min(source_vars_total_freq,100)
 
-        v2_total_freq = math.ceil(sum([x.freq*100 for x in v2]))
-        v2_total_freq = min(v2_total_freq,100)
+        target_vars_total_freq = math.ceil(sum([x.freq*100 for x in target_vars]))
+        target_vars_total_freq = min(target_vars_total_freq,100)
 
-        if v1 and v2:
-            v1_changes = ", ".join([x.change for x in v1])
-            v2_changes = ", ".join([x.change for x in v2])
-            note = f"Loss of function mutation(s) detected in mmpL5 ({v2_changes}) which may abrogate the effect of the genetically linked mmpR5 mutation(s) ({v1_changes})."
+        print(source_vars_total_freq,target_vars_total_freq)
+        if source_vars and target_vars:
+            source_vars_changes = ", ".join([x.change for x in source_vars])
+            target_vars_changes = ", ".join([x.change for x in target_vars])
+            note = f"Mutation(s) detected in {source_vars[0].gene_name} ({source_vars_changes}) which may abrogate the effect of the genetically linked {target_vars[0].gene_name} mutation(s) ({target_vars_changes})."
             
-            if v2_total_freq==100:
-                inactivate_drug_resistance(v1)
-            freq_diff = v1_total_freq - v2_total_freq
-            if freq_diff > 10:
-                note += f" However, the combined frequency of the mmpR5 mutation(s) is {freq_diff}% higher than the mmpL5 mutation(s), indicating a potential resistant subpopulation. Please consult the raw data for more information."
+            if source_vars_total_freq>=source_inactivation_freq_cutoff:
+                inactivate_drug_resistance(target_vars)
+            freq_diff = target_vars_total_freq - source_vars_total_freq 
+            if freq_diff > target_escape_freq_cutoff:
+                note += f" However, the combined frequency of the {target_vars[0].gene_name} mutation(s) is {freq_diff}% higher than the {source_vars[0].gene_name} mutation(s), indicating a potential resistant subpopulation."
             
+            note += " Please consult the raw data for more information."
             args.notes.append(note)
 
+def apply_epistasis_rule(args: argparse.Namespace, variants: List[Variant], parameters: dict):
+    logging.debug(f"Applying epistasis rule with parameters: {parameters}")
+    epistasisRule().process_variants(args=args,variants=variants,**parameters)
 
-class eisWHORule(ProfilePlugin):
-    """
-    Epistasis rule for mmpL5/mmpR5
-    """
+# class MmpR5WHORule(ProfilePlugin):
+#     """
+#     Epistasis rule for mmpL5/mmpR5
+#     """
 
-    def process_variants(self,args,variants: List[Variant]):
-        """Generic variant processing method"""
+#     def process_variants(self,args,variants: List[Variant]):
+#         """Generic variant processing method"""
+
+#         v1 = search_variant(variants,drug='bedaquiline',gene_name='mmpR5')
+#         v2 = search_variant(variants,gene_name='mmpL5',type='LoF')
+
+#         v1_total_freq = math.ceil(sum([x.freq*100 for x in v1]))
+#         v1_total_freq = min(v1_total_freq,100)
 
-        v1 = search_variant(variants,drug=['kanamycin','amikacin'],gene_name='eis')
-        v2 = search_variant(variants,gene_name='eis',type='LoF')
+#         v2_total_freq = math.ceil(sum([x.freq*100 for x in v2]))
+#         v2_total_freq = min(v2_total_freq,100)
 
-        v1_total_freq = math.ceil(sum([x.freq*100 for x in v1]))
-        v1_total_freq = min(v1_total_freq,100)
+#         if v1 and v2:
+#             v1_changes = ", ".join([x.change for x in v1])
+#             v2_changes = ", ".join([x.change for x in v2])
+#             note = f"Loss of function mutation(s) detected in mmpL5 ({v2_changes}) which may abrogate the effect of the genetically linked mmpR5 mutation(s) ({v1_changes})."
+            
+#             if v2_total_freq==100:
+#                 inactivate_drug_resistance(v1)
+#             freq_diff = v1_total_freq - v2_total_freq
+#             if freq_diff > 10:
+#                 note += f" However, the combined frequency of the mmpR5 mutation(s) is {freq_diff}% higher than the mmpL5 mutation(s), indicating a potential resistant subpopulation. Please consult the raw data for more information."
+            
+#             args.notes.append(note)
+
+
+# class eisWHORule(ProfilePlugin):
+#     """
+#     Epistasis rule for mmpL5/mmpR5
+#     """
 
-        v2_total_freq = math.ceil(sum([x.freq*100 for x in v2]))
-        v2_total_freq = min(v2_total_freq,100)
+#     def process_variants(self,args,variants: List[Variant]):
+#         """Generic variant processing method"""
 
-        if v1 and v2:
-            v1_changes = ", ".join([x.change for x in v1])
-            v2_changes = ", ".join([x.change for x in v2])
-            note = f"Loss of function mutation(s) detected in eis ({v2_changes}) which may abrogate the effect of the genetically linked eis promoter mutation(s) ({v1_changes})."
+#         v1 = search_variant(variants,drug=['kanamycin','amikacin'],gene_name='eis')
+#         v2 = search_variant(variants,gene_name='eis',type='LoF')
 
-            if v2_total_freq==100:
-                inactivate_drug_resistance(v1)
-            freq_diff = v1_total_freq - v2_total_freq
-            if freq_diff > 10:
-                note += f" However, the combined frequency of the eis promoter mutation(s) is {freq_diff}% higher than the eis coding mutation(s), indicating a potential resistant subpopulation. Please consult the raw data for more information."
+#         v1_total_freq = math.ceil(sum([x.freq*100 for x in v1]))
+#         v1_total_freq = min(v1_total_freq,100)
+
+#         v2_total_freq = math.ceil(sum([x.freq*100 for x in v2]))
+#         v2_total_freq = min(v2_total_freq,100)
+
+#         if v1 and v2:
+#             v1_changes = ", ".join([x.change for x in v1])
+#             v2_changes = ", ".join([x.change for x in v2])
+#             note = f"Loss of function mutation(s) detected in eis ({v2_changes}) which may abrogate the effect of the genetically linked eis promoter mutation(s) ({v1_changes})."
+
+#             if v2_total_freq==100:
+#                 inactivate_drug_resistance(v1)
+#             freq_diff = v1_total_freq - v2_total_freq
+#             if freq_diff > 10:
+#                 note += f" However, the combined frequency of the eis promoter mutation(s) is {freq_diff}% higher than the eis coding mutation(s), indicating a potential resistant subpopulation. Please consult the raw data for more information."
             
-            args.notes.append(note)
+#             args.notes.append(note)
 
 
 class SetConfidence(ProfilePlugin):

Original file line number	Diff line number	Diff line change
`@@ -43,6 +43,7 @@`
`43`	`43`	`"spoligotype_spacers": "tbdb.spoligotype_spacers.txt",`
`44`	`44`	`"spoligotype_annotations": "tbdb.spoligotype_list.csv",`
`45`	`45`	`"bedmask": "tbdb.mask.bed",`
	`46`	`+ "rules": "tbdb.rules.yml",`
`46`	`47`	`"barcode": "tbdb.barcode.bed"`
`47`	`48`	`}`
`48`	`49`	`}`