fix: Split intersection with target region and normalize calls in sep…

…arate rules. (#90)

* Fix: Split intersection with target region and normalize calls in separate rules.

* fix: bedtools outputs vcf, not vcf.gz, add pipe

Co-authored-by: Johannes Köster <johannes.koester@tu-dortmund.de>

---------

Co-authored-by: Johannes Köster <johannes.koester@tu-dortmund.de>

workflow/rules/common.smk

@@ -93,7 +93,7 @@ def get_plot_cov_labels():  # TODO check if ever used anywhere
     def label(name):
         lower, upper = get_cov_interval(name)
         if upper:
-            return f"{lower}-{upper - 1}"
+            return f"{lower}-{upper-1}"
         return f"≥{lower}"
 
     return {name: label(name) for name in low_coverages}
@@ -260,6 +260,13 @@ def get_target_regions(wildcards):
             return []
 
 
+def intersect_calls(wildcards):
+    if get_target_regions(wildcards) == []:
+        return False
+    else:
+        return True
+
+
 def get_target_regions_intersect_statement(wildcards, input):
     if input.target:
         return f"bedtools intersect -a /dev/stdin -b {input.target} |"

workflow/rules/eval.smk

@@ -91,12 +91,30 @@ rule remove_non_pass:
         "v3.3.6/bio/bcftools/view"
 
 
-rule normalize_calls:
+rule intersect_calls_with_target_regions:
     input:
         bcf="results/filtered-variants/{callset}.bcf",
+        regions=get_target_regions,
+    output:
+        pipe("results/normalized-variants/{callset}_intersected.vcf"),
+    log:
+        "logs/intersect-calls/{callset}.log",
+    conda:
+        "../envs/tools.yaml"
+    shell:
+        "(bedtools intersect -b {input.regions} -a "
+        "<(bcftools view {input.bcf}) -wa -f 1.0 -header > {output}) 2> {log}"
+
+
+rule normalize_calls:
+    input:
+        calls=branch(
+            intersect_calls,
+            then="results/normalized-variants/{callset}_intersected.vcf",
+            otherwise="results/filtered-variants/{callset}.bcf",
+        ),
         ref="resources/reference/genome.fasta",
         ref_index="resources/reference/genome.fasta.fai",
-        regions=get_target_regions,
     output:
         "results/normalized-variants/{callset}.vcf.gz",
     params:
@@ -106,9 +124,7 @@ rule normalize_calls:
     conda:
         "../envs/tools.yaml"
     shell:
-        "(bedtools intersect -b {input.regions} -a "
-        "<(bcftools view {input.bcf}) -wa -f 1.0 -header | "
-        "bcftools norm {params.extra} --fasta-ref {input.ref} | "
+        "(bcftools norm {params.extra} --fasta-ref {input.ref} {input.calls} | "
         "bcftools view -Oz > {output}) 2> {log}"