species = {
    'Salmonella-enterica': {
        'readsID': '2492428',
        'reference_url': 'http://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000988525.2_ASM98852v2/GCA_000988525.2_ASM98852v2_genomic.fna.gz'
    },
    'Staphylococcus-aureus': {
       'readsID': '1274026',
       'reference_url': 'http://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000013425.1_ASM1342v1/GCA_000013425.1_ASM1342v1_genomic.fna.gz'
    }
}

THREADS=2
TEMP='./tmp'

FINAL_FILES = [specie+'.vcf' for specie in species]

rule all:
    input: FINAL_FILES

# === DOWNLOAD ===

rule download_sra:
    output: '{specie}.sra'
    log: 'logs/download_sra/{specie}.log'
    benchmark: 'benchmarks/download_sra/{specie}.txt'
    run:
        readsID = species[wildcards.specie]['readsID']
        shell('''
            bionode-ncbi download sra {readsID} > {TEMP} && rm {TEMP};
            cp {readsID}/*.sra {output} && rm -rf {readsID};
        ''')

rule download_reference:
    output: '{specie}.genomic.fna.gz'
    run:
        reference_url = species[wildcards.specie]['reference_url']
        shell('curl -s {reference_url} -o {output}')

# === EXTRACT/DECOMPRESS ===

rule extract_sra:
    input: '{specie}.sra'
    output: ['{specie}_'+num+'.fastq.gz' for num in ['1', '2']]
    shell: 'fastq-dump --split-files --skip-technical --gzip {input}'

rule decompress_reference:
    input: '{specie}.genomic.fna.gz'
    output: '{specie}.genomic.fna'
    shell: 'bgzip -d {input} --stdout > {output}'

# === MAPPING ===

rule index_reference:
    input: '{specie}.genomic.fna.gz'
    output: ['{specie}.genomic.fna.gz.'+suffix for suffix in ['amb', 'ann', 'bwt', 'pac', 'sa']]
    shell: 'bwa index {input}'

rule bwa_mem:
    input:
        ref = '{specie}.genomic.fna.gz',
        sample = ['{specie}_'+num+'.fastq.gz' for num in ['1', '2']],
        index_files = ['{specie}.genomic.fna.gz.'+suffix for suffix in ['amb', 'ann', 'bwt', 'pac', 'sa']]
    log: 'logs/bwa_mem/{specie}.log'
    output: '{specie}.sam'
    threads: THREADS
    # https://bitbucket.org/snakemake/snakemake-wrappers/src/821062b43423a44c384df3faeb5e0714e0f86205/bio/bwa/mem/wrapper.py?at=master&fileviewer=file-view-default
    wrapper: '0.0.8/bio/bwa_mem'

rule bam:
    input: '{specie}.sam'
    output: '{specie}.bam'
    threads: THREADS
    # https://bitbucket.org/snakemake/snakemake-wrappers/src/821062b43423a44c384df3faeb5e0714e0f86205/bio/samtools/sort/wrapper.py?at=master&fileviewer=file-view-default
    wrapper: '0.0.8/bio/samtools_sort'

rule bai:
    input: '{specie}.bam'
    output: '{specie}.bam.bai'
    # https://bitbucket.org/snakemake/snakemake-wrappers/src/821062b43423a44c384df3faeb5e0714e0f86205/bio/samtools/index/wrapper.py?at=master&fileviewer=file-view-default
    # breaks with v1.2, see: https://www.biostars.org/p/102370/
    # wrapper: '0.0.8/bio/samtools_index'
    shell: 'samtools index {input}'

rule call:
    input:
        bam = '{specie}.bam',
        alignment_index = '{specie}.bam.bai',
        reference = '{specie}.genomic.fna'
    output: '{specie}.vcf'
    shell: 'samtools mpileup -uf {input.reference} {input.bam} | bcftools call -c - > {output}'

# === CLEAN ===

rule clean:
    run:
        for specie in species:
            shell('rm -rf {specie}*; rm -rf logs')
