BWA MEM SAMBLASTER¶
Map reads using bwa mem, mark duplicates by samblaster and sort and index by sambamba.
Example¶
This wrapper can be used in the following way:
rule bwa_mem:
input:
reads=["reads/{sample}.1.fastq", "reads/{sample}.2.fastq"],
idx=multiext("genome", ".amb", ".ann", ".bwt", ".pac", ".sa"),
output:
bam="mapped/{sample}.bam",
index="mapped/{sample}.bam.bai",
log:
"logs/bwa_mem_sambamba/{sample}.log",
params:
extra=r"-R '@RG\tID:{sample}\tSM:{sample}'",
sort_extra="", # Extra args for sambamba.
threads: 8
wrapper:
"v2.6.0/bio/bwa/mem-samblaster"
Note that input, output and log file paths can be chosen freely.
When running with
snakemake --use-conda
the software dependencies will be automatically deployed into an isolated environment before execution.
Software dependencies¶
bwa=0.7.17
sambamba=1.0
samblaster=0.1.26
Authors¶
- Christopher Schröder
Code¶
__author__ = "Christopher Schröder"
__copyright__ = "Copyright 2020, Christopher Schröder"
__email__ = "christopher.schroeder@tu-dortmund.de"
__license__ = "MIT"
from os import path
from snakemake.shell import shell
# Extract arguments.
extra = snakemake.params.get("extra", "")
sort_extra = snakemake.params.get("sort_extra", "")
samblaster_extra = snakemake.params.get("samblaster_extra", "")
index = snakemake.input.get("index", "")
if isinstance(index, str):
index = path.splitext(snakemake.input.idx)[0]
else:
index = path.splitext(snakemake.input.idx[0])[0]
log = snakemake.log_fmt_shell(stdout=False, stderr=True)
# Check inputs/arguments.
if not isinstance(snakemake.input.reads, str) and len(snakemake.input.reads) not in {
1,
2,
}:
raise ValueError("input must have 1 (single-end) or " "2 (paired-end) elements")
shell(
"(bwa mem"
" -t {snakemake.threads}"
" {extra}"
" {index}"
" {snakemake.input.reads}"
" | samblaster"
" {samblaster_extra}"
" | sambamba view -S -f bam /dev/stdin"
" -t {snakemake.threads}"
" | sambamba sort /dev/stdin"
" -t {snakemake.threads}"
" -o {snakemake.output.bam}"
" {sort_extra}"
") {log}"
)