diff --git a/workflow/envs/bowtie.yaml b/workflow/envs/bowtie.yaml deleted file mode 100644 index 84c3d3d69..000000000 --- a/workflow/envs/bowtie.yaml +++ /dev/null @@ -1,6 +0,0 @@ -channels: - - conda-forge - - bioconda -dependencies: - - bowtie =1.3 - - samtools =1.12 \ No newline at end of file diff --git a/workflow/envs/bwa.yaml b/workflow/envs/bwa.yaml deleted file mode 100644 index 9e71f2a5a..000000000 --- a/workflow/envs/bwa.yaml +++ /dev/null @@ -1,5 +0,0 @@ -channels: - - conda-forge - - bioconda -dependencies: - - bwa =0.7.17 diff --git a/workflow/resources/datavzrd/variant-calls-template.datavzrd.yaml b/workflow/resources/datavzrd/variant-calls-template.datavzrd.yaml index a526cf894..e31446c37 100644 --- a/workflow/resources/datavzrd/variant-calls-template.datavzrd.yaml +++ b/workflow/resources/datavzrd/variant-calls-template.datavzrd.yaml @@ -224,7 +224,7 @@ views: scale: ordinal revel: optional: true - plot: + plot: heatmap: scale: "linear" domain: [0.0, 1.0] @@ -232,14 +232,16 @@ views: - 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 diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index caf07760e..13c713ae5 100644 --- a/workflow/rules/common.smk +++ b/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"): diff --git a/workflow/rules/primers.smk b/workflow/rules/primers.smk index 1929ef077..56a599672 100644 --- a/workflow/rules/primers.smk +++ b/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: