"""
All rules/logic related to converting filetypes from uncompressed to compressed
(big as UCSC calls it). E.g. a from .peak -> .bigpeak
"""


rule bedgraphish_to_bedgraph:
    """
    Convert the bedgraph-ish file generated by genrich into a true bedgraph file.
    """
    input:
        expand("{result_dir}/genrich/{{assembly}}-{{sample}}.bdgish", **config),
    output:
        bedgraph=temp(expand("{result_dir}/genrich/{{assembly}}-{{sample}}.bedgraph", **config)),
    log:
        expand("{log_dir}/bedgraphish_to_bedgraph/{{assembly}}-{{sample}}.log", **config),
    benchmark:
        expand("{benchmark_dir}/bedgraphish_to_bedgraph/{{assembly}}-{{sample}}.log", **config)[0]
    shell:
        """
        splits=$(grep -Pno "([^\/]*)(?=\.bam)" {input})
        splits=($splits)
        lncnt=$(wc -l {input} | echo "$(grep -Po ".*(?=\ )")":)
        splits+=($lncnt)

        counter=1
        for split in "${{splits[@]::${{#splits[@]}}-1}}";
        do
            filename=$(grep -Po "(?<=:).*" <<< $split);
            if [[ $filename =~ {wildcards.sample} ]]; then
                startnr=$(grep -Po ".*(?=\:)" <<< $split);
                endnr=$(grep -Po ".*(?=\:)" <<< ${{splits[counter]}});

                lines="NR>=startnr && NR<=endnr {{ print \$1, \$2, \$3, \$4 }}"
                lines=${{lines/startnr/$((startnr + 2))}}
                lines=${{lines/endnr/$((endnr - 1))}}

                awk "$lines" {input} > {output}
            fi
            ((counter++))
        done
        """


def find_bedgraph(wildcards):
    if wildcards.peak_caller == "genrich":
        suffix = ".bedgraph"
    elif wildcards.peak_caller == "macs2":
        suffix = "_treat_pileup.bdg"
    else:
        suffix = ".bedgraph"

    return f"{config['result_dir']}/{wildcards.peak_caller}/{wildcards.assembly}-{wildcards.sample}{suffix}"


rule bedgraph_bigwig:
    """
    Convert a bedgraph file into a bigwig.
    """
    input:
        bedgraph=find_bedgraph,
        sizes=rules.get_genome_support_files.output.sizes,
    output:
        out=expand("{result_dir}/{{peak_caller}}/{{assembly}}-{{sample}}.bw", **config),
        tmp=temp(expand("{result_dir}/{{peak_caller}}/{{assembly}}-{{sample}}.bedgraphtmp", **config)),
    log:
        expand("{log_dir}/bedgraph_bigwig/{{peak_caller}}/{{assembly}}-{{sample}}.log", **config),
    benchmark:
        expand("{benchmark_dir}/bedgraphish_to_bedgraph/{{assembly}}-{{sample}}-{{peak_caller}}.log", **config)[0]
    conda:
        "../envs/ucsc.yaml"
    shell:
        """
        awk -v OFS='\\t' '{{print $1, $2, $3, $4}}' {input.bedgraph} | sed '/experimental/d' |
        bedSort /dev/stdin {output.tmp} > {log} 2>&1
        bedGraphToBigWig {output.tmp} {input.sizes} {output.out} >> {log} 2>&1
        """


def get_bigpeak_columns(wildcards):
    if get_peak_ftype(wildcards.peak_caller) == "narrowPeak":
        return 10
    if get_peak_ftype(wildcards.peak_caller) == "broadPeak":
        if len(TREPS_FROM_BREP[(wildcards.sample, wildcards.assembly)]) == 1:
            return 9
        return 12
    logger.error(
        f"Something went wrong inferring the correct column for bigpeak for peak type {get_peak_ftype(wildcards.peak_caller)}. "
        "Please make an issue on github if this is unexpected behaviour!"
    )
    os._exit(1)  # noqa


def get_bigpeak_type(wildcards):
    if get_peak_ftype(wildcards.peak_caller) == "narrowPeak":
        return "bed6+4"
    if get_peak_ftype(wildcards.peak_caller) == "broadPeak":
        if len(TREPS_FROM_BREP[(wildcards.sample, wildcards.assembly)]) == 1:
            return "bed6+3"
        return "bed12"
    logger.error(
        f"Something went wrong inferring the correct bed format for bigpeaks for peak type {get_peak_ftype(wildcards.peak_caller)}. "
        "Please make an issue on github if this is unexpected behaviour!"
    )


def get_bigpeak_schema(wildcards):
    if get_peak_ftype(wildcards.peak_caller) == "narrowPeak":
        return f"{config['rule_dir']}/../schemas/bignarrowPeak.as"
    if get_peak_ftype(wildcards.peak_caller) == "broadPeak":
        if len(TREPS_FROM_BREP[(wildcards.sample, wildcards.assembly)]) == 1:
            return f"{config['rule_dir']}/../schemas/bigbroadPeak.as"
        return f"{config['rule_dir']}/../schemas/bigBed.as"
    logger.error(
        f"Something went wrong inferring the correct schema for bigpeaks for peak type {get_peak_ftype(wildcards.peak_caller)}. "
        "Please make an issue on github if this is unexpected behaviour!"
    )


rule peak_bigpeak:
    """
    Convert a narrowpeak file into a bignarrowpeak file.
    https://genome-source.gi.ucsc.edu/gitlist/kent.git/tree/master/src/hg/lib/
    """
    input:
        narrowpeak=expand("{result_dir}/{{peak_caller}}/{{assembly}}-{{sample}}_peaks.{{peak}}", **config),
        sizes=rules.get_genome_support_files.output.sizes,
    output:
        out=expand("{result_dir}/{{peak_caller}}/{{assembly}}-{{sample}}.big{{peak}}", **config),
        tmp=temp(expand("{result_dir}/{{peak_caller}}/{{assembly}}-{{sample}}.tmp.{{peak}}", **config)),
    log:
        expand("{log_dir}/narrowpeak_bignarrowpeak/{{peak_caller}}/{{assembly}}-{{sample}}-{{peak}}.log", **config),
    benchmark:
        expand("{benchmark_dir}/bedgraphish_to_bedgraph/{{assembly}}-{{sample}}-{{peak_caller}}-{{peak}}.log", **config)[
            0
        ]
    conda:
        "../envs/ucsc.yaml"
    params:
        schema=lambda wildcards: get_bigpeak_schema(wildcards),
        type=lambda wildcards: get_bigpeak_type(wildcards),
        columns=lambda wildcards: get_bigpeak_columns(wildcards),
    shell:
        """
        # keep first 10 columns, idr adds extra columns we do not need for our bigpeak
        cut -d$'\t' -f 1-{params.columns} {input.narrowpeak} |
        awk -v OFS="\t" '{{$5=$5>1000?1000:$5}} {{print}}' | 
        bedSort /dev/stdin {output.tmp} > {log} 2>&1

        bedToBigBed -type={params.type} -as={params.schema} {output.tmp} \
        {input.sizes} {output.out} >> {log} 2>&1
        """


if WORKFLOW != "rna_seq":
    rule bam_bigwig:
        """
        Convert a bam file into a bigwig file.
        """
        input:
            bam=expand("{final_bam_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bam", **config),
            bai=expand("{final_bam_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bam.bai", **config),
        output:
            expand("{bigwig_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bw", **config),
        params:
            flags=config["deeptools_bamcoverage"],
        wildcard_constraints:
            sorting=config["bam_sort_order"] if config.get("bam_sort_order") else "",
        log:
            expand("{log_dir}/bam_bigwig/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.log", **config),
        benchmark:
            expand(
                "{benchmark_dir}/bam_bigwig/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.benchmark.txt", **config
            )[0]
        conda:
            "../envs/deeptools.yaml"
        threads: 16
        resources:
            deeptools_limit=lambda wildcards, threads: threads,
        shell:
            """
            bamCoverage --bam {input.bam} --outFileName {output} --numberOfProcessors {threads} {params.flags} --verbose >> {log} 2>&1
            """

else:
    rule bam_bigwig:
        """
        Convert a bam file into a bigwig file (or two).
        Can output strand specific bams. The reverse bam is created "secretly" to avoid checkpoints.
        """
        input:
            bam=expand("{final_bam_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bam",**config),
            bai=expand("{final_bam_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bam.bai",**config),
            report=rules.infer_strandedness.output,
        output:
            expand("{bigwig_dir}/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.bw",**config),
        params:
            strandedness=lambda wildcards, input: get_strandedness(input.report[0]),
            flags=config["deeptools_bamcoverage"],
        wildcard_constraints:
            sorting=config["bam_sort_order"] if config.get("bam_sort_order") else "",
        log:
            expand("{log_dir}/bam_bigwig/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.log",**config),
        benchmark:
            expand("{benchmark_dir}/bam_bigwig/{{assembly}}-{{sample}}.{{sorter}}-{{sorting}}.benchmark.txt",**config)[0]
        conda:
            "../envs/deeptools.yaml"
        threads: 16
        resources:
            deeptools_limit=lambda wildcards, threads: threads,
        shell:
            """
            if [[ "{params.strandedness}" == "no" ]]; then
              echo "generating unstranded bam" > {log}
              bamCoverage --bam {input.bam} --outFileName {output} --numberOfProcessors {threads} {params.flags} --verbose >> {log} 2>&1
            else
              echo "generating stranded bams" > {log}
			  out={output}
              rev=$(echo "${{out/{wildcards.sorting}.bw/{wildcards.sorting}.rev.bw}}")
              rm -f $rev  # in case of a rerun
              bamCoverage --bam {input.bam} --outFileName $rev     --filterRNAstrand reverse --numberOfProcessors {threads} {params.flags} --verbose >> {log} 2>&1
              bamCoverage --bam {input.bam} --outFileName {output} --filterRNAstrand forward --numberOfProcessors {threads} {params.flags} --verbose >> {log} 2>&1
            fi
            """
