Snakefile "ymp/rules/bowtie2.rules"
Env(name="bowtie2", base="bioconda", packages=["bowtie2", "samtools"])
BT2IDX_SUFFIXES = "1.bt2 2.bt2 3.bt2 4.bt2 rev.1.bt2 rev.2.bt2".split()
with Stage("index_bowtie2") as S:
S.doc("""
>>> ymp make toy.ref_genome.index_bowtie2
""")
rule bowtie2_index:
message:
"Bowtie2: Indexing {input}"
input:
"{:prev:}/{target}.fasta.gz"
output:
expand("{{:this:}}/{{target}}.{ext}", ext=BT2IDX_SUFFIXES)
params:
bt2_base="{:this:}/{target}"
threads:
8
log:
"{:this:}/{target}.btbuild.log"
benchmark:
"benchmarks/bowtie2_index/{:this:}/{target}.txt"
conda:
"bowtie2"
shell:
"bowtie2-build-s"
" {input}"
" {params.bt2_base}"
" --threads {threads}"
" >& {log}"
rule bowtie2_index_all:
message:
"Bowtie2: Indexing done"
output:
touch("{:this:}/all_targets.stamp")
input:
expand("{{:this:}}/{{:targets:}}.{ext}", ext=BT2IDX_SUFFIXES)
with Stage("map_bowtie2") as S:
S.doc("""
Map reads using Bowtie2
>>> ymp make toy.ref_genome.index_bowtie2.map_bowtie2
>>> ymp make toy.assemble_megahit.index_bowtie2.map_bowtie2
>>> ymp make toy.group_Subject.assemble_megahit.index_bowtie2.map_bowtie2
>>> ymp make mpic.ref_ssu.index_bowtie2.map_bowtie2
""")
rule bowtie2_map:
message:
"Bowtie2: Mapping {input.fq[0]} to {params.bt2_base}"
input:
fq = "{:prev:}/{:target:}.{:pairnames:}.fq.gz",
index = expand("{{:prev:}}/{{:target:}}.{ext}", ext=BT2IDX_SUFFIXES)
output:
bam = temp( "{:this:}/{target}.bam")
log:
"{:this:}/{target}.log"
benchmark:
"benchmarks/bowtie2_map/{:this:}/{target}.txt"
params:
bt2_base = lambda wc, input: input.index[0][:-len(BT2IDX_SUFFIXES[0])-1],
maxins = 800, # default 500
mem = icfg.mem("8g"),
reads = "-1 {input.fq[0]} -2 {input.fq[1]}"
threads:
16
conda:
"bowtie2"
shell:
"bowtie2 "
" -x {params.bt2_base}"
" {params.reads}"
" -X {params.maxins} "
" -p {threads} "
" 2>{log} "
" | samtools view -b -o {output.bam} -"
rule bowtie2_map_SE: # ymp: extends bowtie2_map
input:
fq = ["{:prev:}/{:target:}.{:pairnames[0]:}.fq.gz"]
params:
reads = "-U {input.fq[0]}"
rule bowtie2_all:
message:
"Bowtie2: Completed"
output:
touch("{:this:}/all_targets.stamp")
input:
"{:this:}/{:targets:}.bam"