From e6217a12782c5edc3b734487c63c379c8ba055ad Mon Sep 17 00:00:00 2001
From: "Jalees A. Nasir" <jalees_nasir@hotmail.com>
Date: Tue, 11 Jun 2024 16:58:23 +0000
Subject: [PATCH] add breseq parameters

---
 Snakefile    | 22 +++++++++++++++++++---
 signalexe.py |  4 ++++
 2 files changed, 23 insertions(+), 3 deletions(-)

diff --git a/Snakefile b/Snakefile
index 1a7e583..5be2bce 100644
--- a/Snakefile
+++ b/Snakefile
@@ -45,7 +45,8 @@ validate(config, 'resources/config.schema.yaml')
 samples = pd.read_table(config['samples'], sep=',')
 validate(samples, 'resources/sample.schema.yaml')
 
-# manual assignment of breseq reference
+# manual assignment of breseq parameters
+# breseq reference, minor variant thresholds
 try:
     if os.path.exists(config['breseq_reference']):
         breseq_ref = config['breseq_reference']
@@ -54,6 +55,20 @@ try:
 except TypeError:
     breseq_ref = ""
 
+if "polymorphism_variant_coverage" in config.keys():
+    breseq_cov = config['polymorphism_variant_coverage']
+    if breseq_cov == "":
+        breseq_cov = 2
+else:
+    breseq_cov = 2
+
+if "polymorphism_frequency" in config.keys():
+    breseq_freq = config['polymorphism_frequency']
+    if breseq_freq == "":
+        breseq_freq = 0.05
+else:
+    breseq_freq = 0.05
+
 # set output directory
 exec_dir = os.getcwd()
 workdir: os.path.abspath(config['result_dir'])
@@ -588,10 +603,12 @@ rule run_breseq:
         "{sn}/benchmarks/{sn}_run_breseq.benchmark.tsv"
     params:
         ref = os.path.join(exec_dir, breseq_ref),
+        cov = breseq_cov,
+        freq = breseq_freq,
         outdir = '{sn}/breseq'
     shell:
         """
-        breseq --reference {params.ref} --num-processors {threads} --polymorphism-prediction --brief-html-output --output {params.outdir} {input} > {log} 2>&1 || touch {output}
+        breseq --reference {params.ref} --num-processors {threads} --polymorphism-prediction --polymorphism-minimum-variant-coverage-each-strand {params.cov} --polymorphism-frequency-cutoff {params.freq} --brief-html-output --output {params.outdir} {input} > {log} 2>&1 || touch {output}
         """
 
 ################## Based on https://github.com/jts/ncov2019-artic-nf/blob/be26baedcc6876a798a599071bb25e0973261861/modules/illumina.nf ##################
@@ -731,7 +748,6 @@ rule run_kraken2:
 
 ##################################  Based on scripts/quast.sh   ####################################
 
-
 rule run_quast:
     threads: 1
     conda: 'conda_envs/assembly_qc.yaml'
diff --git a/signalexe.py b/signalexe.py
index fe05911..d3fdca9 100755
--- a/signalexe.py
+++ b/signalexe.py
@@ -184,6 +184,10 @@ def write_config_file(run_name, config_file, data_directory, opt_tasks):
 run_breseq: {opt_tasks[0]}
 # Used as --reference argument to 'breseq'
 breseq_reference: "{data_directory}/MN908947.3.gbk"
+# Used as --polymorphism-minimum-variant-coverage-each-strand, --polymorphism-frequency-cutoff arguments
+# Parameters needed to determine thresholds for minor variant detection
+polymorphism_variant_coverage: 2
+polymorphism_frequency: 0.05
 
 # run freebayes for variant and consensus calling (as well as ivar)
 run_freebayes: {opt_tasks[1]}