DEV: Add VOG count threshold fix #5

reneo/__main__.py

@@ -161,7 +161,7 @@ def run_options(func):
             "--evalue",
             default=1e-10,
             required=False,
-            help="maximum e-value for phrog annotations",
+            help="maximum e-value for vog annotations",
             type=float,
             show_default=True,
         ),
@@ -173,6 +173,14 @@ def run_options(func):
             type=float,
             show_default=True,
         ),
+        click.option(
+            "--nvogs",
+            default=10,
+            required=False,
+            help="minimum number of vogs to consider a component",
+            type=int,
+            show_default=True,
+        ),
         click.option(
             "--covtol",
             default=100,

reneo/config/config.yaml

@@ -25,6 +25,7 @@ reneo:
     mgfrac: 0.2
     evalue: 1E-10
     hmmscore: 50
+    nvogs: 10
     covtol: 100
     alpha: 1.2
     koverage_args:

reneo/workflow/rules/reneo.smk

@@ -24,6 +24,7 @@ rule run_reneo:
         mgfrac = config['mgfrac'],
         evalue = config['evalue'],
         hmmscore = config['hmmscore'],
+        nvogs = config['nvogs'],
         covtol = config['covtol'],
         alpha = config['alpha'],
         output = RESDIR,

reneo/workflow/scripts/reneo.py

@@ -1360,6 +1360,9 @@ def main(**kwargs):
     kwargs["logger"].info(
         f"Minimum hmm score for vog annotations: {kwargs['hmmscore']}"
     )
+    kwargs["logger"].info(
+        f"Minimum number of vogs to consider a component: {kwargs['nvogs']}"
+    )
     kwargs["logger"].info(
         f"Coverage tolerance for extending subpaths: {kwargs['covtol']}"
     )
@@ -1622,6 +1625,7 @@ def main(**kwargs):
         mgfrac=float(snakemake.params.mgfrac),
         evalue=float(snakemake.params.evalue),
         hmmscore=float(snakemake.params.hmmscore),
+        nvogs = int(snakemake.params.nvogs),
         covtol=float(snakemake.params.covtol),
         alpha=float(snakemake.params.alpha),
         output=snakemake.params.output,

reneo/workflow/scripts/reneo_utils/component_utils.py

@@ -41,7 +41,7 @@ def get_components(**kwargs):
 
             if (
                 vogs_present
-                # and unitig_names[unitig] in circular
+                and len(vogs_found) >= kwargs['nvogs']
                 and kwargs["edges_lengths"][kwargs["unitig_names"][unitig]] > kwargs["minlength"]
             ):
                 pruned_vs[i] = component