feat: map primers by bwa and update datavzrd template (#275)

* feat: replace bowtie

* set min score

* fmt

* cleanup

* cleanup

* set primer params

* fmt

workflow/resources/datavzrd/variant-calls-template.datavzrd.yaml

@@ -224,22 +224,24 @@ views:
                 scale: ordinal
           revel:
             optional: true
-            plot: 
+            plot:
               heatmap:
                 scale: "linear"
                 domain: [0.0, 1.0]
                 range:
                   - white
                   - "#ff5555"
           "regex('.+: allele frequency')":
-            plot: 
-              ticks:
+            plot:
+              heatmap:
                 scale: "linear"
-                domain: [0.0, 1.0]
-                aux-domain-columns:
-                  - "regex('.+: allele frequency')"
+                domain: [0.0, 0.05, 1.0]
+                range:
+                  - white
+                  - "#d3e9f8"
+                  - "#1f77b4"
           "regex('.+: read depth')":
-            plot: 
+            plot:
               ticks:
                 scale: "linear"
                 aux-domain-columns:
@@ -252,6 +254,16 @@ views:
                 range:
                   - white
                   - "#1f77b4"
+          "gnomad genome af":
+            optional: true
+            plot:
+              heatmap:
+                scale: "linear"
+                domain: [0.0, 0.01, 1.0]
+                range:
+                  - white
+                  - "#adebad"
+                  - "#29a329"
           gene:
             display-mode: hidden
           feature:
@@ -328,14 +340,16 @@ views:
                   - "#ec9b00"
                   - "#ecca00"
           "regex('.+: allele frequency')":
-            plot: 
-              ticks:
+            plot:
+              heatmap:
                 scale: "linear"
-                domain: [0.0, 1.0]
-                aux-domain-columns:
-                  - "regex('.+: allele frequency')"
+                domain: [0.0, 0.05, 1.0]
+                range:
+                  - white
+                  - "#d3e9f8"
+                  - "#1f77b4"
           "regex('.+: read depth')":
-            plot: 
+            plot:
               ticks:
                 scale: "linear"
                 aux-domain-columns:
@@ -348,6 +362,16 @@ views:
                 range:
                   - white
                   - "#1f77b4"
+          "gnomad genome af":
+            optional: true
+            plot:
+              heatmap:
+                scale: "linear"
+                domain: [0.0, 0.01, 1.0]
+                range:
+                  - white
+                  - "#adebad"
+                  - "#29a329"
           id:
             optional: true
             display-mode: hidden

workflow/rules/common.smk

@@ -927,15 +927,9 @@ def get_dgidb_datasources():
     return ""
 
 
-def get_bowtie_insertsize():
-    if config["primers"]["trimming"].get("library_length", 0) != 0:
-        return "-X {}".format(config["primers"]["trimming"].get("library_length"))
-    return ""
-
-
 def get_filter_params(wc):
     if isinstance(get_panel_primer_input(wc.panel), list):
-        return "-b -f 2"
+        return "-b -F 12"
     return "-b -F 4"
 
 
@@ -945,10 +939,28 @@ def get_single_primer_flag(wc):
     return ""
 
 
-def format_bowtie_primers(wc, primers):
-    if isinstance(primers, list):
-        return "-1 {r1} -2 {r2}".format(r1=primers[0], r2=primers[1])
-    return primers
+def get_shortest_primer_length(primers):
+    primers = primers if isinstance(primers, list) else [primers]
+    # set to 32 to match bwa-mem default value considering offset of 2
+    min_length = 32
+    for primer_file in primers:
+        with open(primer_file, "r") as p:
+            min_primer = min(
+                [len(p.strip()) for i, p in enumerate(p.readlines()) if i % 2 == 1]
+            )
+            min_length = min(min_length, min_primer)
+    return min_length
+
+
+def get_primer_extra(wc, input):
+    extra = rf"-R '@RG\tID:{wc.panel}\tSM:{wc.panel}' -L 100"
+    min_primer_len = get_shortest_primer_length(input.reads)
+    # Check if shortest primer is below default values
+    if min_primer_len < 32:
+        extra += f" -T {min_primer_len - 2}"
+    if min_primer_len < 19:
+        extra += f" -k {min_primer_len}"
+    return extra
 
 
 def get_datavzrd_data(impact="coding"):

workflow/rules/primers.smk

@@ -44,40 +44,22 @@ rule trim_primers:
         "fgbio TrimPrimers -H -i {input.bam} -p {input.primers} -s {params.sort_order} {params.single_primer} -o {output.trimmed} &> {log}"
 
 
-rule bowtie_build:
+rule map_primers:
     input:
-        genome,
-    output:
-        directory("resources/bowtie_build/"),
-    params:
-        prefix=f"resources/bowtie_build/{genome_name}",
-    log:
-        "logs/bowtie/build.log",
-    conda:
-        "../envs/bowtie.yaml"
-    cache: True
-    shell:
-        "mkdir {output} & "
-        "bowtie-build {input} {params.prefix} &> {log}"
-
-
-rule bowtie_map:
-    input:
-        primers=lambda w: get_panel_primer_input(w.panel),
-        idx="resources/bowtie_build",
+        reads=lambda wc: get_panel_primer_input(wc.panel),
+        idx=rules.bwa_index.output,
     output:
         "results/primers/{panel}_primers.bam",
-    params:
-        primers=lambda wc, input: format_bowtie_primers(wc, input.primers),
-        prefix=f"resources/bowtie_build/{genome_name}",
-        insertsize=get_bowtie_insertsize(),
-        primer_format=lambda wc, input: "-f" if input_is_fasta(input.primers) else "",
     log:
-        "logs/bowtie/{panel}_map.log",
-    conda:
-        "../envs/bowtie.yaml"
-    shell:
-        "bowtie {params.primers} -x {params.prefix} {params.insertsize} {params.primer_format} -S | samtools view -b - > {output} 2> {log}"
+        "logs/bwa_mem/{panel}.log",
+    params:
+        extra=lambda wc, input: get_primer_extra(wc, input),
+        sorting="none",  # Can be 'none', 'samtools' or 'picard'.
+        sort_order="queryname",  # Can be 'queryname' or 'coordinate'.
+        sort_extra="",  # Extra args for samtools/picard.
+    threads: 8
+    wrapper:
+        "v2.13.0/bio/bwa/mem"
 
 
 rule filter_unmapped_primers: