#!/usr/bin/env snakemake
# -*- coding: utf-8 -*-
import os
import glob
import random
from pathlib import Path


from ikiss import IKISS, dico_tool
from ikiss.global_variables import RANDOM_SEED, MAX_CONTIGS_REPORT, MAX_GENES_REPORT, MAX_RULES_REPORT, JACCARD_PATTERN_LIMIT, JACCARD_PAIRS, MAX_PAIRPLOT_VARIABLES, MAX_BIPLOT_LABELS, MAX_MANHATTAN_DIMENSIONS
import click

ikiss_obj = IKISS(dico_tool, workflow=workflow, config=config)
tools_config = ikiss_obj.tools_config
#cluster_config = ikiss_obj.cluster_config

#print(ikiss_obj.export_use_yaml)
# print for debug:
#print(ikiss_obj)
#print(tools_config)
#exit()

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

# Getting paths on usefully variables
output_dir = config['DATA']['OUTPUT']
fastq_dir = config['DATA']['FASTQ']
reference_file = config['PARAMS']['MAPPING_KMERS']['REF']
gff_file = config['PARAMS']['INTERSECT']['GFF']
feature = config['PARAMS']['INTERSECT']['FEATURE']
samples_file = config['PARAMS']['KMERS_MODULE']['SAMPLES_FILE']
pheno_file = config['PARAMS']['LFMM']['PHENOTYPE_FILE']

basename_reference = Path(reference_file).stem

reference_assembly = config['PARAMS']['ASSEMBLY_KMERS']['REF']
basename_reference2 = Path(reference_assembly).stem

FASTQ, = glob_wildcards(f"{config['DATA']['FASTQ']}{{fastq}}{ikiss_obj.fastq_files_ext}")
#SAMPLE, = glob_wildcards(f"{config['DATA']['FASTQ']}{{sample}}_R1{ikiss_obj.fastq_files_ext}")
SAMPLE = list(ikiss_obj.samples.keys())
PHENO = list(ikiss_obj.phenotype.keys())

# check tools version
if not Path(f"{output_dir}versions.csv").exists():
    click.secho("Check if tools are available and their version before run wokflow")
    ikiss_obj.tools_version_to_df(csv_file=f"{ikiss_obj.snakemake_scripts}/report_template/versions.csv",
        active_tools_list=["KMC", "KMERS_GWAS", "BWA-MEM2", "BWA", "FLAGSTATS", "SAMTOOLS", "SEQTK", "BEDTOOLS", "MERGETAGS", "R"],
        output_file=f"{output_dir}versions.csv")

def manipulation_pages():
    """the pages report.py writes in REPORT/MANIPULATION for this run

    Same conditions as the templates it fills: a page is there only when the module that gives it
    its data ran. They are what the rule render_manipulation renders.
    """
    pages = []
    if 'KMERS_MODULE' in ikiss_obj.tools_activated:
        pages.append("jaccard_manipulation")
    if 'SNMF' in ikiss_obj.diversity_method:
        pages.append("snmf_manipulation")
    if 'GENOME_OFFSET' in ikiss_obj.tools_activated:
        pages.append("genomeoffset_manipulation")
    # the pages of the methods draw the first two dimensions, this one draws them all
    if 'MAPPING_KMERS' in ikiss_obj.tools_activated and ikiss_obj.method:
        pages.append("manhattan_manipulation")
    return pages


def output_final(wildcards):
    dico_final = {
        "fastq_table": rules.fastq_stats.output.fastq_table,
        "kmer_module": f"{output_dir}3.TABLE2BED/",
    }
    #if self.config['WORKFLOW']['KMERS_MODULE']) and self.config['WORKFLOW']['INTERSECT']:
    if ('KMERS_MODULE' in ikiss_obj.tools_activated):
        dico_final.update({
        "kmers_table" : rules.kmers_table.output.kmers_table,
        # how the kmers are shared between the samples, and how alike each group is
        "kmers_by_sample" : rules.kmers_stats.output.by_sample,
        "kmers_sharing" : rules.kmers_stats.output.sharing,
        "kmers_pangenome" : rules.kmers_stats.output.pangenome,
        "kmers_jaccard" : rules.kmers_stats.output.jaccard,
        "kmers_jaccard_pairs" : rules.kmers_stats.output.jaccard_pairs,
        # template filled with the paths of the run: the jaccard again, another way, without
        # running iKISS a second time
        "jaccard_template" : f"{ikiss_obj.snakemake_scripts}/report_template/jaccard_template.qmd",
        })
        if not 'PCADAPT' in ikiss_obj.tools_activated and not 'LFMM' in ikiss_obj.tools_activated:
            if 'MAPPING_KMERS' in ikiss_obj.tools_activated:
                dico_final.update({
                    "global_bam" : rules.samtools_merge.output.combined,
                    "global_bam_stats": rules.samtools_merge.output.stats
                    })
                if 'INTERSECT' in ikiss_obj.tools_activated:
                    dico_final.update({
                    "global_matrix_with_annotation": rules.merging_annotations_and_binary_matrix.output.annotations_and_binary_matrix_info
                    })

    if 'SEX_DETECTION' in ikiss_obj.tools_activated:
        dico_final.update({
            "sex_kmers": rules.sex_kmers.output.sex_kmers,
            "contigs_sex_kmers": rules.sex_contigs.output.contigs_sex,
        })
        # the tagged contigs are annotated only when their alignment on the reference and the
        # feature of the GFF both exist
        if 'INTERSECT' in ikiss_obj.tools_activated and config['PARAMS']['ASSEMBLY_KMERS']['MAPPING_CONTIGS']:
            dico_final.update({
                "sex_by_feature": rules.sex_intersect.output.by_feature,
                "sex_intersect_stats": rules.sex_intersect.output.stats,
            })
    if 'GENOME_OFFSET' in ikiss_obj.tools_activated:
        dico_final.update({
            "genome_offset": rules.genome_offset_models.output.offsets,
            "genome_offset_by_model": expand(rules.genome_offset.output.offset, model=list(ikiss_obj.climate_models)),
            # template filled with the paths of the run: the plots of the offset drawn again from
            # the tables, and another K without editing the configuration
            "genomeoffset_template": f"{ikiss_obj.snakemake_scripts}/report_template/genomeoffset_template.qmd",
        })
    if 'SNMF' in ikiss_obj.diversity_method:
        dico_final.update({
            "method_diversity": expand(rules.merge_diversity_method.output.ok, diversity_method = ikiss_obj.diversity_method),
            # read by report.py through params: declared here so that editing the template
            # rebuilds the filled copy of REPORT/MANIPULATION
            "snmf_template": f"{ikiss_obj.snakemake_scripts}/report_template/snmf_template.qmd",
        })
    if 'PCADAPT' in ikiss_obj.method or 'LFMM' in ikiss_obj.method:
        dico_final.update({
        "method_kmers": expand(rules.merge_method.output.outliers_combined, method=ikiss_obj.method),
        })
        if 'MAPPING_KMERS' in ikiss_obj.tools_activated:
            dico_final.update({
                "outliers_and_mapping" : expand(rules.outliers_position.output.csv, method=ikiss_obj.method),
                "stats" : expand(rules.outliers_position.output.stats, method=ikiss_obj.method),
                "bam" : expand(rules.mapping_kmers_outliers.output.sortedbam, method=ikiss_obj.method),
                "by_chrom" : expand(rules.outliers_position.output.by_chrom, method=ikiss_obj.method),
                "manhattan_template": f"{ikiss_obj.snakemake_scripts}/report_template/manhattan_template.qmd",
            })
            if 'INTERSECT' in ikiss_obj.tools_activated:
                dico_final.update({
                    "stats_by_method" : expand(rules.stats_intersect_and_outliers.output.stats_all, method=ikiss_obj.method),
                    "stats_all_kmers": expand(rules.stats_intersect_allkmers.output.stats_all),
                })
        if 'ASSEMBLY_KMERS' in ikiss_obj.tools_activated:
            dico_final.update({
                "outliers_assembly": expand(rules.mergetags.output.assembled_outliers, method=ikiss_obj.method),
                "outliers_csv": expand(rules.mergetags.output.assembled_csv, method=ikiss_obj.method),
            })
            if config['PARAMS']['ASSEMBLY_KMERS']['MAPPING_CONTIGS']:
                dico_final.update({
                    "sortedbam" : expand(rules.mapping_contigs.output.sortedbam, method=ikiss_obj.method),
                })
                if config['PARAMS']['ASSEMBLY_KMERS']['MANHATTAN_CONTIGS']:
                    dico_final.update({
                        "contigs_position": expand(rules.contigs_position.output.csv, method=ikiss_obj.method),
                        "contigs_by_chrom": expand(rules.contigs_position.output.by_chrom, method=ikiss_obj.method),
                    })
                if 'INTERSECT' in ikiss_obj.tools_activated:
                    dico_final.update({
                        "stats_contigs" : expand(rules.intersect_and_contigs.output.stats_contigs, method=ikiss_obj.method),
                })
    return dico_final

rule final:
    input:
        f"{output_dir}REPORT/BOOK/ikiss_report/index.html"

###################### rules
rule kmers_gwas_per_sample:
    """
    kmers_gwas_per_sample automates the process of counting k-mers for each sample, both as canonical and non-canonical, and then combines the results to produce a file with k-mers and their strand information
    """
    threads: 1
    input:
        forw = f"{fastq_dir}{{sample}}_R1{ikiss_obj.fastq_files_ext}"
    params:
        rev = "" if not ikiss_obj.reverse else f"{fastq_dir}{{sample}}_R2{ikiss_obj.fastq_files_ext}",
        name = f"{{sample}}",
        kmer_size = config['PARAMS']['KMERS_MODULE']['KMER_SIZE'],
        dir = f"{output_dir}1.KMERS_MODULE/{{sample}}"
    output:
        kmers_file = f"{output_dir}1.KMERS_MODULE/{{sample}}/{{sample}}_kmers_with_strand"
    log:
        output = f"{output_dir}LOGS/1.KMERS_MODULE/{{sample}}/{{sample}}_KMERS_MODULE.o",
        error = f"{output_dir}LOGS/1.KMERS_MODULE/{{sample}}/{{sample}}_KMERS_MODULE.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_gwas_per_sample/{{sample}}_KMERS_MODULE.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            forward : {input.forw}
        params:
            reverse : {params.rev}
        output:
            kmers_file: {output.kmers_file}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["KMERS_GWAS"]
    shell:
        """
        cd {params.dir}
        # creating txt files path
        realpath {input.forw} {params.rev} > {params.name}_files.txt
        # calculate canonical and not canonical by each sample
        kmc_v3 -t{threads} -k{params.kmer_size} -ci2 @{params.name}_files.txt {params.name}_kmc3_canon ./ 1> {log.error} 2> {log.output}
        kmc_v3 -t{threads} -k{params.kmer_size} -ci0 -b @{params.name}_files.txt {params.name}_kmc3_all ./ 1>> {log.error} 2>> {log.output}
        #combine 2 runs
        kmers_add_strand_information -c {params.name}_kmc3_canon -n {params.name}_kmc3_all -k {params.kmer_size} -o {params.name}_kmers_with_strand 1>> {log.error} 2>> {log.output}
        # the two KMC databases are read by kmers_add_strand_information and by nothing else:
        # once the kmers with their strand are written they are only taking room on the disk
        rm -f {params.name}_kmc3_canon.kmc_pre {params.name}_kmc3_canon.kmc_suf {params.name}_kmc3_all.kmc_pre {params.name}_kmc3_all.kmc_suf
        """

