if config["illumina_se"] is None and config["Haplotyper"] is True:
    rule map_BAM_to_SAM_pe:
        input:
            MAP_DIR + "{sample}.sorted.bam"
        output:
            MAP_SAM + "{sample}.sorted.sam"
        params:
            scripts_dir = config["scripts_dir"]
        shell:
            "samtools view -h -o {output} {input}"


elif config["illumina_pe"] is None and config["Haplotyper"] is True:
    rule map_BAM_to_SAM_se:
        input:
            MAP_DIR + "{sample}.sorted.bam"
        output:
            MAP_SAM + "{sample}.sorted.sam"
        shell:
            "samtools view -h -o {output} {input}"

else:
	pass

##########################################################################################

if config["illumina_se"] is None and config["Haplotyper"] is True:
    rule make_label:
        input:
        	expand(MAP_SAM + "{sample}.sorted.sam", sample = SAMPLES_PE)
        output:
            MICROHAPLOT + "labels.txt"
        params:
            scripts_dir = config["scripts_dir"]
        shell:
            "{params.scripts_dir}/Make_labels.sh {output}"

elif config["illumina_pe"] is None and config["Haplotyper"] is True:
    rule make_label:
        input:
            expand(MAP_SAM + "{sample}.sorted.sam", sample = SAMPLES_SE)
        output:
            MICROHAPLOT + "labels.txt"
        params:
            scripts_dir = config["scripts_dir"]
        shell:
            "{params.scripts_dir}/Make_labels.sh {output}"

else:
	pass

##########################################################################################

if config["Haplotyper"] is True:
    rule make_SAM_link:
        """
        Make symlinks to run "microhaplot"
        """
        input:
            MAP_SAM + "{sample}.sorted.sam"
        output:
            MICROHAPLOT + "{sample}.sorted.sam"
        log:
            HAPLOTYPE_DOC + "make_SAM_link.{sample}.log"
        threads:
            1
        shell:
            "ln -s $(readlink -f {input}) {output} 2> {log} "

else:
    pass 


########################################################################################


if config["Haplotyper"] is True and config["Validation"] is None:
    rule make_VCF_link:
        """
        Make symlinks to run "microhaplot"
        """
        input:
            vcf = FINAL_VARIANTS + "filtered_snps.recode.vcf"
        output:
            vcf = MICROHAPLOT + "snps.vcf"
        log:
            HAPLOTYPE_DOC + "make_VCF_link.log"
        threads:
            1
        shell:
            "ln -s $(readlink -f {input.vcf}) {output.vcf} 2> {log} "

elif config["Haplotyper"] is True and config["Validation"] is not None:
    rule make_VCF_link:
        """
        Make symlinks to run "microhaplot"
        """
        input:
            vcf = FINAL_VARIANTS + "Overlapping_SNPs.vcf"
        output:
            vcf = MICROHAPLOT + "snps.vcf"
        log:
            HAPLOTYPE_DOC + "make_VCF_link.log"
        threads:
            1
        shell:
            "ln -s $(readlink -f {input.vcf}) {output.vcf} 2> {log} "

else:
	pass

rule map_bwa_index_reference:
    input:
        RAW_DIR + "genome.fasta"
    output:
        RAW_DIR + "genome.fasta.bwt"
    log:
        MAP_DOC + "bwa_index_reference.log"
    benchmark:
        MAP_DOC + "bwa_index_reference.json"
    threads:
        4
    shell:
        "bwa index {input} &> {log}"

#########################################################################################################

if config["illumina_se"] is None and config["Flash"] is False:
   rule map_bwa_sample_pe:
        input:
            reference = RAW_DIR + "genome.fasta",
            index = RAW_DIR + "genome.fasta.bwt",
            forward = QC_DIR + "{sample}.final.1.fq.gz",
            reverse = QC_DIR + "{sample}.final.2.fq.gz"
        output:
            temp(MAP_DIR + "{sample}.unsorted.bam")
        params:
            rg="@RG\tID:{sample}\tSM:{sample}\tPL:Illumina",
            bwa_params = config["BWA-Mapping"]["BWA_mem_params"],
            samtools_params = config["Samtools_mapping_params"]["samtools_view_params"]
        log:
            MAP_DOC + "bwa_{sample}_pe.log"
        benchmark:
            MAP_DOC + "bwa_{sample}_pe.json"
        shell:
            "(bwa mem -R '{params.rg}' {params.bwa_params} {input.reference} {input.forward} {input.reverse} | "
            "samtools view {params.samtools_params} - "
            "> {output}) "
            "2> {log}"

elif config["illumina_se"] is None and config["Flash"] is True:
   rule map_bwa_sample_pe_flashed:
        input:
            reference = RAW_DIR + "genome.fasta",
            index = RAW_DIR + "genome.fasta.bwt",
            Flashed_read = FLASHED_READS + "{sample}.flashed.extendedFrags.fastq.gz"
        output:
            temp(MAP_DIR + "{sample}.unsorted.bam")
        params:
            rg="@RG\tID:{sample}\tSM:{sample}\tPL:Illumina",
            bwa_params = config["BWA-Mapping"]["BWA_mem_params"],
            samtools_params = config["Samtools_mapping_params"]["samtools_view_params"]
        log:
            MAP_DOC + "bwa_{sample}.log"
        benchmark:
            MAP_DOC + "bwa_{sample}.json"
        shell:
            "(bwa mem -R '{params.rg}' {params.bwa_params} {input.reference} {input.Flashed_read} | "
            "samtools view {params.samtools_params} - "
            "> {output}) "
            "2> {log}"


elif config["illumina_pe"] is None:
    rule map_bwa_sample_se:
        input:
            reference = RAW_DIR + "genome.fasta",
            index = RAW_DIR + "genome.fasta.bwt",
            single = QC_DIR + "{sample}.se.final.fq.gz"
        output:
            temp(MAP_DIR + "{sample}.unsorted.bam")
        params:
            rg="@RG\tID:{sample}\tSM:{sample}\tPL:Illumina",
            bwa_params = config["BWA-Mapping"]["BWA_mem_params"],
            samtools_params = config["Samtools_mapping_params"]["samtools_view_params"]
        log:
            MAP_DOC + "bwa_{sample}.log"
        benchmark:
            MAP_DOC + "bwa_{sample}.json"
        shell:
            "(bwa mem -R '{params.rg}' {params.bwa_params} {input.reference} {input.single} | "
            "samtools view {params.samtools_params} - "
            "> {output}) "
            "2> {log}"

#########################################################################################################


rule map_sort_sample:
    input:
        MAP_DIR + "{sample}.unsorted.bam"
    output:
        temp(MAP_DIR + "{sample}.Rmdup.bam")
    log:
        MAP_DOC + "sort_{sample}.log"
    benchmark:
        MAP_DOC + "sort_{sample}.json"
    threads:
        4
    shell:
        "samtools sort -n -@ {threads} "
            "-T $(mktemp --dry-run) "
            "-O bam {input} "
        "> {output} "
        "2> {log}"

#########################################################################################################


rule map_fixmate_sample:
    input:
        MAP_DIR + "{sample}.Rmdup.bam"
    output:
        temp(MAP_DIR + "{sample}.fixmate.bam")
    params:
        fixmate_params = config["SamTools_Fixmate"]["fixmate_params"]
    log:
        MAP_DOC + "fixmate_{sample}.log"
    benchmark:
        MAP_DOC + "fixmate_{sample}.json"
    shell:
        "samtools fixmate {params.fixmate_params} {input} {output} "




#########################################################################################################


rule map_sort_2_sample:
    input:
        MAP_DIR + "{sample}.fixmate.bam"
    output:
        temp(MAP_DIR + "{sample}.fixmate_sorted.bam")
    log:
        MAP_DOC + "sort_2_{sample}.log"
    benchmark:
        MAP_DOC + "sort_2_{sample}.json"
    threads:
        4
    shell:
        "samtools sort -@ {threads} "
            "-T $(mktemp --dry-run) "
            "-O bam {input} "
        "> {output} "
        "2> {log}"


#########################################################################################################


rule map_markdup_sample:
    input:
        MAP_DIR + "{sample}.fixmate_sorted.bam"
    output:
        BAM = MAP_DIR + "{sample}.sorted.bam",
        BAI = MAP_DIR + "{sample}.sorted.bam.bai"
    params:
        markdup_params = config["Marking_Duplicates"]["markdup_params"]
    log:
        MAP_DOC + "rmdup_{sample}.log"
    benchmark:
        MAP_DOC + "rmdup_{sample}.json"
    shell:
        """
        samtools markdup {params.markdup_params} {input} {output.BAM}
        samtools index {output.BAM} {output.BAI}
        """

#########################################################################################################


rule map_stats_sample:
    input:
        MAP_DIR + "{sample}.sorted.bam"
    output:
        FINAL_DOCS_GENERALSTATS + "{sample}.generalStats"
    log:
        MAP_DOC + "stats_{sample}.log"
    benchmark:
        MAP_DOC + "stats_{sample}.json"
    shell:
        "samtools stats {input} | "
        "grep ^SN | cut -f 2- > {output}"


