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"