rule kmers_to_use:
    """
    kmers_to_use prepares and filters a list of k-mers found in a samples
    """
    threads: 1
    input:
        samples_list = expand({rules.kmers_gwas_per_sample.output.kmers_file}, sample=SAMPLE)
    params:
        dir = f"{output_dir}2.KMERS_TABLE/",
        kmer_size = config['PARAMS']['KMERS_MODULE']['KMER_SIZE'],
        mac = config['PARAMS']['KMERS_MODULE']['MAC'],
        p = config['PARAMS']['KMERS_MODULE']['P'],
        kmers_list_path = f"{output_dir}2.KMERS_TABLE/kmers_list_paths.txt"
    output:
        kmers_to_use = f"{output_dir}2.KMERS_TABLE/kmers_to_use",
        # written by the same command, next to the list: how many samples carry each kmer, and
        # the kmers the strand filter rejected. kmers_sharing.py reads both
        shareness = f"{output_dir}2.KMERS_TABLE/kmers_to_use.shareness",
        no_pass = f"{output_dir}2.KMERS_TABLE/kmers_to_use.no_pass_kmers"
    log:
        output = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_TO_USE.o",
        error = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_TO_USE.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_to_use/KMERS_TO_USE.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            samples_list : {input.samples_list}
        output:
            kmers_to_use: {output.kmers_to_use}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        cd {params.dir}
        # create file with paths to kmers_with_strand before merging
        realpath {input.samples_list} > samplesTMP.txt;
        awk -F '/' '{{print $_"\t"$NF"TMP" }}' samplesTMP.txt | sed 's/_kmers_with_strandTMP//g' - > {params.kmers_list_path}
        rm samplesTMP.txt
        # calculate kmers to use
        list_kmers_found_in_multiple_samples -l {params.kmers_list_path} -k {params.kmer_size} --mac {params.mac} -p {params.p} -o {output.kmers_to_use} 1> {log.error} 2> {log.output}
        """

rule kmers_table:
    """
    create_kmers_table build a kmer table using the filtered kmers of several samples
    """
    threads: 1
    input:
        kmers_to_use = rules.kmers_to_use.output.kmers_to_use
    params:
        dir = f"{output_dir}2.KMERS_TABLE/",
        kmer_size = config['PARAMS']['KMERS_MODULE']['KMER_SIZE'],
        kmers_list_path = f"{output_dir}2.KMERS_TABLE/kmers_list_paths.txt",
        kmers_table_name = f"kmers_table"
    output:
        kmers_table = f"{output_dir}2.KMERS_TABLE/kmers_table.table"
    log:
        output = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_TABLE.o",
        error = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_TABLE.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_table/KMERS_TABLE.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            kmers_to_use : {input.kmers_to_use}
        output:
            kmers_table: {output.kmers_table}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["KMERS_GWAS"]
    shell:
        """
        cd {params.dir}
        # create the kmer table
        build_kmers_table -l {params.kmers_list_path} -k {params.kmer_size} -a {input.kmers_to_use} -o {params.kmers_table_name} 1> {log.error} 2> {log.output}
        """

##### BED WC
checkpoint kmers_table_to_bed:
    """
    kmers_table_to_bed converts a kmer table to bed file format
    """
    threads: 1
    input:
        kmers_table = rules.kmers_table.output.kmers_table,
    params:
        dir = f"{output_dir}3.TABLE2BED/",
        kmers_table_name = f"{output_dir}2.KMERS_TABLE/kmers_table",
        pheno = samples_file,
        #pheno = pheno_file,
        kmer_size = config['PARAMS']['KMERS_MODULE']['KMER_SIZE'],
        mac = config['PARAMS']['KMERS_MODULE']['MAC'],
        maf = config['PARAMS']['KMERS_MODULE']['MAF'],
        nb_kmers_in_bed = config['PARAMS']['KMERS_MODULE']['B'],
        bed_name = f"output_file"
    output:
        bed = directory(f"{output_dir}3.TABLE2BED/")
    log:
        output = f"{output_dir}3.TABLE2BED/log/TABLE2BED.o",
        error = f"{output_dir}3.TABLE2BED/log/TABLE2BED.e",
    benchmark:
        f"{output_dir}BENCHMARK/kmers_table_to_bed/TABLE2BED.txt"
    message:
        """
        Launching CHECKPOINT {rule}
        threads: {threads}
        input:
            kmers_table : {input.kmers_table}
        output:
            kmers_table: {output.bed}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["KMERS_GWAS"]
    shell:
        """
        cd {params.dir}
        # obtain the plink binary files (bed, fam, bim) 
        kmers_table_to_bed -t {params.kmers_table_name} -k {params.kmer_size} --maf {params.maf} --mac {params.mac} -b {params.nb_kmers_in_bed} -o {params.bed_name} --phentype_file {params.pheno} 1>> {log.error} 2>> {log.output}
        """

rule kmers_stats:
    """
    kmers_stats says how the kmers are shared between the samples (core, shell, cloud, and how the
    pangenome would grow) and how alike the samples of a group are (Jaccard). Nothing is counted
    again: the sharing comes from the histogram list_kmers_found_in_multiple_samples wrote, the
    Jaccard from the presence/absence of the bed files.
    """
    threads: 1
    input:
        shareness = rules.kmers_to_use.output.shareness,
        bed = rules.kmers_table_to_bed.output.bed
    params:
        table_dir = f"{output_dir}2.KMERS_TABLE/",
        logs = expand(rules.kmers_gwas_per_sample.log.error, sample=SAMPLE),
        mac = config['PARAMS']['KMERS_MODULE']['MAC'],
        soft_core = config['PARAMS']['KMERS_MODULE']['SOFT_CORE'],
        samples_file = samples_file,
        pattern_limit = JACCARD_PATTERN_LIMIT,
        pairs = JACCARD_PAIRS,
        seed = RANDOM_SEED
    output:
        by_sample = f"{output_dir}2.KMERS_TABLE/kmers_by_sample.tsv",
        sharing = f"{output_dir}2.KMERS_TABLE/kmers_sharing.tsv",
        pangenome = f"{output_dir}2.KMERS_TABLE/kmers_pangenome.tsv",
        jaccard = f"{output_dir}2.KMERS_TABLE/kmers_jaccard.tsv",
        jaccard_pairs = f"{output_dir}2.KMERS_TABLE/kmers_jaccard_pairs.tsv"
    log:
        output = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_STATS.o",
        error = f"{output_dir}LOGS/2.KMERS_TABLE/KMERS_STATS.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_stats/KMERS_STATS.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            shareness : {input.shareness}
            bed : {input.bed}
        output:
            sharing : {output.sharing}
            jaccard : {output.jaccard}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (python3 {ikiss_obj.snakemake_scripts}/kmers_sharing.py -t {params.table_dir} \
            -l {params.logs} --mac {params.mac} --soft-core {params.soft_core} \
            -o {params.table_dir}
        python3 {ikiss_obj.snakemake_scripts}/kmers_jaccard.py -b {input.bed} \
            -s {params.samples_file} --pattern-limit {params.pattern_limit} \
            --pairs {params.pairs} --seed {params.seed} \
            -o {output.jaccard} -d {output.jaccard_pairs}
        ) 1>{log.output} 2>{log.error}
        """

rule extract_kmers_from_bed:
    """
    extract_kmers_from_bed takes kmers sequences from a bed and creates a fasta format 
    """
    threads: 1
    input:
        bed = f"{output_dir}3.TABLE2BED/{{bed}}.bed"
    params:
        dir = f"{output_dir}4.EXTRACT_FASTA/",
        tmp = f"{output_dir}3.TABLE2BED/{{bed}}.fasta",
    output:
        fasta = f"{output_dir}4.EXTRACT_FASTA/{{bed}}.fasta.gz"
    log:
        output = f"{output_dir}LOGS/4.EXTRACT_FASTA/{{bed}}_EXTRACT_FASTA.o",
        error = f"{output_dir}LOGS/4.EXTRACT_FASTA/{{bed}}_EXTRACT_FASTA.e",
    benchmark:
        f"{output_dir}BENCHMARK/extract_kmers_from_bed/{{bed}}_EXTRACT_FASTA.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
        output:
            fasta: {output.fasta}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        python3 {ikiss_obj.snakemake_scripts}/kmer-bed2fasta.py -b {input.bed} 
        gzip {params.tmp} 
        mv {params.tmp}.gz {output.fasta} ) 1> {log.error} 2> {log.output}
        """

rule index_ref:
    """
    index_ref for bwa-mem2 or bwa 
    """
    threads: 4
    input:
        ref = reference_file
    params:
        ref_dir = f"{output_dir}REF",
        index_options = config['PARAMS']['MAPPING_KMERS']['INDEX_OPTIONS'],
        index_type = f"bwa-mem2 " if config['PARAMS']['MAPPING_KMERS']['MODE'] == "bwa-mem2" else "bwa",
    output:
        # TODO: manage fasta extensions
        new_ref = f"{output_dir}REF/{basename_reference}.fasta",
        index_tag = f"{output_dir}REF/{basename_reference}.fasta.bwt.2bit.64" if config['PARAMS']['MAPPING_KMERS']['MODE'] == "bwa-mem2" else f"{output_dir}REF/{basename_reference}.fasta.sa"
    log:
        output = f"{output_dir}LOGS/{basename_reference}_INDEXING.o",
        error = f"{output_dir}LOGS/{basename_reference}_INDEXING.e",
    benchmark:
        f"{output_dir}BENCHMARK/index_ref/{basename_reference}_INDEXING.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            ref: {input.ref}
        output:
            new_ref : {output.new_ref}
            index_tag : {output.index_tag}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
    shell:
        """
        (mkdir -p {params.ref_dir}
        cd {params.ref_dir}
        cp {reference_file} {output.new_ref} 
        {params.index_type} index {params.index_options} {output.new_ref} ) 1> {log.error} 2> {log.output}
        """

rule index_ref_to_assembly:
    """
    index_ref_to_assembly for bwa-mem2
    """
    threads: 4
    input:
        ref = reference_assembly
    params:
        ref_dir = f"{output_dir}REF_ASSEMBLY",
    output:
        new_ref = f"{output_dir}REF_ASSEMBLY/{basename_reference2}.fasta",
        index_tag = f"{output_dir}REF_ASSEMBLY/{basename_reference2}.fasta.bwt.2bit.64"
    log:
        output = f"{output_dir}LOGS/{basename_reference2}_ASSEMBLY_INDEXING.o",
        error = f"{output_dir}LOGS/{basename_reference2}_ASSEMBLY_INDEXING.e",
    benchmark:
        f"{output_dir}BENCHMARK/index_ref_to_assembly/{basename_reference2}_ASSEMBLY_INDEXING.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            ref: {input.ref}
        output:
            new_ref : {output.new_ref}
            index_tag : {output.index_tag}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
    shell:
        """
        (mkdir -p {params.ref_dir}
        cd {params.ref_dir}
        cp {input.ref} {output.new_ref} 
        bwa-mem2 index {output.new_ref}) 1> {log.error} 2> {log.output}
        """

rule mapping_kmers:
    """
    mapping_kmers
    """
    threads: 4
    input:
        fasta = f"{output_dir}4.EXTRACT_FASTA/{{bed}}.fasta.gz",
        ref = rules.index_ref.output.new_ref,
        index_tag = rules.index_ref.output.index_tag
    params:
        dir = f"{output_dir}8.MAPPING_KMERS/",
        sai = f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}.sai",
        options = config['PARAMS']['MAPPING_KMERS']['OPTIONS'],
        mode = config['PARAMS']['MAPPING_KMERS']['MODE'],
    output:
        # only the stats is kept: the filtered bams are merged by samtools_merge
        sortedbam = temp(f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_sorted.bam"),
        bai = temp(f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_sorted.bam.bai"),
        stats = f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_sorted.bam.stats",
    log:
        output = f"{output_dir}LOGS/8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_MAPPING.o",
        error = f"{output_dir}LOGS/8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_MAPPING.e",
    benchmark:
        f"{output_dir}BENCHMARK/mapping_kmers/{{bed}}_vs_{basename_reference}_MAPPING.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            fasta: {input.fasta}
        output:
            bam : {output.sortedbam}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
        tools_config["ENV-MODULES"]["SAMTOOLS"]
    shell:
        """
        # piped into samtools sort: the unsorted bam is never written
        (cd {params.dir} 
        if [[ {params.mode} == bwa-mem2 ]]; then 
             bwa-mem2 mem {params.options} -t {threads} {input.ref} {input.fasta} | samtools sort -o {output.sortedbam};
        fi
        if [[ {params.mode} == bwa-aln ]]; then
             bwa aln {params.options} -t {threads} {input.ref} {input.fasta} > {params.sai};
             bwa samse {input.ref} {params.sai} {input.fasta} | samtools sort -o {output.sortedbam};
             rm -f {params.sai};
        fi
        samtools index {output.sortedbam}
        samtools stats {output.sortedbam} > {output.sortedbam}.stats) 1> {log.error} 2> {log.output}
        """

rule filter_bam:
    """
    filter_bam
    """
    threads: 4
    input:
        sortedbam = f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_sorted.bam",
    params:
        flag = config['PARAMS']['MAPPING_KMERS']['FILTER_FLAG'],
        qual = config['PARAMS']['MAPPING_KMERS']['FILTER_QUAL'],
        dir = f"{output_dir}8.MAPPING_KMERS/",
    output:
        # merged by samtools_merge, then removed
        filterbam = temp(f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_FMQ.bam"),
    log:
        output = f"{output_dir}LOGS/8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_FMQ.o",
        error = f"{output_dir}LOGS/8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_FMQ.e",
    benchmark:
        f"{output_dir}BENCHMARK/filter_bam/{{bed}}_vs_{basename_reference}_FMQ.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            sortedbam: {input.sortedbam}
        output:
            filterbam : {output.filterbam}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["SAMTOOLS"]
    shell:
        """
        (cd {params.dir}
        samtools view -bh -F {params.flag} -q {params.qual} {input.sortedbam} > {output.filterbam}
        ) 1> {log.error} 2> {log.output}
        """

###################################### BED WC

def aggregate_bim(wildcards):
    """bim files of every bed: one line per kmer, used for the kmers_total stat"""
    bed_dir = checkpoints.kmers_table_to_bed.get(**wildcards).output.bed
    beds = sorted(glob_wildcards(os.path.join(bed_dir, "{bed}.bed")).bed)
    return expand(f"{output_dir}3.TABLE2BED/{{bed}}.bim", bed=beds)

def aggregate_mapping_stats(wildcards):
    """samtools stats of the per-bed mapping of all kmers, used for the kmers_mapped stat"""
    bed_dir = checkpoints.kmers_table_to_bed.get(**wildcards).output.bed
    beds = sorted(glob_wildcards(os.path.join(bed_dir, "{bed}.bed")).bed)
    return expand(f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_sorted.bam.stats", bed=beds)

def aggregate_bam_to_samtools_merge(wildcards):
    checkpoint_output = checkpoints.kmers_table_to_bed.get(**wildcards).output[0]
    bed=glob_wildcards(os.path.join(checkpoint_output,"{bed}.bed")).bed
    return expand(f"{output_dir}8.MAPPING_KMERS/{{bed}}_vs_{basename_reference}_FMQ.bam", bed=bed)

# merge bam wc
rule samtools_merge:
    """
    list kmers positions files
    """
    threads: 1
    input:
        list = aggregate_bam_to_samtools_merge,
    params:
        dir = f"{output_dir}8.MAPPING_KMERS/"
    output:
        combined = f"{output_dir}9.MERGE_BAM/kmer_position_samtools_merge.bam",
        stats = f"{output_dir}9.MERGE_BAM/kmer_position_samtools_merge.stats",
        bam_list = temp(f"{output_dir}9.MERGE_BAM/bam_files.txt")
    log:
        output = f"{output_dir}LOGS/9.MERGE_BAM/KMERPOSITION_SAMTOOLS_MERGE.o",
        error = f"{output_dir}LOGS/9.MERGE_BAM/KMERPOSITION_SAMTOOLS_MERGE.e",
    benchmark:
        f"{output_dir}BENCHMARK/samtools_merge/KMERPOSITION_SAMTOOLS_MERGE.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            list : {input.list}
        output:
            combined: {output.combined}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir} 
        printf '%s\\n' {input.list} > {output.bam_list}
        samtools merge -@ {threads} -b {output.bam_list} {output.combined}
        samtools flagstats -@ {threads} -O tsv {output.combined} > {output.stats}  ) 1> {log.error} 2> {log.output}
        """
###################################################################################"
###################################################################################"

##### SEGMENT WC
checkpoint split_bed:
    """
    from a bed, random kmers in several list before to PCA
    """
    threads: 1
    input:
        bed = f"{output_dir}3.TABLE2BED/{{bed}}.bed",
        fasta = f"{output_dir}4.EXTRACT_FASTA/{{bed}}.fasta.gz",
    params:
        nb_kmers_in_bed = config['PARAMS']['KMERS_MODULE']['SPLIT_LIST_SIZE'],
        name = f"{{bed}}",
        min_lenght =  config['PARAMS']['KMERS_MODULE']['MIN_LIST_SIZE'],
    output:
        dir = directory(f"{output_dir}5.RANGES/{{bed}}")
    log:
        output = f"{output_dir}LOGS/5.RANGES/{{bed}}_RANGES.o",
        error = f"{output_dir}LOGS/5.RANGES/{{bed}}_RANGES.e",
    benchmark:
        f"{output_dir}BENCHMARK/split_bed/{{bed}}_RANGES.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
            fasta: {input.fasta}
        params:
            name : {params.name} 
        output:
            dir : {output.dir}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        VAR=`zcat {input.fasta} | wc -l - | cut -d' ' -f1 `;
        FILE_LENGHT=$((${{VAR}}/2));
        python3 {ikiss_obj.snakemake_scripts}/split_bed.py --list_length {params.nb_kmers_in_bed} --file_length $FILE_LENGHT --min_length {params.min_lenght} --output-name {params.name} --output-dir {output.dir} --seed {RANDOM_SEED} 1> {log.error} 2> {log.output}
        """

#### Segments: all (bed, segment) pairs and one seeded draw shared by PCADAPT, LFMM and SNMF
# Evaluated lazily from input/params functions, i.e. once the checkpoints kmers_table_to_bed and
# split_bed are done (the Snakefile top level runs before them). The draw depends only on the
# segments and RANDOM_SEED, so a relaunch gives the same plots and does not rerun pcadapt/lfmm.
def all_segments():
    """Sorted (bed, segment) pairs produced by the split_bed checkpoints"""
    bed_dir = checkpoints.kmers_table_to_bed.get().output.bed
    beds = sorted(glob_wildcards(os.path.join(bed_dir, "{bed}.bed")).bed)
    pairs = []
    for bed in beds:
        segment_dir = checkpoints.split_bed.get(bed=bed).output.dir
        segments = glob_wildcards(os.path.join(segment_dir, "{segment}.txt")).segment
        pairs += [(bed, segment) for segment in sorted(segments, key=int)]
    return pairs

def drawn_segments():
    """Seeded permutation of all segments: each method takes its NB_PLOTS first ones"""
    pairs = all_segments()
    return random.Random(RANDOM_SEED).sample(pairs, len(pairs))

def checkpoint_outputs(wildcards):
    """Outputs of the checkpoints read by drawn_segments(): snakemake requires them as inputs
    of any rule whose params depend on a checkpoint (plot_flag)"""
    bed_dir = checkpoints.kmers_table_to_bed.get().output.bed
    beds = sorted(glob_wildcards(os.path.join(bed_dir, "{bed}.bed")).bed)
    return [bed_dir] + [checkpoints.split_bed.get(bed=bed).output.dir for bed in beds]

def plot_flag(nb_plots):
    def flag(wildcards):
        return ' --plot T ' if (wildcards.bed, wildcards.segment) in drawn_segments()[:nb_plots] else ' --plot F '
    return flag

rule pcadapt:
    """
    pca using a segment
    """
    threads: 1
    input:
        kmer_list_file = f"{output_dir}5.RANGES/{{bed}}/{{segment}}.txt",
        segments = checkpoint_outputs,
    params:
        k = config['PARAMS']['PCADAPT']['K'],
        bed = f"{output_dir}3.TABLE2BED/{{bed}}.bed",
        bim = f"{output_dir}3.TABLE2BED/{{bed}}.bim",
        fam = f"{output_dir}3.TABLE2BED/{{bed}}.fam",
        segment_file = f"{{bed}}_{{segment}}",
        plotting = plot_flag(ikiss_obj.times_pcadapt)
    output:
        pvalues = f"{output_dir}6.PCADAPT/{{bed}}_{{segment}}_PCADAPT_pvalues.csv",
    log:
        output = f"{output_dir}LOGS/6.PCADAPT/{{bed}}_{{segment}}_PCADAPT.o",
        error = f"{output_dir}LOGS/6.PCADAPT/{{bed}}_{{segment}}_PCADAPT.e",
    benchmark:
        f"{output_dir}BENCHMARK/pcadapt/{{bed}}_{{segment}}_PCADAPT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            kmer_list_file : {input.kmer_list_file}
        params:
            bed : {params.bed}
            k :  {params.k}
        output:
            pvalues: {output.pvalues}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"]
    shell:
        """
        Rscript {ikiss_obj.snakemake_scripts}/pcadapt.R -x {params.k} -e {params.bed} -i {params.bim} -a {params.fam} --bedname {wildcards.bed} --segment {wildcards.segment} -o {output.pvalues} -k {input.kmer_list_file} {params.plotting} 1>{log.output} 2>{log.error}
        """

def aggregate_segments(wildcards):
    return [f"{output_dir}6.{wildcards.method}/{bed}_{segment}_{wildcards.method}_pvalues.csv"
            for bed, segment in all_segments()]

rule merge_method:
    """
    merging under selection kmers detected by pcadapt or lfmm
    """
    threads: 1
    input:
        agg_segments = aggregate_segments
    params:
        method = f"{{method}}",
        dir = f"{output_dir}6.{{method}}/",
        correction = lambda wildcards: config['PARAMS'][wildcards.method]['CORRECTION'],
        alpha = lambda wildcards: config['PARAMS'][wildcards.method]['ALPHA'],
        n_expected = lambda wildcards, input: len(input.agg_segments)
    output:
        pvalues_combined = f"{output_dir}7.MERGED_{{method}}/merged_{{method}}_pvalues.csv",
        outliers_combined = f"{output_dir}7.MERGED_{{method}}/merged_{{method}}_outliers.csv",
        outliers_fasta = f"{output_dir}7.MERGED_{{method}}/merged_{{method}}_outliers.fasta.gz"
    log:
        output = f"{output_dir}LOGS/7.MERGED_{{method}}/MERGED_{{method}}.o",
        error = f"{output_dir}LOGS/7.MERGED_{{method}}/MERGED_{{method}}.e"
    benchmark:
        f"{output_dir}BENCHMARK/merge_method/MERGED_{{method}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            pvalues of {params.n_expected} segments in {params.dir}
        params:
            correction : {params.correction} < {params.alpha} (genome-wide)
        output:
            pvalues_combined: {output.pvalues_combined}
            outliers_combined : {output.outliers_combined}
            outliers_fasta : {output.outliers_fasta}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"]
    shell:
        """
        Rscript {ikiss_obj.snakemake_scripts}/merge_method.R --method {params.method} --input-dir {params.dir} --n-expected {params.n_expected} --correction {params.correction} --alpha {params.alpha} --out-pvalues {output.pvalues_combined} --out-outliers {output.outliers_combined} --out-fasta {output.outliers_fasta} 1>{log.output} 2>{log.error}
        """

# input function for rule aggregate, return paths to all files produced by the checkpoint 'somestep'
def aggregate_segments_snmf(wildcards):
    return [f"{output_dir}6.{wildcards.diversity_method}/{bed}_{segment}_{wildcards.diversity_method}/kmer.geno"
            for bed, segment in drawn_segments()[:ikiss_obj.times_div]]

rule merge_diversity_method:
    """
    launching SNMF in a nb of segments defined by aggregate_segment_snmf
    """
    threads: 1
    input:
        agg_segments = aggregate_segments_snmf
    params:
        method = f"{{diversity_method}}",
        dir = f"{output_dir}6.{{diversity_method}}/"
    output:
        ok = f"{output_dir}6.{{diversity_method}}/OK.txt",
    log:
        output = f"{output_dir}LOGS/6.{{diversity_method}}/{{diversity_method}}.o",
        error = f"{output_dir}LOGS/6.{{diversity_method}}/{{diversity_method}}.e"
    benchmark:
        f"{output_dir}BENCHMARK/merge_diversity_method/{{diversity_method}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            segments : {input.agg_segments}
        output:
            ok: {output.ok}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (if [[ {params.method} == SNMF ]]; then
            mkdir -p {params.dir} ; cd {params.dir} 
            echo 'OK' > {output.ok}
        fi
        ) 2>{log.error}
        """

rule mapping_kmers_outliers:
    """
    mapping_kmers outliers
    """
    threads: 4
    input:
        fasta = rules.merge_method.output.outliers_fasta,
        ref = rules.index_ref.output.new_ref,
        index_tag = rules.index_ref.output.index_tag
    params:
        dir = f"{output_dir}10.OUTLIERS_{{method}}_POSITION",
        sai = f"{output_dir}10.OUTLIERS_{{method}}_POSITION/{{method}}_vs_{basename_reference}.sai",
        options = config['PARAMS']['MAPPING_KMERS']['OPTIONS'],
        mode = config['PARAMS']['MAPPING_KMERS']['MODE'],
    output:
        sortedbam = f"{output_dir}10.OUTLIERS_{{method}}_POSITION/{{method}}_vs_{basename_reference}_sorted.bam",
    log:
        output = f"{output_dir}LOGS/10.OUTLIERS_{{method}}_POSITION/{{method}}_vs_{basename_reference}_MAPPING.o",
        error = f"{output_dir}LOGS/10.OUTLIERS_{{method}}_POSITION/{{method}}_vs_{basename_reference}_MAPPING.e",
    benchmark:
        f"{output_dir}BENCHMARK/mapping_kmers_outliers/{{method}}_vs_{basename_reference}_MAPPING.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            fasta: {input.fasta}
        output:
            bam : {output.sortedbam}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
        tools_config["ENV-MODULES"]["SAMTOOLS"]
    shell:
        """
        # piped into samtools sort: the unsorted bam is never written
        (cd {params.dir} 
        if [[ {params.mode} == bwa-mem2 ]]; then 
             bwa-mem2 mem {params.options} -t {threads} {input.ref} {input.fasta} | samtools sort -o {output.sortedbam};
        fi
        if [[ {params.mode} == bwa-aln ]]; then
             bwa aln {params.options} -t {threads} {input.ref} {input.fasta} > {params.sai};
             bwa samse {input.ref} {params.sai} {input.fasta} | samtools sort -o {output.sortedbam};
             rm -f {params.sai};
        fi
        samtools index {output.sortedbam}
        samtools stats {output.sortedbam} > {output.sortedbam}.stats) 1> {log.error} 2> {log.output}
        """

rule outliers_position:
    """
    position of the outlier kmers on the reference, read from the bam of mapping_kmers_outliers
    """
    threads: 1
    input:
        outliers = rules.merge_method.output.outliers_combined,
        bam = rules.mapping_kmers_outliers.output.sortedbam,
        bim = aggregate_bim,
        mapping_stats = aggregate_mapping_stats,
        merge_flagstats = rules.samtools_merge.output.stats,
    output:
        csv = f"{output_dir}10.OUTLIERS_{{method}}_POSITION/outliers_with_position.csv",
        stats = f"{output_dir}10.OUTLIERS_{{method}}_POSITION/stats.txt",
        by_chrom = f"{output_dir}10.OUTLIERS_{{method}}_POSITION/outliers_by_chrom.tsv",
    params:
        method = f"{{method}}",
        flag = config['PARAMS']['MAPPING_KMERS']['FILTER_FLAG'],
        qual = config['PARAMS']['MAPPING_KMERS']['FILTER_QUAL'],
    log:
        output = f"{output_dir}LOGS/10.OUTLIERS_{{method}}_POSITION/OUTLIERS_POSITION.o",
        error = f"{output_dir}LOGS/10.OUTLIERS_{{method}}_POSITION/OUTLIERS_POSITION.e"
    benchmark:
        f"{output_dir}BENCHMARK/outliers_position/OUTLIERS_{{method}}_POSITION.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            outliers : {input.outliers}
            bam : {input.bam}
        params:
            filters : F{params.flag} Q{params.qual}
        output:
            csv : {output.csv}
            stats : {output.stats}
            by_chrom : {output.by_chrom}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        python3 {ikiss_obj.snakemake_scripts}/outliers_and_position.py --outliers {input.outliers} --bam {input.bam} \
            --bim {input.bim} --mapping-stats {input.mapping_stats} --merge-flagstats {input.merge_flagstats} \
            --flag {params.flag} --qual {params.qual} \
            --output {output.csv} --output-stats {output.stats} --output-by-chrom {output.by_chrom} 1> {log.output} 2> {log.error}
        """

rule extracting_features_from_gff:
    """
    extracting_features_from_gff
    """
    threads: 1
    input:
        gff = gff_file,
    params:
        feature = feature
    output:
        gff_feature = f"{output_dir}12.GFF_FEATURES/extracted.gff",
    log:
        output = f"{output_dir}LOGS/12.GFF_FEATURES/GFF_FEATURES.o",
        error = f"{output_dir}LOGS/12.GFF_FEATURES/GFF_FEATURES.e"
    benchmark:
        f"{output_dir}BENCHMARK/extracting_features_from_gff/GFF_FEATURES.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            gff : {input.gff}
        params:
            feature : {params.feature}
        output:
            gff_feature : {output.gff_feature}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        # the character class rather than the backslash-s of grep: python 3.12 warns about that
        # unknown escape sequence in the shell string, and both match the same thing
        (grep -e "{params.feature}[[:space:]]" {input.gff} > {output.gff_feature} || \
            {{ echo "no line of {input.gff} holds the feature {params.feature} (PARAMS/INTERSECT/FEATURE)" >&2; exit 1; }}
        ) 2>{log.error}
        """

rule kmers_bedtools_intersect:
    """
    intersect bedtools (bam vs gff) 
    """
    threads: 1
    input:
        bam = rules.samtools_merge.output.combined,
        gff = rules.extracting_features_from_gff.output.gff_feature,
    params:
        dir = f"{output_dir}13.KMERS_INTERSECT"
    output:
        bed = f"{output_dir}13.KMERS_INTERSECT/kmers_intersect_annotation.bed"
    log:
        output = f"{output_dir}LOGS/13.KMERS_INTERSECT/BEDTOOLS_INTERSECT.o",
        error = f"{output_dir}LOGS/13.KMERS_INTERSECT/BEDTOOLS_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_bedtools_intersect/BEDTOOLS_INTERSECT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bam : {input.bam}
            gff : {input.gff}
        output:
            bed : {output.bed}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        bedtools intersect -abam {input.bam} -b {input.gff} -bed -header -wb -wa > {output.bed}) 1> {log.error} 2> {log.output}
        """

rule kmers_bedtools_intersect_method:
    """
    intersect bedtools bam from method and annotation
    """
    threads: 1
    input:
        bam = rules.mapping_kmers_outliers.output.sortedbam,
        gff = rules.extracting_features_from_gff.output.gff_feature,
    params:
        dir = f"{output_dir}13.KMERS_INTERSECT_{{method}}"
    output:
        bed = f"{output_dir}13.KMERS_INTERSECT_{{method}}/{{method}}_kmers_intersect_annotation.bed"
    log:
        output = f"{output_dir}LOGS/13.KMERS_INTERSECT_{{method}}/BEDTOOLS_INTERSECT.o",
        error=f"{output_dir}LOGS/13.KMERS_INTERSECT_{{method}}/BEDTOOLS_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/kmers_bedtools_intersect_method/BEDTOOLS_INTERSECT__{{method}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bam : {input.bam}
            gff : {input.gff}
        output:
            bed : {output.bed}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        bedtools intersect -abam {input.bam} -b {input.gff} -bed -header -wb -wa > {output.bed}) 1> {log.error} 2> {log.output}
        """

#### -------------------extraction matrix presence absence

checkpoint extract_fasta_after_intersect_all_kmers:
    """
    extract_fasta_after_intersect_all_kmers
    """
    threads: 1
    input:
        bed = rules.kmers_bedtools_intersect.output.bed,
    params:
        fasta = f"{output_dir}15.KMERS_IN_GENES_MATRIX/kmers_into_feature.fasta",
        nbkmers = f"{config['PARAMS']['KMERS_MODULE']['SPLIT_LIST_SIZE']}"
    output:
        dir = directory(f"{output_dir}15.KMERS_IN_GENES_MATRIX/RANGES",)
    log:
        output = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/EXTRACTING_KMERS_INTO_FEATURE.o",
        error = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/EXTRACTING_KMERS_INTO_FEATURE.e"
    benchmark:
        f"{output_dir}BENCHMARK/extract_fasta_after_intersect_all_kmers/EXTRACTING_KMERS_INTO_FEATURE.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
        output:
            dir :  {output.dir} 
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (mkdir {output.dir} ; cd {output.dir}
        cut -f4 {input.bed} | awk -F '_' '{{print $4}}' - > {params.fasta} 
        split -l {params.nbkmers} {params.fasta} kmers_into_genes_
        for file in kmers_into_genes*; do mv $file $file.txt; done
        #rm {params.fasta}
        ) 2> {log.output}
        """

rule filter_presence_absence_matrix:
    """
    filter_presence_absence_matrix using kmers located only into genes
    """
    threads: 1
    input:
        fasta = f"{output_dir}15.KMERS_IN_GENES_MATRIX/RANGES/{{txt}}.txt",
    params:
        kmers_table = f"{output_dir}2.KMERS_TABLE/kmers_table",
        dir = f"{output_dir}15.KMERS_IN_GENES_MATRIX/MATRIX"
    output:
        matrix = f"{output_dir}15.KMERS_IN_GENES_MATRIX/MATRIX/{{txt}}.matrix"
    log:
        output = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/MATRIX_{{txt}}.o",
        error = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/MATRIX{{txt}}.e"
    benchmark:
        f"{output_dir}BENCHMARK/filter_presence_absence_matrix/KMERS_IN_GENES_MATRIX_{{txt}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            fasta : {input.fasta}
            kmers_table : {params.kmers_table}
        output:
            matrix : {output.matrix}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        filter_kmers -t {params.kmers_table} -k {input.fasta}  -o {output.matrix} ) 1> {log.error} 2> {log.output}
        """

def run_filter_presence_absence(wildcards):
    checkpoint_output = checkpoints.extract_fasta_after_intersect_all_kmers.get(**wildcards).output[0]
    txt=glob_wildcards(os.path.join(checkpoint_output,"{txt}.txt")).txt
    return expand(rules.filter_presence_absence_matrix.output.matrix, txt=txt)

rule aggregate_presence_absence_matrix:
    """
    aggregate_presence_absence_matrix
    """
    threads: 1
    input:
        list_matrix = run_filter_presence_absence
    params:
        dir = f"{output_dir}15.KMERS_IN_GENES_MATRIX"
    output:
        global_matrix = f"{output_dir}15.KMERS_IN_GENES_MATRIX/global_presence_absence_matrix.txt"
    log:
        output = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/AGGREGATED_MATRIX.o",
        error = f"{output_dir}LOGS/15.KMERS_IN_GENES_MATRIX/AGGREGATED_MATRIX.e"
    benchmark:
        f"{output_dir}BENCHMARK/aggregate_presence_absence_matrix/AGGREGATED_MATRIX.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            list_matrix : {input.list_matrix}
        output:
            global_matrix : {output.global_matrix}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        head -n1 {input.list_matrix[0]} > head.txt
        cat head.txt > {output.global_matrix};
        for file in {input.list_matrix}; do
            do
               grep -v ^"kmer" $file >> {output.global_matrix};
            done
        rm head.txt
        ) 2> {log.output}
        """

rule merging_annotations_and_binary_matrix:
    """
     merging presence_absence_matrix with annotations 
    """
    threads: 1
    input:
        matrix = f"{output_dir}15.KMERS_IN_GENES_MATRIX/global_presence_absence_matrix.txt",
        bed = rules.kmers_bedtools_intersect.output.bed
    params:
        dir = f"{output_dir}16.MATRIX_AND_ANNOTATION",
        samples_dico = ikiss_obj.samples
    output:
        annotations_and_binary_matrix_info = f"{output_dir}16.MATRIX_AND_ANNOTATION/annotations_and_binary_matrix_info.csv",
        stats_by_sample = f"{output_dir}16.MATRIX_AND_ANNOTATION/occurrences_by_sample.csv",
        stats= f"{output_dir}16.MATRIX_AND_ANNOTATION/occurrences_by_group.csv"
    log:
        output = f"{output_dir}LOGS/16.MATRIX_AND_ANNOTATION/MATRIX_AND_ANNOTATION.o",
        error = f"{output_dir}LOGS/16.MATRIX_AND_ANNOTATION/MATRIX_AND_ANNOTATION.e"
    benchmark:
        f"{output_dir}BENCHMARK/merging_annotations_and_binary_matrix/MATRIX_AND_ANNOTATION/.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            matrix : {input.matrix} 
            bed : {input.bed}  
        output:
             annotations_and_binary_matrix_info : {output.annotations_and_binary_matrix_info},
             stats : {output.stats}
        log:
            output: {log.output}
            error: {log.error}
        """
    run:
        import pandas as pd
        # importing data
        matrix_df = pd.read_csv(f"{input.matrix}", delimiter='\t', header="infer")
        bed_df = pd.read_csv(f"{input.bed}", delimiter='\t', header=None)
        # cleaning dataframe before merging
        def cleaning_str(seq):
            return seq.split('_')[3]

        bed_df['sequence'] = bed_df[3].apply(cleaning_str)
        bed_df[3]=bed_df['sequence']
        bed_df.rename(columns={3:"kmer"}, inplace=True)
        bed_df.drop(columns=['sequence'], inplace=True)
        # merging dataframes
        results_df = matrix_df.merge(bed_df, on=['kmer'])
        # writing results
        results_df.to_csv(f"{output.annotations_and_binary_matrix_info}", index=False, sep='\t')
        # calculate stats of presence and absence by group given in the samples_file by user
        df = pd.DataFrame(columns=['sample', 'group', 'presence', 'absence'])
        sample = []
        group = []
        presence = []
        absence = []
        for k, j in ikiss_obj.samples.items():
            if k in results_df.columns:
                sample.append(k)
                group.append(j)
                presence.append(len(results_df[results_df[k] == 1]))
                absence.append(len(results_df[results_df[k] == 0]))
        df['sample'] = sample
        df['group'] = group
        df['presence'] = presence
        df['absence'] = absence
        df.to_csv(f"{output.stats_by_sample}", index=False, sep='\t')
        occ_by_group = df.groupby(group).sum()
        occ_by_group.to_csv(f"{output.stats}", index=True, sep='\t')

## =============================== LFMM

rule get_pca_from_phenotype:
    """
    doing a pca analysis using phenotype file with variables done by user
    """
    threads: 6
    input:
        phenotype = pheno_file,
    output:
        pca_variance = f"{output_dir}6.LFMM_PHENO/PCA_from_phenotype.csv",
    params:
        jupyter = f"{output_dir}6.LFMM_PHENO/PCA_from_phenotype.ipynb",
        html = f"{output_dir}6.LFMM_PHENO/PCA_from_phenotype.html",
        # the sex is not a phenotype to reduce: it stays out of the PCA and is added back next to
        # the components, so LFMM tests it and SEX_DETECTION can read its pvalues
        exclude_column = f"-e {config['PARAMS']['SEX_DETECTION']['COLUMN']}" if 'SEX_DETECTION' in ikiss_obj.tools_activated else "",
        max_pairplot = MAX_PAIRPLOT_VARIABLES,
        max_labels = MAX_BIPLOT_LABELS,
    log:
        output = f"{output_dir}LOGS/6.LFMM_PHENO/PCA_FROM_PHENOTYPE_LFMM.o",
        error = f"{output_dir}LOGS/6.LFMM_PHENO/PCA_FROM_PHENOTYPE_LFMM.e",
    benchmark:
        f"{output_dir}BENCHMARK/get_pca_from_phenotype/PCA_FROM_PHENOTYPE_LFMM.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            phenotype: {input.phenotype}
        output:
            pca_variance: {output.pca_variance}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["KMERS_GWAS"]
    shell:
        """
        # the three commands are redirected together: the script only writes the notebook, what
        # fails is its execution, and its traceback used to go to the job output instead of the log
        (python3 {ikiss_obj.snakemake_scripts}/pca_from_phenotype.py -p {input.phenotype} -o {params.jupyter} {params.exclude_column} --max-pairplot {params.max_pairplot} --max-labels {params.max_labels}
        jupyter nbconvert --execute --inplace {params.jupyter} --ExecutePreprocessor.timeout=6000
        jupyter nbconvert --execute {params.jupyter} --no-input --ExecutePreprocessor.timeout=6000 --to=html
        ) 1>{log.output} 2>{log.error}
        """

rule lfmm:
    """
    lfmm using a segment
    """
    threads: 4
    input:
        #kmer_list_file = create_segments,
        kmer_list_file = f"{output_dir}5.RANGES/{{bed}}/{{segment}}.txt",
        phenotype = pheno_file if not config['PARAMS']['LFMM']['PHENOTYPE_PCA_ANALYSIS'] else rules.get_pca_from_phenotype.output.pca_variance,
        segments = checkpoint_outputs,
    params:
        k = config['PARAMS']['LFMM']['K'],
        alpha = config['PARAMS']['LFMM']['ALPHA'],  # threshold line on the plots only
        bed = f"{output_dir}3.TABLE2BED/{{bed}}.bed",
        bim = f"{output_dir}3.TABLE2BED/{{bed}}.bim",
        fam = f"{output_dir}3.TABLE2BED/{{bed}}.fam",
        segment_file = f"{{bed}}_{{segment}}",
        plotting = plot_flag(ikiss_obj.times_lfmm)
    output:
        pvalues = f"{output_dir}6.LFMM/{{bed}}_{{segment}}_LFMM_pvalues.csv",
    log:
        output = f"{output_dir}LOGS/6.LFMM/{{bed}}_{{segment}}_LFMM.o",
        error = f"{output_dir}LOGS/6.LFMM/{{bed}}_{{segment}}_LFMM.e",
    benchmark:
        f"{output_dir}BENCHMARK/lfmm/{{bed}}_{{segment}}_LFMM.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            kmer_list_file : {input.kmer_list_file}
            phenotype : {input.phenotype}
        params:
            bed : {params.bed}
            k: {params.k}
        output:
            pvalues: {output.pvalues}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"]
    shell:
        """
        Rscript {ikiss_obj.snakemake_scripts}/lfmm.R -x {params.k} -e {params.bed} -i {params.bim} -a {params.fam} --bedname {wildcards.bed} --segment {wildcards.segment} -o {output.pvalues} -k {input.kmer_list_file} -p {input.phenotype} -f {params.alpha} {params.plotting} 1>{log.output} 2>{log.error}
        """

################### DIVERSITY MODE
rule snmf:
    """
    SNMF using segments
    """
    threads: 4
    input:
        kmer_list_file = f"{output_dir}5.RANGES/{{bed}}/{{segment}}.txt",
    params:
        dir = f"{output_dir}LOGS/6.SNMF/",
        k_min = config['PARAMS']['SNMF']['K_MIN'],
        k_max= config['PARAMS']['SNMF']['K_MAX'],
        repetitions = config['PARAMS']['SNMF']['REPETITIONS'],
        bed = f"{output_dir}3.TABLE2BED/{{bed}}.bed",
        segment_file = f"{{bed}}_{{segment}}",
        name_project = f"{{bed}}_{{segment}}_SNMF",
    output:
        snmf_out = f"{output_dir}6.SNMF/{{bed}}_{{segment}}_SNMF/kmer.geno",
        #snmf_pdf= f"{output_dir}6.SNMF/{{bed}}_{{segment}}_SNMF/kmer.snmf.pdf"
    log:
        output = f"{output_dir}LOGS/6.SNMF/{{bed}}_{{segment}}_SNMF.o",
        error = f"{output_dir}LOGS/6.SNMF/{{bed}}_{{segment}}_SNMF.e",
    benchmark:
        f"{output_dir}BENCHMARK/snmf/{{bed}}_{{segment}}_SNMF.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            kmer_list_file : {input.kmer_list_file}
        output:
            out: {output.snmf_out}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"]
    shell:
        """
        cd {params.dir}
        Rscript {ikiss_obj.snakemake_scripts}/snmf.R -e {params.bed} -k {input.kmer_list_file} -r {params.repetitions} --kmin {params.k_min} --kmax {params.k_max} -j 4 -o {output.snmf_out} -t {threads} -n {params.name_project} --plot T 1>{log.output} 2>{log.error}
        """



################### ASSEMBLY_KMERS

def outliers_to_mergetags(wildcards):
    if 'LFMM' in wildcards.method :
        return rules.merge_method.output.outliers_combined
    elif 'PCADAPT' in wildcards.method :
        return rules.merge_method.output.outliers_combined

rule mergetags:
    """
    assembly significant kmers obtained by lfmm  by mergeTags
    https://github.com/Transipedia/dekupl-mergeTags
    """
    threads: 1
    input:
        merged_pvalues = outliers_to_mergetags,
    params:
        method = f"{{method}}",
        dir = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}",
        kmer_size = config['PARAMS']['KMERS_MODULE']['KMER_SIZE'],
        min_overlap = config['PARAMS']['ASSEMBLY_KMERS']['OVERLAP_SIZE'],
        contig_size = config['PARAMS']['ASSEMBLY_KMERS']['FILTER_CONTIG_SIZE'],
        tmp4mergetags = temp(f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/tmp4mergetags.csv"),
    output:
        assembled_outliers = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/outliers_{{method}}_mergetags.fasta",
        assembled_csv= f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/outliers_{{method}}_mergetags.csv",
    log:
        output = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/OUTLIERS_MERGETAGS.o",
        error = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/OUTLIERS_MERGETAGS.e",
    benchmark:
        f"{output_dir}BENCHMARK/mergetags/OUTLIERS_{{method}}_MERGETAGS.txt",
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            merged_pvalues : {input.merged_pvalues}
        output:
            assembled_outliers : {output.assembled_outliers}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (
        if [[ {params.method} == LFMM ]]; then
            mkdir -p {params.dir}; cd {params.dir} 
            awk '{{print $1\"\\t"$4\"\\t"$4\"\\t\"$4\"\\t\"$3}}' {input.merged_pvalues} > {params.tmp4mergetags}
            mergeTags -k 31 -m 15 -n  {params.tmp4mergetags} > {output.assembled_csv}
            awk '{{ if (NR>1 && length($2)>={params.contig_size}) print \">mergeTags_\"length($2)\"_\"$3\"\\n\"$2}}' {output.assembled_csv} > {output.assembled_outliers}
            rm {params.tmp4mergetags}
        fi

        if [[ {params.method} == PCADAPT ]]; then
            mkdir -p {params.dir}; cd {params.dir} 
            awk '{{print $1\"\\t"$5\"\\t"$5\"\\t\"$5\"\\t\"$3}}' {input.merged_pvalues} > {params.tmp4mergetags}
            mergeTags -k 31 -m 15 -n  {params.tmp4mergetags} > {output.assembled_csv}
            awk '{{ if (NR>1 && length($2)>={params.contig_size}) print \">mergeTags_\"length($2)\"_\"$3\"\\n\"$2}}' {output.assembled_csv} > {output.assembled_outliers}
            rm {params.tmp4mergetags}
        fi
        )1> {log.error} 2> {log.output}
        """


##################### MAPPING CONTIGS

def fasta_to_map(wildcards):
    if 'LFMM' in wildcards.method :
        #return rules.mergetags_lfmm.output.assembled_outliers
        return f"{output_dir}11.ASSEMBLY_OUTLIERS_LFMM/outliers_LFMM_mergetags.fasta",
    elif 'PCADAPT' in wildcards.method :
        #return rules.mergetags_pcadapt.output.assembled_outliers
        return f"{output_dir}11.ASSEMBLY_OUTLIERS_PCADAPT/outliers_PCADAPT_mergetags.fasta"


def get_ref(wildcards):
    if config['PARAMS']['MAPPING_KMERS']:
        if (basename_reference == basename_reference2):
            if config['PARAMS']['MAPPING_KMERS']['MODE'] == "bwa-mem2" :
                return (rules.index_ref.output.new_ref)
            elif config['PARAMS']['MAPPING_KMERS']['MODE'] == "bwa-aln":
                return (rules.index_ref_to_assembly.output.new_ref)
        else :
            return (rules.index_ref_to_assembly.output.new_ref)
    else :
            return (rules.index_ref_to_assembly.output.new_ref)

rule mapping_contigs:
    """
    mapping contigs vs ref
    """
    threads: 4
    input:
        fasta = fasta_to_map,
        #ref = rules.index_ref_to_assembly.output.new_ref,
        #index_tag = rules.index_ref_to_assembly.output.index_tag
        ref = get_ref
    params:
        dir = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/",
        options = config['PARAMS']['ASSEMBLY_KMERS']['MAPPING_OPTIONS'],
        mode = "bwa-mem2",
    output:
        sortedbam = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/contigs_{{method}}_vs_{basename_reference}.sorted.bam",
    log:
        output = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/contigs_{{method}}_vs_{basename_reference}.MAPPING.o",
        error = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/contigs_{{method}}_vs_{basename_reference}.MAPPING.e",
    benchmark:
        f"{output_dir}BENCHMARK/mapping_contigs/{{method}}_vs_{basename_reference}_MAPPING.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            fasta: {input.fasta}
            ref: {input.ref}
        output:
            bam : {output.sortedbam}
        log:
            output: {log.output}
            error: {log.error} 
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
        tools_config["ENV-MODULES"]["SAMTOOLS"]
    shell:
        """
        # piped into samtools sort: the unsorted bam is never written
        (cd {params.dir} 
        bwa-mem2 mem {params.options} -t {threads} {input.ref} {input.fasta} | samtools sort -o {output.sortedbam};
        samtools index {output.sortedbam}
        samtools idxstats {output.sortedbam} > {output.sortedbam}.idxstats 
        samtools stats {output.sortedbam} > {output.sortedbam}.stats) 1> {log.error} 2> {log.output}
        """

rule contigs_position:
    """
    position of the assembled contigs on the reference and number of outlier kmers they carry,
    for the contig manhattan plot of the report
    """
    threads: 1
    input:
        bam = rules.mapping_contigs.output.sortedbam,
        contigs = rules.mergetags.output.assembled_csv,
    params:
        flag = config['PARAMS']['MAPPING_KMERS']['FILTER_FLAG'],
        qual = config['PARAMS']['MAPPING_KMERS']['FILTER_QUAL'],
    output:
        csv = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/contigs_with_position.tsv",
        by_chrom = f"{output_dir}11.ASSEMBLY_OUTLIERS_{{method}}/contigs_by_chrom.tsv",
    log:
        output = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/CONTIGS_POSITION.o",
        error = f"{output_dir}LOGS/11.ASSEMBLY_OUTLIERS_{{method}}/CONTIGS_POSITION.e"
    benchmark:
        f"{output_dir}BENCHMARK/contigs_position/CONTIGS_{{method}}_POSITION.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bam : {input.bam}
            contigs : {input.contigs}
        output:
            csv : {output.csv}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        python3 {ikiss_obj.snakemake_scripts}/contigs_and_position.py --contigs {input.contigs} --bam {input.bam} --flag {params.flag} --qual {params.qual} --output {output.csv} --output-by-chrom {output.by_chrom} 1>{log.output} 2>{log.error}
        """

rule contigs_bedtools_intersect:
    """
    intersect contigs with bedtools (bam vs gff) 
    """
    threads: 1
    input:
        bam = rules.mapping_contigs.output.sortedbam,
        gff = rules.extracting_features_from_gff.output.gff_feature,
    params:
        dir = f"{output_dir}13.CONTIGS_INTERSECT_{{method}}"
    output:
        bed = f"{output_dir}13.CONTIGS_INTERSECT_{{method}}/{{method}}_contigs_intersect_annotation.bed"
    log:
        output = f"{output_dir}LOGS/13.CONTIGS_INTERSECT_{{method}}/BEDTOOLS_CONTIGS_{{method}}_INTERSECT.o",
        error = f"{output_dir}LOGS/13.CONTIGS_INTERSECT_{{method}}/BEDTOOLS_CONTIGS_{{method}}_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/contigs_bedtools_intersect/BEDTOOLS_CONTIGS_INTERSECT_{{method}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bam : {input.bam}
            gff : {input.gff}
        output:
            bed : {output.bed}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        bedtools intersect -abam {input.bam} -b {input.gff} -bed -header -wb -wa > {output.bed})1> {log.error} 2> {log.output}
        """

##################### INTERSECT CONTIGS

rule intersect_and_contigs:
    """
    stats after intersection (bam vs gff) and contigs {using kmers detected by method}
    """
    threads: 1
    input:
        bed = rules.contigs_bedtools_intersect.output.bed
    params:
        dir = f"{output_dir}14.STATS_INTERSECT/CONTIGS_{{method}}_INTERSECT",
        mapq = config["PARAMS"]["INTERSECT"]["FILTER_MAPQ_STATS"],
        nbinfeature = config["PARAMS"]["INTERSECT"]["FILTER_NB_CONTIGS_BY_FEATURE"]
    output:
        stats_contigs = f"{output_dir}14.STATS_INTERSECT/CONTIGS_{{method}}_INTERSECT/intersect_stats/nb_feature_by_gene_filter.csv",
    log:
        output = f"{output_dir}LOGS/14.STATS_INTERSECT/CONTIGS_{{method}}_INTERSECT/CONTIGS_{{method}}_INTERSECT.o",
        error = f"{output_dir}LOGS/14.STATS_INTERSECT/CONTIGS_{{method}}_INTERSECT/CONTIGS_{{method}}_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/intersect_and_contigs/STATS_CONTIGS_{{method}}_INTERSECT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
        params:
            nbinfeature : {params.nbinfeature}
        output:
            stats_contigs : {output.stats_contigs}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        mkdir -p intersect_stats 
        python3 {ikiss_obj.snakemake_scripts}/stats_intersect.py --bed {input.bed} --mapq {params.mapq} --nbinfeature {params.nbinfeature} --outdir {params.dir}/intersect_stats )1> {log.error} 2> {log.output}
        """


##################### INTERSECT KMERS

rule stats_intersect_and_outliers:
    """
    stats after intersection (bam vs gff) and outliers {using kmers detected by method}
    """
    threads: 1
    input:
        bed = rules.kmers_bedtools_intersect_method.output.bed
    params:
        dir = f"{output_dir}14.STATS_INTERSECT/OUTLIERS_{{method}}_INTERSECT",
        mapq= config["PARAMS"]["INTERSECT"]["FILTER_MAPQ_STATS"],
        nbinfeature = config["PARAMS"]["INTERSECT"]["FILTER_NB_KMERS_BY_FEATURE"]
    output:
        stats_all = f"{output_dir}14.STATS_INTERSECT/OUTLIERS_{{method}}_INTERSECT/intersect_stats/nb_feature_by_gene_filter.csv",
    log:
        output = f"{output_dir}LOGS/14.STATS_INTERSECT/OUTLIERS_{{method}}_INTERSECT/OUTLIERS_{{method}}_INTERSECT.o",
        error = f"{output_dir}LOGS/14.STATS_INTERSECT/OUTLIERS_{{method}}_INTERSECT/OUTLIERS_{{method}}_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/stats_intersect_and_outliers/STATS_OUTLIERS_{{method}}_INTERSECT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
        params:
            nbinfeature : {params.nbinfeature}
        output:
            stats_all : {output.stats_all}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        mkdir -p intersect_stats
        python3 {ikiss_obj.snakemake_scripts}/stats_intersect.py --bed {input.bed} --mapq {params.mapq} --nbinfeature {params.nbinfeature} --outdir {params.dir}/intersect_stats ) 1> {log.error} 2> {log.output}
        """


rule stats_intersect_allkmers:
    """
    stats after intersection (bam vs gff) and outliers {using kmers detected by method}
    """
    threads: 1
    input:
        bed = rules.kmers_bedtools_intersect.output.bed
    params:
        dir = f"{output_dir}14.STATS_INTERSECT/ALLKMERS_INTERSECT",
        mapq= config["PARAMS"]["INTERSECT"]["FILTER_MAPQ_STATS"],
        nbinfeature = config["PARAMS"]["INTERSECT"]["FILTER_NB_KMERS_BY_FEATURE"]
    output:
        stats_all = f"{output_dir}14.STATS_INTERSECT/ALLKMERS_INTERSECT/intersect_stats/nb_feature_by_gene_filter.csv",
    log:
        output = f"{output_dir}LOGS/14.STATS_INTERSECT/ALLKMERS_INTERSECT/ALLKMERS_INTERSECT.o",
        error = f"{output_dir}LOGS/14.STATS_INTERSECT/ALLKMERS_INTERSECT/ALLKMERS_INTERSECT.e",
    benchmark:
        f"{output_dir}BENCHMARK/stats_intersect_allkmers/STATS_ALLKMERS_INTERSECT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            bed : {input.bed}
        params :
            nbinfeature : {params.nbinfeature}
        output:
            stats_all : {output.stats_all}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (cd {params.dir}
        mkdir -p intersect_stats;
        python3 {ikiss_obj.snakemake_scripts}/stats_intersect.py --bed {input.bed} --mapq {params.mapq} --nbinfeature {params.nbinfeature} --outdir {params.dir}/intersect_stats ) 1> {log.error} 2> {log.output}
        """


################### SEX_DETECTION

rule sex_kmers:
    """
    kmers associated with the sex column of the phenotype file: their pvalues come from LFMM, their
    presence/absence pattern from the kmers table, and they are classified MALE/FEMALE/NONE by
    their frequency in each sex
    """
    threads: 1
    input:
        pvalues = f"{output_dir}7.MERGED_LFMM/merged_LFMM_pvalues.csv",
        kmers_table = rules.kmers_table.output.kmers_table,
        # the file LFMM was given, so the sex column is the one it tested
        phenotype = pheno_file if not config['PARAMS']['LFMM']['PHENOTYPE_PCA_ANALYSIS'] else rules.get_pca_from_phenotype.output.pca_variance,
    params:
        dir = f"{output_dir}17.SEX_DETECTION/",
        # filter_kmers wants the table without its .table extension
        table = f"{output_dir}2.KMERS_TABLE/kmers_table",
        column = config['PARAMS']['SEX_DETECTION']['COLUMN'],
        male = config['PARAMS']['SEX_DETECTION']['MALE'],
        female = config['PARAMS']['SEX_DETECTION']['FEMALE'],
        min_freq = config['PARAMS']['SEX_DETECTION']['MIN_FREQ'],
        max_freq = config['PARAMS']['SEX_DETECTION']['MAX_FREQ'],
        correction = config['PARAMS']['LFMM']['CORRECTION'],
        alpha = config['PARAMS']['LFMM']['ALPHA'],
    output:
        kmer_list = temp(f"{output_dir}17.SEX_DETECTION/sex_kmers.kmerlist"),
        filtered_table = temp(f"{output_dir}17.SEX_DETECTION/sex_kmers_table.txt"),
        sex_kmers = f"{output_dir}17.SEX_DETECTION/sex_kmers.tsv",
        male_fasta = f"{output_dir}17.SEX_DETECTION/male_kmers.fasta",
        female_fasta = f"{output_dir}17.SEX_DETECTION/female_kmers.fasta",
    log:
        output = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_KMERS.o",
        error = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_KMERS.e"
    benchmark:
        f"{output_dir}BENCHMARK/sex_kmers/SEX_KMERS.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            pvalues : {input.pvalues}
            phenotype : {input.phenotype}
        params:
            column : {params.column} ({params.male} = male, {params.female} = female)
        output:
            sex_kmers : {output.sex_kmers}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"],
        tools_config["ENV-MODULES"]["KMERS_GWAS"]
    shell:
        """
        (mkdir -p {params.dir}
        Rscript {ikiss_obj.snakemake_scripts}/sex_pvalues.R --pvalues {input.pvalues} --column {params.column} --correction {params.correction} --alpha {params.alpha} --output {output.kmer_list}
        if [ -s {output.kmer_list} ]; then
            filter_kmers -t {params.table} -k {output.kmer_list} -o {output.filtered_table}
        else
            # no kmer associated with the sex: filter_kmers stops on an empty list
            : > {output.filtered_table}
        fi
        python3 {ikiss_obj.snakemake_scripts}/sex_kmers.py --table {output.filtered_table} --phenotype {input.phenotype} --column {params.column} --male {params.male} --female {params.female} --min-freq {params.min_freq} --max-freq {params.max_freq} --output {output.sex_kmers} --output-male-fasta {output.male_fasta} --output-female-fasta {output.female_fasta}
        ) 1>{log.output} 2>{log.error}
        """

rule sex_contigs:
    """
    map the male and the female kmers on the contigs assembled from the LFMM outliers and count
    how many of each sex every contig carries
    """
    threads: 4
    input:
        male_fasta = rules.sex_kmers.output.male_fasta,
        female_fasta = rules.sex_kmers.output.female_fasta,
        contigs = f"{output_dir}11.ASSEMBLY_OUTLIERS_LFMM/outliers_LFMM_mergetags.fasta",
    params:
        dir = f"{output_dir}17.SEX_DETECTION/",
        # the contigs are indexed in 17.SEX_DETECTION, the index files do not go next to the
        # contigs of 11.ASSEMBLY_OUTLIERS_LFMM
        contigs = f"{output_dir}17.SEX_DETECTION/contigs_LFMM.fasta",
        options = config['PARAMS']['MAPPING_KMERS']['OPTIONS'],
        mode = config['PARAMS']['MAPPING_KMERS']['MODE'],
        flag = config['PARAMS']['MAPPING_KMERS']['FILTER_FLAG'],
        min_freq = config['PARAMS']['SEX_DETECTION']['MIN_FREQ'],
        max_freq = config['PARAMS']['SEX_DETECTION']['MAX_FREQ'],
    output:
        male_bam = f"{output_dir}17.SEX_DETECTION/male_kmers_vs_contigs.sorted.bam",
        female_bam = f"{output_dir}17.SEX_DETECTION/female_kmers_vs_contigs.sorted.bam",
        male_coverage = f"{output_dir}17.SEX_DETECTION/male_kmers_vs_contigs.coverage.txt",
        female_coverage = f"{output_dir}17.SEX_DETECTION/female_kmers_vs_contigs.coverage.txt",
        contigs_sex = f"{output_dir}17.SEX_DETECTION/contigs_sex_kmers.tsv",
    log:
        output = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_CONTIGS.o",
        error = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_CONTIGS.e"
    benchmark:
        f"{output_dir}BENCHMARK/sex_contigs/SEX_CONTIGS.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            male_fasta : {input.male_fasta}
            female_fasta : {input.female_fasta}
            contigs : {input.contigs}
        output:
            contigs_sex : {output.contigs_sex}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["BWAMEM2"],
        tools_config["ENV-MODULES"]["BWA"],
        tools_config["ENV-MODULES"]["SAMTOOLS"]
    shell:
        """
        (cd {params.dir}
        cp -f {input.contigs} {params.contigs}
        if [[ {params.mode} == bwa-mem2 ]]; then
            bwa-mem2 index {params.contigs}
        else
            bwa index {params.contigs}
        fi
        for sex in male female; do
            kmers={params.dir}${{sex}}_kmers.fasta
            bam={params.dir}${{sex}}_kmers_vs_contigs.sorted.bam
            # only the mapped kmers are counted, one alignment per kmer
            if [[ {params.mode} == bwa-mem2 ]]; then
                bwa-mem2 mem {params.options} -t {threads} {params.contigs} ${{kmers}} | samtools view -bh -F {params.flag} | samtools sort -o ${{bam}};
            else
                bwa aln {params.options} -t {threads} {params.contigs} ${{kmers}} > ${{kmers}}.sai;
                bwa samse {params.contigs} ${{kmers}}.sai ${{kmers}} | samtools view -bh -F {params.flag} | samtools sort -o ${{bam}};
                rm -f ${{kmers}}.sai;
            fi
            samtools index ${{bam}}
            samtools coverage ${{bam}} > {params.dir}${{sex}}_kmers_vs_contigs.coverage.txt
        done
        python3 {ikiss_obj.snakemake_scripts}/sex_contigs.py --male-coverage {output.male_coverage} --female-coverage {output.female_coverage} --min-freq {params.min_freq} --max-freq {params.max_freq} --output {output.contigs_sex}
        ) 1>{log.output} 2>{log.error}
        """

rule sex_intersect:
    """
    annotation of the contigs tagged MALE or FEMALE: where they align on the reference is
    intersected with the feature of the GFF, and the features each sex hits are counted
    """
    threads: 1
    input:
        contigs_sex = rules.sex_contigs.output.contigs_sex,
        bam = f"{output_dir}11.ASSEMBLY_OUTLIERS_LFMM/contigs_LFMM_vs_{basename_reference}.sorted.bam",
        gff = rules.extracting_features_from_gff.output.gff_feature,
    params:
        feature = config['PARAMS']['INTERSECT']['FEATURE'],
        # the same mapping quality as the other statistics of the intersection: a contig placed
        # nowhere in particular is not evidence that a gene carries it
        mapq = config['PARAMS']['INTERSECT']['FILTER_MAPQ_STATS'],
    output:
        bed = f"{output_dir}17.SEX_DETECTION/sex_contigs.bed",
        intersect = f"{output_dir}17.SEX_DETECTION/sex_contigs_intersect_annotation.bed",
        by_feature = f"{output_dir}17.SEX_DETECTION/sex_contigs_by_feature.tsv",
        stats = f"{output_dir}17.SEX_DETECTION/sex_intersect_stats.tsv",
    log:
        output = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_INTERSECT.o",
        error = f"{output_dir}LOGS/17.SEX_DETECTION/SEX_INTERSECT.e"
    benchmark:
        f"{output_dir}BENCHMARK/sex_intersect/SEX_INTERSECT.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            contigs_sex : {input.contigs_sex}
            bam : {input.bam}
            gff : {input.gff} ({params.feature})
        output:
            by_feature : {output.by_feature}
            stats : {output.stats}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        # like the other bedtools rules: tools_path.yaml has no env module for bedtools
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        (python3 {ikiss_obj.snakemake_scripts}/sex_contigs_bed.py --bam {input.bam} --sex {input.contigs_sex} --min-mapq {params.mapq} --output {output.bed}
        if [ -s {output.bed} ]; then
            bedtools intersect -a {output.bed} -b {input.gff} -wa -wb > {output.intersect}
        else
            # no contig tagged MALE or FEMALE in this run: bedtools has nothing to intersect
            : > {output.intersect}
        fi
        python3 {ikiss_obj.snakemake_scripts}/sex_intersect.py --bed {output.intersect} --sex {input.contigs_sex} --by-feature {output.by_feature} --stats {output.stats}
        ) 1>{log.output} 2>{log.error}
        """