#########################################################################################################


rule map_idxstats_sample:
    input:
        MAP_DIR + "{sample}.sorted.bam"
    output:
        idxstats = FINAL_DOCS_IDXSTATS + "{sample}.idxstats",
        flagstat = FINAL_DOCS_FLAGSTAT + "{sample}.flagstat"
    benchmark:
        MAP_DOC + "index_{sample}_idxstats.json"
    shell:
        """
        samtools idxstats {input} > {output.idxstats}
        samtools flagstat {input} > {output.flagstat}
        """


#########################################################################################################



rule map_QualiMap_bam_sample:
    input:
        MAP_DIR + "{sample}.sorted.bam"
    output:
        sample_header = temp(MAP_DIR + "{sample}.header.sam"),
        new_sample_header = temp(MAP_DIR + "{sample}.header.new.sam"),
        new_bam = temp(MAP_DIR + "{sample}.new.bam"),
        new_bam_index = temp(MAP_DIR + "{sample}.new.bam.bai"),
        qualimap_outfile = FINAL_DOCS_QUALIMAP + "{sample}.qualimap"
    benchmark:
        MAP_DOC + "qualimap_{sample}.json"
    shell:
        """
        samtools view -H {input} > {output.sample_header}
        head -n -3 {output.sample_header} > {output.new_sample_header}
        samtools reheader {output.new_sample_header} {input} > {output.new_bam}
        samtools index {output.new_bam} {output.new_bam_index}
        qualimap bamqc -bam {output.new_bam} -outdir {output.qualimap_outfile} -outfile {output.qualimap_outfile} -outformat html
        """

#########################################################################################################

if config["illumina_se"] is None:
    rule map_MultiQC_Samtools_sample_pe:
        input:
            idxstats = expand(FINAL_DOCS_IDXSTATS + "{sample}.idxstats", sample = SAMPLES_PE),
            flagstat = expand(FINAL_DOCS_FLAGSTAT + "{sample}.flagstat", sample = SAMPLES_PE)
        output:
            MultiQC_idxstats = FINAL_DOCS_MULTIQC_IDXSTATS,
            MultiQC_flagstat = FINAL_DOCS_MULTIQC_FLAGSTAT
        shell:
            """
            multiqc -m samtools {input.idxstats} -o {output.MultiQC_idxstats}
            multiqc -m samtools {input.flagstat} -o {output.MultiQC_flagstat}
            """

elif config["illumina_pe"] is None:
    rule map_MultiQC_Samtools_se:
        input:
            idxstats = expand(FINAL_DOCS_IDXSTATS + "{sample}.idxstats", sample = SAMPLES_SE),
            flagstat = expand(FINAL_DOCS_FLAGSTAT + "{sample}.flagstat", sample = SAMPLES_SE)
        output:
            MultiQC_idxstats = FINAL_DOCS_MULTIQC_IDXSTATS,
            MultiQC_flagstat = FINAL_DOCS_MULTIQC_FLAGSTAT
        shell:
            """
            multiqc -m samtools {input.idxstats} -o {output.MultiQC_idxstats}
            multiqc -m samtools {input.flagstat} -o {output.MultiQC_flagstat}
            """

#########################################################################################################

if config["illumina_se"] is None and config["QC_Tests"] is True:
    rule map_MultiQC_Qualimap_sample_pe:
        input:
            QUALIMAP = expand(FINAL_DOCS_QUALIMAP + "{sample}.qualimap", sample = SAMPLES_PE)
        output:
            MULTIQC_QUALIMAP = FINAL_DOCS_MULTIQC_QUALIMAP
        shell:
            """
            multiqc -m qualimap {input.QUALIMAP} -o {output.MULTIQC_QUALIMAP}
            """

elif config["illumina_pe"] is None and config["QC_Tests"] is True:
    rule map_MultiQC_Qualimap_sample_se:
        input:
            QUALIMAP = expand(FINAL_DOCS_QUALIMAP + "{sample}.qualimap", sample = SAMPLES_SE)
        output:
            MULTIQC_QUALIMAP = FINAL_DOCS_MULTIQC_QUALIMAP
        shell:
            """
            multiqc -m qualimap {input.QUALIMAP} -o {output.MULTIQC_QUALIMAP}
            """
else:
	pass

#########################################################################################################

if config["illumina_se"] is None:
    rule map_merge_bam_pe:
        input:
            bam = expand(MAP_DIR + "{sample}.sorted.bam",
                         sample = SAMPLES_PE),
            bai = expand(MAP_DIR + "{sample}.sorted.bam.bai",
                         sample = SAMPLES_PE)
        output:
            BAM_1 = MAP_DIR + "merged.bam",
            BAM_2 = BAM_FILES + "merged.bam",
            BAI_1 = BAM_FILES + "merged.bam.bai",
            BAI_2 = MAP_DIR + "merged.bam.bai",
            MD5_1 = MAP_DIR + "mergedMD5.txt",
            MD5_2 = BAM_FILES + "mergedMD5.txt"
        log:
            MAP_DOC + "merge_bam_pe.log"
        benchmark:
            MAP_DOC + "merge_bam_pe.json"
        shell:
            """
            samtools merge {output.BAM_1} {input.bam} > {log} 2>&1
            cp {output.BAM_1} {output.BAM_2}
            md5sum {output.BAM_1} > {output.MD5_1}
            md5sum {output.BAM_1} > {output.MD5_2}
            samtools index {output.BAM_1} {output.BAI_1}
            cp {output.BAI_1} {output.BAI_2}
            """

elif config["illumina_pe"] is None:
    rule map_merge_bam_se:
        input:
            bam = expand(MAP_DIR + "{sample}.sorted.bam",
                         sample = SAMPLES_SE),
            bai = expand(MAP_DIR + "{sample}.sorted.bam.bai",
                         sample = SAMPLES_SE)
        output:
            BAM_1 = MAP_DIR + "merged.bam",
            BAM_2 = BAM_FILES + "merged.bam",
            BAI_1 = BAM_FILES + "merged.bam.bai",
            BAI_2 = MAP_DIR + "merged.bam.bai",
            MD5_1 = MAP_DIR + "mergedMD5.txt",
            MD5_2 = BAM_FILES + "mergedMD5.txt"
        log:
            MAP_DOC + "merge_bam_se.log"
        benchmark:
            MAP_DOC + "merge_bam_se.json"
        shell:
            """
            samtools merge {output.BAM_1} {input.bam} > {log} 2>&1
            cp {output.BAM_1} {output.BAM_2}
            md5sum {output.BAM_1} > {output.MD5_1}
            md5sum {output.BAM_1} > {output.MD5_2}
            samtools index {output.BAM_1} {output.BAI_1}
            cp {output.BAI_1} {output.BAI_2}
            """

#########################################################################################################


if config["illumina_se"] is None and config["Flash"] is False:
    rule map_samtools_properly_paired_pe:
        input:
            MAP_DIR + "merged.bam"
        output:
            BAM_1 = MAP_DIR + "merged_properly_paired.bam",
            BAM_2 = BAM_FILES + "merged_properly_paired.bam"
        params:
            samtools_params = config["Samtools_mapping_params"]["samtools_view_params"],
            samtools_filtering_params = config["Samtools_read_filtering"]["samtools_view_filtering"]
        log:
            MAP_DOC + "properly_paired_merged.log"
        benchmark:
            MAP_DOC + "properly_paired_merged.json"
        shell:
            """
			samtools view {params.samtools_params} {params.samtools_filtering_params} {input} > {output.BAM_1}
            cp {output.BAM_1} {output.BAM_2}
            """

elif config["illumina_se"] is None and config["Flash"] is True:
    pass


elif config["illumina_pe"] is None:
    pass


#########################################################################################################

if config["illumina_se"] is None:
    rule properly_paired_index_pe:
        input:
            MAP_DIR + "merged_properly_paired.bam"
        output:
            BAI_1 = MAP_DIR + "merged_properly_paired.bam.bai",
            BAI_2 = BAM_FILES + "merged_properly_paired.bam.bai"
        log:
            MAP_DOC + "index_merged_properly_paired.log"
        benchmark:
            MAP_DOC + "index_merged_properly_paired.json"
        shell:
            """
            samtools index {input} {output.BAI_1}
            cp {output.BAI_1} {output.BAI_2}
            """

elif config["illumina_se"] is None and config["Flash"] is True:
    pass

elif config["illumina_pe"] is None:
    pass

#########################################################################################################


if config["illumina_se"] is None:
    rule map:
        input:
            MAP_DIR + "merged.bam.bai",
            expand([FINAL_DOCS_IDXSTATS + "{sample}.idxstats", FINAL_DOCS_FLAGSTAT + "{sample}.flagstat", FINAL_DOCS_GENERALSTATS + "{sample}.generalStats"],
                sample = SAMPLES_PE)

elif config["illumina_pe"] is None:
    rule map:
        input:
            MAP_DIR + "merged.bam.bai",
            expand([FINAL_DOCS_IDXSTATS + "{sample}.idxstats", FINAL_DOCS_FLAGSTAT + "{sample}.flagstat", FINAL_DOCS_GENERALSTATS + "{sample}.generalStats"],
                sample = SAMPLES_SE)

if config["illumina_se"] is None and config["QC_Tests"] is True:
    rule fastqc_raw_pe:
        input:
            RAW_DIR + "{sample}_{pair}.fq.gz",
        output:
            fastQC_raw + "{sample}.raw.{pair}.fastqc"
        shell:
            """
            fastqc --outdir {fastQC_raw} {input} > {output}
            """

elif config["illumina_pe"] is None and config["QC_Tests"] is True:
    rule fastqc_raw_se:
        input:
            RAW_DIR + "{sample}.se.fq.gz"
        output:
            fastQC_raw + "{sample}.raw.fastqc"
        shell:
            """
            fastqc --outdir {fastQC_raw} {input} > {output}
            """
else:
    pass


############################################################################################################
#"""Run trimmomatic in paired end mode to eliminate Illumina adaptors and remove
 #low quality regions and reads."""
############################################################################################################

if config["illumina_se"] is None:
    rule qc_trimmomatic_pe:
        """
        Run trimmomatic in paired end mode to eliminate Illumina adaptors and
        remove low quality regions and reads.
        Inputs _1 and _2 are piped through gzip/pigz.
        Outputs _1 and _2 are piped to gzip/pigz (level 9).
        Outputs _3 and _4 are compressed with the builtin compressor from
        Trimmomatic. Further on they are catted and compressed with gzip/pigz
        (level 9).
        Note: The cut -f 1 -d " " is to remove additional fields in the FASTQ
        header. It is done posterior to the trimming since the output comes
        slower than the input is read.
        Number of threads used:
            4 for trimmomatic
            2 for gzip inputs
            2 for gzip outputs
            Total: 8
        """
        input:
            forward = RAW_DIR + "{sample}_1.fq.gz",
            reverse = RAW_DIR + "{sample}_2.fq.gz"
        output:
            forward     = QC_DIR + "{sample}.final.1.fq.gz",
            reverse     = QC_DIR + "{sample}.final.2.fq.gz",
            unpaired_pe    = protected(QC_DIR + "{sample}.final.unpaired.fq.gz")
        params:
            unpaired_1  = QC_DIR + "{sample}_3.fq.gz",
            unpaired_2  = QC_DIR + "{sample}_4.fq.gz",
            adapter     = lambda wildcards: config["illumina_pe"][wildcards.sample]["adapter"],
            phred       = lambda wildcards: config["illumina_pe"][wildcards.sample]["phred"],
            trimmomatic_params = config["Read Trimming"]["trimmomatic_params"]
        log:
            QC_DOC + "trimmomatic_pe_{sample}.log"
        benchmark:
            QC_DOC + "trimmomatic_pe_{sample}.json"
        threads:
            24 # I've been able to work with pigz and 24 trimmomatic threads.
        shell:
            """
            trimmomatic PE \
                -threads {threads} \
                -{params.phred} \
                {input.forward} {input.reverse} \
                {output.forward} {params.unpaired_1} {output.reverse} {params.unpaired_2} \
                ILLUMINACLIP:{params.adapter}:2:30:10 \
                {params.trimmomatic_params} \
                2> {log}

                zcat {params.unpaired_1} {params.unpaired_2} |
                cut -f 1 -d " " |
                pigz -9 > {output.unpaired_pe}

            rm {params.unpaired_1} {params.unpaired_2}
            """

elif config["illumina_pe"] is None:
    rule qc_trimmomatic_se:
        """
        Run trimmomatic in paired end mode to eliminate Illumina adaptors and
        remove low quality regions and reads.
        Inputs _1 and _2 are piped through gzip/pigz.
        Outputs _1 and _2 are piped to gzip/pigz (level 9).
        Outputs _3 and _4 are compressed with the builtin compressor from
        Trimmomatic. Further on they are catted and compressed with gzip/pigz
        (level 9).
        Note: The cut -f 1 -d " " is to remove additional fields in the FASTQ
        header. It is done posterior to the trimming since the output comes
        slower than the input is read.
        Number of threads used:
            4 for trimmomatic
            2 for gzip inputs
            2 for gzip outputs
            Total: 8
        """
        input:
            single = RAW_DIR + "{sample}.se.fq.gz"
        output:
            single	= QC_DIR + "{sample}.se.final.fq.gz"
        params:
            adapter	= lambda wildcards: config["illumina_se"][wildcards.sample]["adapter"],
            phred	= lambda wildcards: config["illumina_se"][wildcards.sample]["phred"],
            trimmomatic_params = config["Read Trimming"]["trimmomatic_params"]
        log:
            QC_DOC + "trimmomatic_se_{sample}.log"
        benchmark:
            QC_DOC + "trimmomatic_se_{sample}.json"
        threads:
            24 # I've been able to work with pigz and 24 trimmomatic threads.
        shell:
            """
            trimmomatic SE \
            	-threads {threads} \
            	-{params.phred} \
            	{input.single} \
            	{output.single} \
            	ILLUMINACLIP:{params.adapter}:2:30:10 \
            	{params.trimmomatic_params} \
            	2> {log}
            """



############################################################################################################

if config["illumina_se"] is None and config["QC_Tests"] is True:
    rule fastqc_pe:
        input:
            QC_DIR + "{sample}.final.{pair}.fq.gz"
        output:
            fastQC + "{sample}.{pair}"
        shell:
            """
            fastqc --outdir {fastQC} {input} > {output}
            """

elif config["illumina_pe"] is None and config["QC_Tests"] is True:
    rule fastqc_se:
        input:
            QC_DIR + "{sample}.se.final.fq.gz"
        output:
            fastQC + "{sample}.se"
        shell:
            """
            fastqc --outdir {fastQC} {input} > {output}
            """
else:
	pass

############################################################################################################

if config["illumina_se"] is None and config["Flash"] is True:
    rule flash_pe_reads:
        input:
            forward = QC_DIR + "{sample}.final.1.fq.gz",
            reverse = QC_DIR + "{sample}.final.2.fq.gz"
        output:
            name = FLASHED_READS + "{sample}.flashed",
            file_name = FLASHED_READS + "{sample}.flashed.extendedFrags.fastq.gz"
        params:
            scripts_dir = config["scripts_dir"],
            flash_params = config["FLASH"]["FLASH_params"]
        threads:
            4
        shell:
            "flash {params.flash_params} {input.forward} {input.reverse} -o {output.name} -z"
            #"{params.scripts_dir}/flash_reads.sh {input.forward} {input.reverse} {output.name}"
            "> {output.name} "   


elif config["illumina_pe"] is None:
    pass


############################################################################################################

if config["illumina_se"] is None:
    rule qc_results_pe:
        '''Generate only the resulting files, not the reports'''
        input:
            pe_files = expand(
                [QC_DIR + "{sample}.final.{pair}.fq.gz", fastQC_raw + "{sample}.raw.{pair}.fastqc", FLASHED_READS + "{sample}.flashed.extendedFrags.fastq.gz", FLASHED_READS + "{sample}.flashed"],
                sample = SAMPLES_PE,
                pair = "1 2".split()
            )

elif config["illumina_pe"] is None:
    rule qc_results_se:
        '''Generate only the resulting files, not the reports'''
        input:
            se_files = expand(
				[QC_DIR + "{sample}.se.final.fq.gz", fastQC_raw + "{sample}.raw_fastqc.zip"],
                sample = SAMPLES_SE
            )

############################################################################################################



if config["illumina_se"] is None and config["QC_Tests"] is True:
    rule raw_results_fastqc_pe:
        """Checkpoint to generate all the links for raw data"""
        input:
            expand(
                fastQC_raw + "{sample}.raw.{pair}.fastqc.{extension}", 
                sample = SAMPLES_PE,
                pair = "1 2".split(),
                extension = "zip html".split()
                )
elif config["illumina_pe"] is None and config["QC_Tests"] is True:
    rule raw_results_fastqc_se:
        """Checkpoint to generate all the links for raw data"""
        input:
            expand(
                fastQC_raw + "{sample}.raw.fastqc.{extension}",
                sample = SAMPLES_SE,
                extension = "zip html".split()
           )

rule variant_call_split_genome:
    """
    Split the reference genome into separate components for parallel variant calling
    """
    input:
        genome = RAW_DIR + "genome.fasta"
    output:
        expand(
            RAW_DIR + "genome.{index}.fasta",
            index = map(lambda n: "%01d" % n, range(0, int(config["variant_call_params"]["split_genome"])))
            )
    threads:
        1
    params:
        split = config["variant_call_params"]["split_genome"]
    log:
        VARIANT_CALL_DOC + "variant_call_split_fasta.log"
    benchmark:
        VARIANT_CALL_DOC + "variant_call_split_fasta.json"
    shell:
        "pyfasta split -n {params.split} {input} 2> {log}"


rule variant_call_make_bed:
    """
    Generate a BED file for each of the split FASTA files
    """
    input:
        RAW_DIR + "genome.{index}.fasta"
    output:
        VARIANT_CALL_DIR + "genome.{index}.bed"
    params:
        scripts_dir = config["scripts_dir"]
    threads:
        2
    log:
        VARIANT_CALL_DOC + "variant_call_make_bed.{index}.log"
    benchmark:
        VARIANT_CALL_DOC + "variant_call_make_bed.{index}.json"
    shell:
        "{params.scripts_dir}/fa_to_bed.sh {input}"
        "> {output} "
        "2> {log}"


#########################################################################################################

if config["illumina_se"] is None and config["Flash"] is False:
    rule variant_call_subset_bam_pe:
        """
        Generate a BAM file for each of the genome partitions
        """
        input:
            bam = MAP_DIR + "merged_properly_paired.bam",
            bai = MAP_DIR + "merged_properly_paired.bam.bai",
            bed = VARIANT_CALL_DIR + "genome.{index}.bed"
        output:
            bam = VARIANT_CALL_DIR + "merged.{index}.bam",
            bai = VARIANT_CALL_DIR + "merged.{index}.bam.bai"
        params:
            scripts_dir = config["scripts_dir"]
        threads:
            2
        log:
            VARIANT_CALL_DOC + "variant_call_make_bam.log"
        benchmark:
            VARIANT_CALL_DOC + "variant_call_make_bam.json"
        shell:
            "samtools view -L {input.bed} -b {input.bam} "
            "> {output.bam} "
            "2> {log} && "
            "samtools index {output.bam}"

elif config["illumina_se"] is None and config["Flash"] is True:
    rule variant_call_flashed_bam_pe:
        """
        Generate a BAM file for each of the genome partitions
        """
        input:
            bam = MAP_DIR + "merged.bam",
            bai = MAP_DIR + "merged.bam.bai",
            bed = VARIANT_CALL_DIR + "genome.{index}.bed"
        output:
            bam = VARIANT_CALL_DIR + "merged.{index}.bam",
            bai = VARIANT_CALL_DIR + "merged.{index}.bam.bai"
        params:
            scripts_dir = config["scripts_dir"]
        threads:
            2
        log:
            VARIANT_CALL_DOC + "variant_call_make_bam.log"
        benchmark:
            VARIANT_CALL_DOC + "variant_call_make_bam.json"
        shell:
            "samtools view -L {input.bed} -b {input.bam} "
            "> {output.bam} "
            "2> {log} && "
            "samtools index {output.bam}"


elif config["illumina_pe"] is None:
    rule variant_call_subset_bam_se:
        """
        Generate a BAM file for each of the genome partitions
        """
        input:
            bam = MAP_DIR + "merged.bam",
            bai = MAP_DIR + "merged.bam.bai",
            bed = VARIANT_CALL_DIR + "genome.{index}.bed"
        output:
            bam = VARIANT_CALL_DIR + "merged.{index}.bam",
            bai = VARIANT_CALL_DIR + "merged.{index}.bam.bai"
        params:
            scripts_dir = config["scripts_dir"]
        threads:
            2
        log:
            VARIANT_CALL_DOC + "variant_call_make_bam.log"
        benchmark:
            VARIANT_CALL_DOC + "variant_call_make_bam.json"
        shell:
            "samtools view -L {input.bed} -b {input.bam} "
            "> {output.bam} "
            "2> {log} && "
            "samtools index {output.bam}"

#########################################################################################################


rule variant_call_run_freebayes:
    """
    Run FreeBayes
    """
    input:
        bam = VARIANT_CALL_DIR + "merged.{index}.bam",
        bai = VARIANT_CALL_DIR + "merged.{index}.bam.bai",
        fasta = RAW_DIR + "genome.{index}.fasta"
    output:
        VARIANT_CALL_DIR + "snps.{index}.vcf"
    params:
        fb_params = config["FreeBayes_variant_calling"]["freebayes_params"]
    threads:
        4
    log:
        VARIANT_CALL_DOC + "freebayes.{index}.log"
    benchmark:
        VARIANT_CALL_DOC + "freebayes.{index}.json"
    shell:
        "freebayes {params.fb_params} -f {input.fasta} {input.bam} > {output}"


#########################################################################################################


rule variant_call_bgzip_vcf:
    """
    Compress the VCF files with bgzip
    """
    input:
        VARIANT_CALL_DIR + "snps.{index}.vcf"
    output:
        VARIANT_CALL_DIR + "snps.{index}.vcf.gz"
    threads:
        1
    log:
        VARIANT_CALL_DOC + "bgzip_{index}.log"
    benchmark:
        VARIANT_CALL_DOC + "bgzip_{index}.json"
    shell:
        "bgzip -c {input} > {output}"


#########################################################################################################


rule variant_call_index_vcf:
    """
    Index the VCF files with tabix
    """
    input:
        VARIANT_CALL_DIR + "snps.{index}.vcf.gz"
    output:
        VARIANT_CALL_DIR + "snps.{index}.vcf.gz.tbi"
    threads:
        2
    log:
        VARIANT_CALL_DOC + "vcf_index_{index}.log"
    benchmark:
        VARIANT_CALL_DOC + "vcf_index_{index}.json"
    shell:
        "tabix {input}"


#########################################################################################################


rule variant_call_merge_vcf:
    """
    Merge all VCF files into a single file with bcftools merge
    """
    input:
        vcf = expand(
            VARIANT_CALL_DIR + "snps.{index}.vcf.gz",
            index = map(lambda n: "%01d" % n, range(0, int(config["variant_call_params"]["split_genome"])))),
	    index = expand(
            VARIANT_CALL_DIR + "snps.{index}.vcf.gz.tbi",
            index = map(lambda n: "%01d" % n, range(0, int(config["variant_call_params"]["split_genome"]))))
    output:
        VARIANTS_1 = VARIANT_CALL_DIR + "raw_snps.vcf.gz",
        VARIANTS_2 = RAW_VARIANTS + "raw_snps.vcf.gz"
    threads:
        10
    log:
        VARIANT_CALL_DOC + "bcftools_merge.log"
    benchmark:
        VARIANT_CALL_DOC + "bcftools_merge.json"
    shell:
        """
        bcftools concat -Oz {input.vcf} > {output.VARIANTS_1} 2> {log}
        cp {output.VARIANTS_1} {output.VARIANTS_2}
        """

#########################################################################################################


rule variant_call_filter_vcf:
    """
    Basic site quality filtering for the VCF file
    """
    input:
        VARIANT_CALL_DIR + "raw_snps.vcf.gz"
    output:
        FILTERED_VARIANTS_1 = VARIANT_CALL_DIR + "filtered_snps.recode.vcf",
        FILTERED_VARIANTS_2 = FINAL_VARIANTS + "filtered_snps.recode.vcf",
        HARDY_WEINGBERG_VARIANTS = FINAL_VARIANTS + "Hardy_Weinberg_values"
    threads:
        2
    params:
        vcftools_params = config["VCFtools"]["variant_filtering"]
    log:
        VARIANT_CALL_DOC + "variant_call_filter_vcf.log"
    benchmark:
        VARIANT_CALL_DOC + "variant_call_filter_vcf.json"
    shell:
        """
        vcftools --gzvcf {input} {params.vcftools_params} -c > {output.FILTERED_VARIANTS_1}
        vcftools --gzvcf {output.FILTERED_VARIANTS_1} --hardy -c > {output.HARDY_WEINGBERG_VARIANTS}
        cp {output.FILTERED_VARIANTS_1} {output.FILTERED_VARIANTS_2}
        """

#########################################################################################################

if config["illumina_se"] is None and config["Validation"] is not None:
    rule variant_validation_vcf:
        """
        Use a previous variant dataset to compare, extract and store common high-value SNPs between the two datasets
        """
        input:
            FILTERED_VARIANTS = FINAL_VARIANTS + "filtered_snps.recode.vcf",
            INPUT_VCF_FILE = config["Validation"]
        output:
            FINAL_VARIANTS + "Overlapping_SNPs"
        log:
            VARIANT_CALL_DOC + "variant_validation.log"
        benchmark:
            VARIANT_CALL_DOC + "variant_validation.json"
        shell:
            """
            vcftools --vcf {input.FILTERED_VARIANTS} --diff {input.INPUT_VCF_FILE} --diff-site -c > {output}
            """

else:
    pass

#########################################################################################################

if config["illumina_se"] is None and config["Validation"] is not None:
    rule intersect_vcf:
        input:
            FILTERED_VARIANTS = FINAL_VARIANTS + "filtered_snps.recode.vcf",
            INPUT_VCF_FILE = config["Validation"]
        output:
            FINAL_VARIANTS + "Overlapping_SNPs.vcf"
        log:
            VARIANT_CALL_DOC + "bedtools_intersect.log"
        benchmark:
            VARIANT_CALL_DOC + "bedtools_intersect.json"
        shell:
            """
            bedtools intersect -a {input.FILTERED_VARIANTS} -b {input.INPUT_VCF_FILE} -header > {output}
            """

else:
    pass