Rename variant calling config options

Also adds demix option to restrict the maximum depth passed to samtools mpileup

Author: ahmig

.test/targets.yaml

@@ -23,11 +23,9 @@ PLOTS:
 PROBLEMATIC_VCF:
   "https://raw.githubusercontent.com/W-L/ProblematicSites_SARS-CoV2/da322c32004f7b16bdaa6a8ee7fd24d56e79d9dc/problematic_sites_sarsCov2.vcf"
 VC:
-  MAX_DEPTH: 0
   MIN_QUALITY: 0
-  IVAR_QUALITY: 0
-  IVAR_FREQ: 0.05
-  IVAR_DEPTH: 5
+  MIN_FREQ: 0.05
+  MIN_DEPTH: 5
 DEMIX:
   MIN_QUALITY: 0
   COV_CUTOFF: 0

config/README.md

@@ -140,14 +140,14 @@ All of the following variables are pre-defined in [config.yaml](/config/config.y
   - `REPS`: number of replicates.
   - `THRESHOLD`: value cutoff for visualization.
 - `VC`: variant calling configuration:
-  - `MAX_DEPTH`: maximum depth at a position for `samtools mpileup` (option `-d`).
-  - `MIN_QUALITY`: minimum base quality for `samtools mpileup` (option `-Q`).
-  - `IVAR_QUALITY`: minimum base quality for `ivar variants` (option `-q`).
-  - `IVAR_FREQ`: minimum frequency threshold for `ivar variants` (option `-t`).
-  - `IVAR_DEPTH`: minimum read depth for `ivar variants` (option `-m`).
+  - `MIN_QUALITY`: minimum base quality for `ivar variants` (option `-q`).
+  - `MIN_FREQ`: minimum frequency threshold for `ivar variants` (option `-t`).
+  - `MIN_DEPTH`: minimum read depth for `ivar variants` (option `-m`).
+  - `MAX_DEPTH`: maximum read depth for `samtools mpileup` (option `-d`).
 - `DEMIX`: demixing configuration (uses [Freyja](https://github.com/andersen-lab/Freyja), see also [its docs](https://andersen-lab.github.io/Freyja/index.html)):
   - `PATHOGEN`: pathogen of interest for `freyja update` (option `--pathogen`); must be 'SARS-CoV-2'.
   - `MIN_QUALITY`: minimum quality for `freyja variants` (option `--minq`).
+  - `MAX_DEPTH`: maximum read depth for `samtools mpileup` (option `-d`).
   - `COV_CUTOFF`: minimum depth to calculate the reported "coverage" (percent of sites with that depth) for `freyja demix` (option `--covcut`).
   - `MIN_ABUNDANCE`: minimum lineage estimated abundance for `freyja demix` (option `--eps`).
   - `CONFIRMED_ONLY`: exclude unconfirmed lineages in `freyja demix` (option `--confirmedonly`).

config/config.yaml

@@ -13,11 +13,10 @@ SHALRT:
   REPS: 1000
   THRESHOLD: 0.80
 VC:
-  MAX_DEPTH: 0
-  MIN_QUALITY: 0
-  IVAR_QUALITY: 20
-  IVAR_FREQ: 0.05
-  IVAR_DEPTH: 30
+  MIN_QUALITY: 20
+  MIN_FREQ: 0.05
+  MIN_DEPTH: 30
+  MAX_DEPTH: 1000000
 ANNOTATION:
   # see: https://pcingola.github.io/SnpEff/adds/VCFannotationformat_v1.0.pdf
   SNPEFF_COLS:
@@ -52,6 +51,7 @@ ANNOTATION:
 DEMIX:
   PATHOGEN: "SARS-CoV-2"
   MIN_QUALITY: 20
+  MAX_DEPTH: 1000000
   COV_CUTOFF: 30
   MIN_ABUNDANCE: 0.0001
   CONFIRMED_ONLY: false

workflow/rules/demix.smk

@@ -28,6 +28,7 @@ rule demix_preprocessing:
         bam = get_input_bam,
         ref_fasta = lambda wildcards: select_mapping_references_fasta()
     params:
+        max_depth = config["DEMIX"]["MAX_DEPTH"],
         minq = config["DEMIX"]["MIN_QUALITY"],
     output:
         depth_file = OUTDIR/"demixing"/"{sample}/{sample}_depth.txt",
@@ -37,7 +38,7 @@ rule demix_preprocessing:
         ivar = LOGDIR / "demix_preprocessing" / "{sample}_ivar.log.txt",
     shell:
         "set -euo pipefail && "
-        "samtools mpileup -aa -A -d 600000 -Q {params.minq} -q 0 -B -f {input.ref_fasta:q} {input.bam:q} >sample.pileup 2>{log.pileup:q} && "
+        "samtools mpileup -aa -A -d {params.max_depth} -Q {params.minq} -q 0 -B -f {input.ref_fasta:q} {input.bam:q} >sample.pileup 2>{log.pileup:q} && "
         "ivar variants -p variants -q {params.minq} -r {input.ref_fasta:q} >{log.ivar:q} 2>&1 <sample.pileup && "
         "cut -f1-4 sample.pileup >{output.depth_file:q} && "
         "mv variants.tsv {output.variants_file:q}"

workflow/rules/report.smk

@@ -81,7 +81,7 @@ rule extract_genbank_regions:
 rule polymorphic_sites_over_time_data:
     conda: "../envs/renv.yaml"
     params:
-        max_alt_freq = 1.0 - config["VC"]["IVAR_FREQ"],
+        max_alt_freq = 1.0 - config["VC"]["MIN_FREQ"],
     input:
         variants = OUTDIR/f"{OUTPUT_NAME}.variants.tsv",
         metadata = config["METADATA"],
@@ -428,11 +428,11 @@ rule report:
         value      = REPORT_DIR_TABLES/"diversity.json",
         stats_lm   = REPORT_DIR_TABLES/"time_signal.json",
         stats_ml   = REPORT_DIR_TABLES/"context_phylogeny.json",
-        table      = OUTDIR/"summary_table.csv",
+        table      = REPORT_DIR_TABLES/"summary_table.csv",
         sum_nv     = REPORT_DIR_TABLES/"nv_panel_summary.json",
     params:
         workflow_version = get_repo_version(BASE_PATH.as_posix(), __version__),
-        min_ivar_freq = config["VC"]["IVAR_FREQ"],
+        min_ivar_freq = config["VC"]["MIN_FREQ"],
         ufboot_reps = config["UFBOOT"]["REPS"],
         shalrt_reps = config["SHALRT"]["REPS"],
         name = config["OUTPUT_NAME"],

workflow/rules/vaf.smk

@@ -4,11 +4,11 @@ rule snps_to_ancestor:
     shadow: "minimal"
     conda: "../envs/var_calling.yaml"
     params:
-        max_depth = config["VC"]["MAX_DEPTH"],
-        min_quality = config["VC"]["MIN_QUALITY"],
-        ivar_quality = config["VC"]["IVAR_QUALITY"],
-        ivar_freq = config["VC"]["IVAR_FREQ"],
-        ivar_depth = config["VC"]["IVAR_DEPTH"]
+        mpileup_depth = config["VC"]["MAX_DEPTH"],
+        mpileup_quality = 0,
+        ivar_quality = config["VC"]["MIN_QUALITY"],
+        ivar_freq = config["VC"]["MIN_FREQ"],
+        ivar_depth = config["VC"]["MIN_DEPTH"],
     input:
         reference_fasta = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
         bam = get_input_bam,
@@ -35,10 +35,10 @@ rule snps_to_ancestor:
         samtools mpileup \
             -aa \
             --ignore-overlaps \
-            -d {params.max_depth} \
+            -d {params.mpileup_depth} \
             --count-orphans \
             --no-BAQ \
-            -Q {params.min_quality} \
+            -Q {params.mpileup_quality} \
             -f renamed_reference.fasta \
             {input.bam} \
             | ivar variants \
@@ -215,4 +215,4 @@ rule pairwise_trajectory_correlation:
     log:
         LOGDIR / "pairwise_trajectory_correlation" / "log.txt"
     script:
-        "../scripts/report/pairwise_trajectory_correlation.R"
+        "../scripts/pairwise_trajectory_correlation.R"