################### GENOME_OFFSET

# one job per line of the climate manifest, read by module.py
CLIMATE_MODELS = ikiss_obj.climate_models

rule genome_offset:
    """
    genomic offset of the samples for one climate scenario, computed with genetic.offset (LEA) on
    the kmers LFMM associated with the phenotype
    """
    threads: 4
    input:
        outliers = f"{output_dir}7.MERGED_LFMM/merged_LFMM_outliers.csv",
        bed_dir = f"{output_dir}3.TABLE2BED/",
        reference = lambda wildcards: CLIMATE_MODELS[wildcards.model]["reference"],
        past = lambda wildcards: CLIMATE_MODELS[wildcards.model]["past"],
        future = lambda wildcards: CLIMATE_MODELS[wildcards.model]["future"],
    params:
        dir = f"{output_dir}18.GENOME_OFFSET/{{model}}/",
        k = config['PARAMS']['GENOME_OFFSET']['K'],
        pca_variance = config['PARAMS']['GENOME_OFFSET']['PCA_VARIANCE'],
        exclude_vars = ",".join(config['PARAMS']['GENOME_OFFSET']['EXCLUDE_VARS'] or []),
    output:
        # the plots of the model (climate PCA, offset, map) are written next to it
        offset = f"{output_dir}18.GENOME_OFFSET/{{model}}/genome_offset.tsv",
    log:
        output = f"{output_dir}LOGS/18.GENOME_OFFSET/{{model}}.o",
        error = f"{output_dir}LOGS/18.GENOME_OFFSET/{{model}}.e"
    benchmark:
        f"{output_dir}BENCHMARK/genome_offset/GENOME_OFFSET_{{model}}.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            outliers : {input.outliers}
            reference : {input.reference}
            future : {input.future}
        params:
            K : {params.k}
        output:
            offset : {output.offset}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    envmodules:
        tools_config["ENV-MODULES"]["R"]
    shell:
        """
        Rscript {ikiss_obj.snakemake_scripts}/genome_offset.R --outliers {input.outliers} --bed-dir {input.bed_dir} --reference {input.reference} --past {input.past} --future {input.future} --model {wildcards.model} --k {params.k} --pca-variance {params.pca_variance} --exclude-vars '{params.exclude_vars}' --outdir {params.dir} 1>{log.output} 2>{log.error}
        """

rule genome_offset_models:
    """
    the offsets of every climate scenario side by side
    """
    threads: 1
    input:
        offsets = expand(f"{output_dir}18.GENOME_OFFSET/{{model}}/genome_offset.tsv", model=list(CLIMATE_MODELS)),
    output:
        offsets = f"{output_dir}18.GENOME_OFFSET/genome_offset_all_models.tsv",
    log:
        output = f"{output_dir}LOGS/18.GENOME_OFFSET/ALL_MODELS.o",
        error = f"{output_dir}LOGS/18.GENOME_OFFSET/ALL_MODELS.e"
    benchmark:
        f"{output_dir}BENCHMARK/genome_offset_models/GENOME_OFFSET_ALL_MODELS.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            offsets : {input.offsets}
        output:
            offsets : {output.offsets}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        python3 {ikiss_obj.snakemake_scripts}/genome_offset_models.py --offsets {input.offsets} --output {output.offsets} 1>{log.output} 2>{log.error}
        """

############################ REPORT #######################################################
# report contains dico_final now

rule fastq_stats:
    """
    run fastq_stats
    """
    threads: 8
    input:
        dir = fastq_dir
    output:
        fastq_table = f"{output_dir}0.FASTQ_STATS/fastq_stats.txt"
    log:
        output = f"{output_dir}LOGS/0.FASTQ_STATS/FASTQ_STATS.o",
        error = f"{output_dir}LOGS/0.FASTQ_STATS/FASTQ_STATS.e"
    benchmark:
        f"{output_dir}BENCHMARK/fastq_stats/FASTQ_STATS.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            kmers_to_use : {input.dir}
        output:
            kmers_table: {output.fastq_table}
        log:
            output: {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        cd {input.dir}
        # one pass over the fastq files: -a adds the quality (Q20, Q30) and the GC content of
        # every file to the lengths, which is what the quality control of the report draws
        seqkit stats -a -T -j {threads} {input.dir}/*{ikiss_obj.fastq_files_ext} -o {output.fastq_table} 1>{log.error}  2> {log.output}
        sed -i 's|{input.dir}/||' {output.fastq_table}
        sed -i 's|{ikiss_obj.fastq_files_ext}||' {output.fastq_table}
        """

rule report_ikiss:
    """
    rule to recovery all files to report
    """
    threads: 1
    input:
        unpack(output_final)
    params:
        yml=f"{ikiss_obj.snakemake_scripts}/report_template/_quarto.yml",
        book_dir=f"{output_dir}REPORT/BOOK/",
        # the templates filled with the paths of the run, outside the quarto project so that
        # rendering them writes next to them and not among the sources of the book
        manipulation_dir=f"{output_dir}REPORT/MANIPULATION/",
        txt_config = f"{output_dir}/config_corrected.yaml",
        workflow_steps = ikiss_obj.tools_activated,
        methods_steps = ikiss_obj.method,
        fastq_stats = rules.fastq_stats.output.fastq_table,
        kmers_matrix_mapping_stats = rules.samtools_merge.output.stats if 'KMERS_MODULE' in ikiss_obj.tools_activated and not 'PCADAPT' in ikiss_obj.tools_activated and not 'LFMM' in ikiss_obj.tools_activated and 'MAPPING_KMERS' in ikiss_obj.tools_activated else "",
        kmers_matrix_intersect_stats = rules.merging_annotations_and_binary_matrix.output.stats_by_sample if 'KMERS_MODULE' in ikiss_obj.tools_activated and not 'PCADAPT' in ikiss_obj.tools_activated and not 'LFMM' in ikiss_obj.tools_activated and 'MAPPING_KMERS' in ikiss_obj.tools_activated and 'INTERSECT' in ikiss_obj.tools_activated else "",
        list_log_kmer_per_sample = expand(rules.kmers_gwas_per_sample.log.error, sample=SAMPLE) if 'KMERS_MODULE' in ikiss_obj.tools_activated else "",
        kmer_table_rep = f"{output_dir}3.TABLE2BED/" if 'KMERS_MODULE' in ikiss_obj.tools_activated else "",
        plots_pcadapt = f"{output_dir}6.PCADAPT" if 'PCADAPT' in ikiss_obj.tools_activated else "",
        plots_lfmm = f"{output_dir}6.LFMM" if 'LFMM' in ikiss_obj.tools_activated else "",
        plots_snmf = f"{output_dir}6.SNMF" if 'SNMF' in ikiss_obj.tools_activated else "",
        phenotype = ikiss_obj.phenotype,
        contig_size = config['PARAMS']['ASSEMBLY_KMERS']['FILTER_CONTIG_SIZE'] if 'ASSEMBLY_KMERS' in ikiss_obj.tools_activated else "",
        ref = reference_file if 'ASSEMBLY_KMERS' in ikiss_obj.tools_activated else "",
        phenotype_pca_html = rules.get_pca_from_phenotype.params.html if config['PARAMS']['LFMM']['PHENOTYPE_PCA_ANALYSIS'] else "",
        # manhattan plots of the outliers: a template filled with the paths of the run, written in
        # REPORT/MANIPULATION so the user can draw every dimension and choose the chromosomes
        manhattan_template = f"{ikiss_obj.snakemake_scripts}/report_template/manhattan_template.qmd",
        positions_by_method = expand(rules.outliers_position.output.csv, method=ikiss_obj.method) if 'MAPPING_KMERS' in ikiss_obj.tools_activated else [],
        by_chrom_by_method = expand(rules.outliers_position.output.by_chrom, method=ikiss_obj.method) if 'MAPPING_KMERS' in ikiss_obj.tools_activated else [],
        contigs_by_method = expand(rules.contigs_position.output.csv, method=ikiss_obj.method) if config['PARAMS']['ASSEMBLY_KMERS']['MANHATTAN_CONTIGS'] and 'ASSEMBLY_KMERS' in ikiss_obj.tools_activated else [],
        contigs_by_chrom_by_method = expand(rules.contigs_position.output.by_chrom, method=ikiss_obj.method) if config['PARAMS']['ASSEMBLY_KMERS']['MANHATTAN_CONTIGS'] and 'ASSEMBLY_KMERS' in ikiss_obj.tools_activated else [],
        plot_chromosomes = config['PARAMS']['MAPPING_KMERS'].get('PLOT_CHROMOSOMES', []),
        max_contigs = MAX_CONTIGS_REPORT,
        # how many p-value columns the manhattan sections of the method pages draw
        max_manhattan = MAX_MANHATTAN_DIMENSIONS,
        # annotation plots: how many genes are drawn, and the filters stats_intersect.py used
        max_genes = MAX_GENES_REPORT,
        # kmers page: how the sharing classes are named, and how the jaccard was computed
        mac = config['PARAMS']['KMERS_MODULE']['MAC'],
        soft_core = config['PARAMS']['KMERS_MODULE']['SOFT_CORE'],
        jaccard_pairs = JACCARD_PAIRS,
        jaccard_template = f"{ikiss_obj.snakemake_scripts}/report_template/jaccard_template.qmd",
        jaccard_pattern_limit = JACCARD_PATTERN_LIMIT,
        jaccard_seed = RANDOM_SEED,
        bed_dir = f"{output_dir}3.TABLE2BED/",
        # resources page: the benchmark files of every job, and how many rules its plots show
        benchmark_dir = f"{output_dir}BENCHMARK",
        max_rules = MAX_RULES_REPORT,
        intersect_mapq = config['PARAMS']['INTERSECT']['FILTER_MAPQ_STATS'] if 'INTERSECT' in ikiss_obj.tools_activated else "",
        intersect_min_kmers = config['PARAMS']['INTERSECT']['FILTER_NB_KMERS_BY_FEATURE'] if 'INTERSECT' in ikiss_obj.tools_activated else "",
        intersect_min_contigs = config['PARAMS']['INTERSECT']['FILTER_NB_CONTIGS_BY_FEATURE'] if 'INTERSECT' in ikiss_obj.tools_activated else "",
        intersect_feature = config['PARAMS']['INTERSECT']['FEATURE'] if 'INTERSECT' in ikiss_obj.tools_activated else "",
        genome_offset_dir = f"{output_dir}18.GENOME_OFFSET" if 'GENOME_OFFSET' in ikiss_obj.tools_activated else "",
        genome_offset_models = list(ikiss_obj.climate_models),
        genome_offset_k = config['PARAMS']['GENOME_OFFSET']['K'],
        # template to draw the offset plots again, and to run genome_offset.R with other settings
        genomeoffset_template = f"{ikiss_obj.snakemake_scripts}/report_template/genomeoffset_template.qmd",
        genome_offset_climate = ikiss_obj.climate_models,
        genome_offset_outliers = f"{output_dir}7.MERGED_LFMM/merged_LFMM_outliers.csv",
        genome_offset_pca_variance = config['PARAMS']['GENOME_OFFSET']['PCA_VARIANCE'],
        genome_offset_exclude = ",".join(config['PARAMS']['GENOME_OFFSET']['EXCLUDE_VARS'] or []),
        sex_column = config['PARAMS']['SEX_DETECTION']['COLUMN'],
        sex_male = config['PARAMS']['SEX_DETECTION']['MALE'],
        sex_female = config['PARAMS']['SEX_DETECTION']['FEMALE'],
        sex_min_freq = config['PARAMS']['SEX_DETECTION']['MIN_FREQ'],
        sex_max_freq = config['PARAMS']['SEX_DETECTION']['MAX_FREQ'],
        # template to go further than the snmf plots of the report
        snmf_template = f"{ikiss_obj.snakemake_scripts}/report_template/snmf_template.qmd",
        samples_file = samples_file,
        k_min = config['PARAMS']['SNMF']['K_MIN'],
        k_max = config['PARAMS']['SNMF']['K_MAX'],
        best_k = config['PARAMS']['SNMF']['BEST_K'],
    output:
        reads_ipynb = f"{output_dir}REPORT/IPYNB/reads.ipynb",
        # the kmers of every sample and the table they are merged into are one page
        kmers_info = f"{output_dir}REPORT/IPYNB/kmers_info.ipynb",
        # time and peak memory of every rule, from the benchmark files
        resources_ipynb = f"{output_dir}REPORT/IPYNB/resources.ipynb",
        yml=f"{output_dir}REPORT/BOOK/_quarto.yml"
    log:
        output = f"{output_dir}LOGS/REPORT/report.o",
        error = f"{output_dir}LOGS/REPORT/report.e",
    benchmark:
        f"{output_dir}BENCHMARK/report_ikiss/Report.txt",
    message:
        """
        Launching {rule}
        threads: {threads}
        input:
            {input}
        output:
            reads_ipynb : {output.reads_ipynb}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    script:
        f"{ikiss_obj.snakemake_scripts}/report.py"

#### QUARTO REPORT

localrules: rule_graph

rule rule_graph:
    """
    run dag on {rule}
    """
    threads: 1
    input:
        #quarto_conf = f"{output_dir}REPORT/BOOK/_quarto.yml",
        config = f"{output_dir}config_corrected.yaml"
    params:
        cmd = ikiss_obj.string_to_dag
    output:
        # .dot rendered by quarto itself (```{{dot}}``` block): no dot binary needed, and the
        # graph ends up as a vector image in the html report
        dag = f"{output_dir}REPORT/BOOK/dag.dot",
    log:
        output = f"{output_dir}REPORT/LOGS/GRAPH.o",
        error = f"{output_dir}REPORT/LOGS/GRAPH.e"
    message:
        """
        making dag ...
        {params.cmd} --configfile {input.config} > {output.dag}
        """
    shell:
        """
        # the command writes the graph on its stdout, but a warning of snakemake or of an executor
        # plugin lands there too on some clusters. It would stay in the middle of the dot file and
        # quarto stops the whole report on `syntax error in file dag.dot`: only what follows the
        # `digraph` line is kept, and an output holding no graph at all fails here, with its log
        ({params.cmd} --configfile {input.config} | sed -n '/^digraph/,$p' 1>{output.dag}) 2>{log.error}
        if [ ! -s {output.dag} ]; then
            echo "no digraph in the output of: {params.cmd} --configfile {input.config}" >>{log.error}
            exit 1
        fi
        """

rule run_report_snakemake:
    """
    run dag on {rule}
    """
    threads: 1
    input:
        #quarto_conf = f"{output_dir}REPORT/BOOK/_quarto.yml",
        config = f"{output_dir}config_corrected.yaml"
    output:
        report_snakemake = f"{output_dir}REPORT/ikiss_report/snake_report.html"
    params:
        cmd = ikiss_obj.string_to_dag
    log:
        output = f"{output_dir}REPORT/LOGS/REPORT-SNAKE.o",
        error = f"{output_dir}REPORT/LOGS/REPORT-SNAKE.e"
    message:
        """
        making report with dag and config file...
        """
    shell:
        """
        cd {ikiss_obj.install_path}/
        {params.cmd} --configfile {input.config} --report {output.report_snakemake} 1>{log.output} 2>{log.error}
        """

rule report_about_workflow:
    """from culebront project: build qmd for tools version dag and config"""
    threads: 1
    input:
        dag = rules.rule_graph.output.dag,
        #report_snakemake= rules.run_report_snakemake.output.report_snakemake
    output:
        about_ipynb = f"{output_dir}REPORT/IPYNB/about_workflow.ipynb"
    params:
        config_yaml = ikiss_obj.export_use_yaml,
        versions = f"{output_dir}versions.csv",
    log:
        output=f"{output_dir}REPORT/LOGS/report_about_workflow_QMD.o",
        error=f"{output_dir}REPORT/LOGS/report_about_workflow_QMD.e"
    message:
        """
        Launching {rule}
        threads : {threads}
        input:
             dag : {input.dag}
             versions : {params.versions}
        log:
            output : {log.output}
            error: {log.error}
        """
    script:
        f"{ikiss_obj.snakemake_scripts}/generate_tools.py"

rule ipynb_convert_qmd:
    """build qmd for tools version dag and config"""
    threads: 1
    input:
        about_ipynb= rules.report_about_workflow.output.about_ipynb,
        reads_ipynb=rules.report_ikiss.output.reads_ipynb,
        kmers_info=rules.report_ikiss.output.kmers_info,
        resources_ipynb=rules.report_ikiss.output.resources_ipynb,
    output:
        reads_qmd= f"{output_dir}REPORT/BOOK/reads.qmd",
        kmers_qmd= f"{output_dir}REPORT/BOOK/kmers_info.qmd",
        resources_qmd= f"{output_dir}REPORT/BOOK/resources.qmd",
        about_qmd= f"{output_dir}REPORT/BOOK/about_workflow.qmd",
    params:
        ipynb_dir = f"{output_dir}REPORT/IPYNB/",
        book_dir = f"{output_dir}REPORT/BOOK/",
    log:
        output=f"{output_dir}REPORT/LOGS/ipynb_convert_qmd.o",
        error=f"{output_dir}REPORT/LOGS/ipynb_convert_qmd.e"
    message:
        """
        Launching {rule}
        threads : {threads}
        input:
             ipynb : {params.ipynb_dir}
        output:
             qmd : {output.about_qmd}
        log:
            output : {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        # every notebook report.py wrote, whatever the activated modules: the qmd goes straight
        # into the book, no copy of the ipynb is left there
        (mkdir -p {params.book_dir}
        for notebook in {params.ipynb_dir}*.ipynb; do
            quarto convert "$notebook" --output {params.book_dir}$(basename "${{notebook%.ipynb}}").qmd
        done
        )1>{log.output} 2>{log.error}
        """
rule build_book:
    """copy template files to build report"""
    threads: 1
    input:
        qmd_about_wf = rules.ipynb_convert_qmd.output.about_qmd,
        qmd_reads = rules.ipynb_convert_qmd.output.reads_qmd,
        qmd_kmers= rules.ipynb_convert_qmd.output.kmers_qmd,
        qmd_resources= rules.ipynb_convert_qmd.output.resources_qmd,
        yml=rules.report_ikiss.output.yml,
        index_qmd=f"{ikiss_obj.snakemake_scripts}/report_template/index.qmd",
        ref_bib=f"{ikiss_obj.snakemake_scripts}/report_template/references.bib",
        ref_qmd=f"{ikiss_obj.snakemake_scripts}/report_template/references.qmd",
        csl=f"{ikiss_obj.snakemake_scripts}/report_template/cellalina.csl",
        logo=f"{ikiss_obj.install_path}/logo_ikiss.png",
        # the logos of the institutes, shown at the bottom of the home page
        logos=[f"{ikiss_obj.install_path}/img/ird_small.png",
               f"{ikiss_obj.install_path}/img/Logo-DIADE-green1.svg"],
    output:
        html = f"{output_dir}REPORT/BOOK/ikiss_report/index.html"
    params:
        book_dir=directory(f"{output_dir}REPORT/BOOK/"),
        input_dir=f"{ikiss_obj.snakemake_scripts}/",
        link_html=f"{output_dir}REPORT/ikiss_report.html"
    log:
        output=f"{output_dir}REPORT/LOGS/build_book.o",
        error=f"{output_dir}REPORT/LOGS/build_book.e"
    message:
        """
        Launching {rule}
        threads : {threads}
        input:
             qmd : {input.yml}
        output:
             about_qmd : {output.html}
        log:
            output : {log.output}
            error: {log.error}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
            (mkdir -p {params.book_dir};
            cp {input.index_qmd} {params.book_dir};
            cp {input.ref_bib} {params.book_dir};
            cp {input.csl} {params.book_dir};
            cp {input.logo} {input.logos} {params.book_dir};
            cp {input.ref_qmd} {params.book_dir};
            cd {params.book_dir};
            #quarto add --no-prompt --profile quarto quarto-ext/lightbox;
            quarto render --to html
            # the target is written relative to REPORT/, so the link survives a move of the
            # output directory (an absolute one points at the machine it was built on)
            ln -s -f BOOK/ikiss_report/index.html {params.link_html}
            ) 1>{log.output} 2>{log.error}
        """

rule render_manipulation:
    """
    render the pages of REPORT/MANIPULATION, running their cells. Nothing of the report depends on
    it: ask for it by name, `ikiss run -c config.yaml --cores 4 render_manipulation`
    """
    threads: 1
    input:
        # report.py filled the templates with the paths of the run
        yml = rules.report_ikiss.output.yml,
    params:
        dir = f"{output_dir}REPORT/MANIPULATION/",
        pages = manipulation_pages(),
    output:
        html = expand(f"{output_dir}REPORT/MANIPULATION/{{page}}.html", page=manipulation_pages()),
        # a run with no page to render still has something to ask for
        ok = touch(f"{output_dir}REPORT/MANIPULATION/rendered.ok"),
    log:
        output = f"{output_dir}REPORT/LOGS/render_manipulation.o",
        error = f"{output_dir}REPORT/LOGS/render_manipulation.e"
    benchmark:
        f"{output_dir}BENCHMARK/render_manipulation/RENDER_MANIPULATION.txt"
    message:
        """
        Launching {rule}
        threads: {threads}
        params:
            pages : {params.pages}
        output:
            html : {output.html}
        """
    container:
        tools_config['APPTAINER']['TOOLS']
    shell:
        """
        # XDG_RUNTIME_DIR: without a writable one quarto stops on `mkdir /run/user/<uid>/jt`.
        # --no-execute-daemon: a kernel kept alive between renders is reused with the environment
        # it was started in, and gives a ModuleNotFoundError on a module that is installed
        (cd {params.dir}
        export XDG_RUNTIME_DIR="$PWD/.runtime"
        mkdir -p "$XDG_RUNTIME_DIR"
        for page in {params.pages}; do
            quarto render $page.qmd --to html --no-execute-daemon
        done
        ) 1>{log.output} 2>{log.error}
        """

