shell:
    """
    set -euo pipefail

    bins=({params.bins})
    outfiles=({params.outfile_no_gz})

    for i in "${{!bins[@]}}"; do 
        bash workflow/scripts/find_binned_average.sh \
             {input.reference}.fai \
             ${{bins[$i]}} \
             {params.covFilter} \
             {input.bed} \
             ${{outfiles[$i]}} # find_binned_average gzips the file, so different than {{output.outfile}}
    done

    """   

shell:
    """
    if [ {params.nome} == "True" ]; then
        biscuit mergecg -N {input.ref} {input.bed} 1> {params.mergecg} 2> {log.mergecg}
    else
        biscuit mergecg {input.ref} {input.bed} 1> {params.mergecg} 2> {log.mergecg}
    fi

    bgzip {params.mergecg} 2> {log.mergecg_gz}
    tabix -p bed {output.mergecg_gz} 2> {log.mergecg_tbi}
    """

shell:
    """
    biscuit vcf2bed -t snp {input.vcf_gz} > {params.snp_bed} 2> {log}
    bgzip {params.snp_bed}
    tabix -p bed {output.snp_bed_gz}
    """

shell:
    """
    if [[ "$(zcat {input.snps} | head -n 1 | wc -l)" == "1" ]]; then
        if [ {params.nome} == "True" ]; then
            biscuit epiread -N -@ {threads} -B <(zcat {input.snps}) {input.ref} {input.bam} | sort -k1,1 -k2,2n > {params.epibed} 2> {log}
        else
            biscuit epiread -@ {threads} -B <(zcat {input.snps}) {input.ref} {input.bam} | sort -k1,1 -k2,2n > {params.epibed} 2> {log}
        fi
    else
        if [ {params.nome} == "True" ]; then
            biscuit epiread -N {input.ref} {input.bam} | sort -k1,1 -k2,2n > {params.epibed} 2> {log}
        else
            biscuit epiread {input.ref} {input.bam} | sort -k1,1 -k2,2n > {params.epibed} 2> {log}
        fi
    fi

    bgzip {params.epibed}
    tabix -p bed {output.epibed_gz}
    """

script:
    '../scripts/plot_covg_uniformity.py'

shell:
    """
    set +o pipefail
    mkdir -p {params.outdir}

    # Unfiltered coverage
    bedtools genomecov -bga -split -ibam {input.bam} | \
    LC_ALL=C sort --parallel={threads} -k1,1 -k2,2n -T {params.outdir} | \
    bedtools intersect -a {params.cpg} -b - -wo -sorted | \
    bedtools groupby -g 1-3 -c 7 -o min | \
    gzip -c > {output.unf}

    # MAPQ >= 40 filtered coverage
    samtools view -q 40 -b {input.bam} | \
    bedtools genomecov -bga -split -ibam stdin | \
    LC_ALL=C sort --parallel={threads} -k1,1 -k2,2n -T {params.outdir} | \
    bedtools intersect -a {params.cpg} -b - -wo -sorted | \
    bedtools groupby -g 1-3 -c 7 -o min | \
    gzip -c > {output.fil}
    """

shell:
    """
    echo -e "Feature\tCpG_Count\tQ40_Reads\tAll_Reads" > {output}

    # All CpGs
    zcat {input.unf} | wc -l | awk '{{ printf("AllCpGs\\t%s", $1) }}' >> {output}
    zcat {input.fil} | awk '$4>0{{ a += 1 }} END{{ printf("\\t%d", a) }}' >> {output}
    zcat {input.unf} | awk '$4>0{{ a += 1 }} END{{ printf("\\t%d\\n", a) }}' >> {output}

    # CpG Island CpGs
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.cgi}) -sorted | wc -l | \
    awk '{{ printf("CGICpGs\\t%s", $1) }}' >> {output}
    bedtools intersect -a {input.fil} -b <(bedtools merge -i {params.cgi}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d", a) }}' >> {output}
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.cgi}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d\\n", a) }}' >> {output}

    # Exonic CpGs
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.exn}) -sorted | wc -l | \
    awk '{{ printf("ExonicCpGs\\t%s", $1) }}' >> {output}
    bedtools intersect -a {input.fil} -b <(bedtools merge -i {params.exn}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d", a) }}' >> {output}
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.exn}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d\\n", a) }}' >> {output}

    # Genic CpGs
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.gen}) -sorted | wc -l | \
    awk '{{ printf("GenicCpGs\\t%s", $1) }}' >> {output}
    bedtools intersect -a {input.fil} -b <(bedtools merge -i {params.gen}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d", a) }}' >> {output}
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.gen}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d\\n", a) }}' >> {output}

    # Repeat-masked CpGs
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.msk}) -sorted | wc -l | \
    awk '{{ printf("RepeatCpGs\\t%s", $1) }}' >> {output}
    bedtools intersect -a {input.fil} -b <(bedtools merge -i {params.msk}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d", a) }}' >> {output}
    bedtools intersect -a {input.unf} -b <(bedtools merge -i {params.msk}) -sorted | \
    awk '$4>0{{ a += 1 }} END{{ printf("\\t%d\\n", a) }}' >> {output}
    """

shell:
    """
    # Total number of CpG islands
    zcat {params.cgi} | wc -l | awk '{{ printf("n_cpg_islands\\t%s\\n", $1) }}' > {output}

    # Number of CpG islands with at least one read with MAPQ>=40 covering 1+, 3+, 5+, and 7+ CpGs in that CpG island
    bedtools intersect -a {input.fil} -b <(bedtools merge -i {params.cgi}) -sorted -wo | \
    awk '$4>0 {{ print $5":"$6"-"$7 }}' | uniq -c | \
    awk '{{ if ($1 >= 1) {{ a += 1; }} if ($1 >= 3) {{ b += 1; }} if ($1 >= 5) {{ c += 1; }} if ($1 >= 7) {{ d += 1; }} }} END {{ printf("%d %d %d %d", a, b, c, d) }}' | \
    awk '{{ printf("one_cpg\\t%s\\nthree_cpgs\\t%s\\nfive_cpgs\\t%s\\nseven_cpgs\\t%s\\n", $1, $2, $3, $4) }}' >> {output}
    """

script:
    '../scripts/plot_cpg_stats_features.py'

script:
    '../scripts/plot_cpg_stats_cgi.py'

shell:
    """
    set +o pipefail;

    echo "Create output directory" 1> {log}
    mkdir -p {params.dir}

    echo "Set reference location" 1> {log}
    # Reference location
    REFLOC={input.ref}.fai

    echo "Retrieve CpG island file" 1> {log}
    # CpG Islands from UCSC on only the canonical chromosomes
    wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/{params.gen}/database/cpgIslandExt.txt.gz | \
    gunzip -c | \
    awk 'BEGIN{{ OFS="\\t"; }}{{ print $2, $3, $4; }}' | \
    awk '{{ if ($1 ~ /^chr[1234567890XYM]{{1,2}}$/) {{ print }} }}' | \
    bedtools sort -i - | \
    gzip -c > {output.cgi}
    echo "CpG island file retrieved" 1> {log}

    # Create middles for finding locations
    bedtools slop \
        -i {output.cgi} \
        -g ${{REFLOC}} \
        -b 2000 | \
    bedtools merge -i - \
    > {params.tmp_sho}

    bedtools slop \
        -i {output.cgi} \
        -g ${{REFLOC}} \
        -b 4000 | \
    bedtools merge -i - \
    > {params.tmp_shl}

    # CpG Open Seas (intervening locations)
    sort -k1,1 -k2,2n ${{REFLOC}} | \
    bedtools complement -L -i {params.tmp_shl} -g - | \
    gzip -c > {output.opn}

    # CpG Shelves (shores +/- 2kb)
    bedtools subtract \
        -a {params.tmp_shl} \
        -b {params.tmp_sho} | \
    gzip -c > {output.shl}

    # CpG Shores (island +/- 2kb)
    bedtools subtract \
        -a {params.tmp_sho} \
        -b {output.cgi} | \
    gzip -c > {output.sho}

    # Clean up CpG temporary files
    rm -f {params.tmp_shl} {params.tmp_sho}
    echo "CpG seascape created and cleaned up" 1> {log}

    # BISCUIT assets (includes CpG and top/bottom 10% GC-content files)
    build_biscuit_QC_assets.pl --outdir {params.dir} --ref {input.ref} 1> {log}

    # Repeat-masked file
    if [[ {params.gen} != "mm9" ]]; then
        wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/{params.gen}/database/rmsk.txt.gz | \
        gunzip -c | \
        awk 'BEGIN{{ OFS="\\t"; }}{{ print $6, $7, $8; }}' | \
        awk '{{ if ($1 ~ /^chr[1234567890XYM]{{1,2}}$/) {{ print }} }}' | \
        bedtools sort -i - | \
        bedtools merge -i - | \
        gzip -c > {output.msk}
    else
        wget -q -nd -r -nH -np --accept-regex 'chr[1234567890XYM]*_rmsk.txt.gz$' \
            --regex-type=posix http://hgdownload.cse.ucsc.edu/goldenpath/mm9/database/

        cat chr*_rmsk.txt.gz | \
        gunzip -c | \
        awk 'BEGIN{{ OFS="\\t"; }}{{ print $6, $7, $8; }}' | \
        bedtools sort -i - | \
        bedtools merge -i - | \
        gzip -c > {output.msk}

        if [[ -f index.html ]]; then rm index.html; fi
        if [[ -f robots.txt ]]; then rm robots.txt; fi
        rm chr*_rmsk.txt.gz
    fi

    # Genic, intergenic, and exon files
    wget -qO- {params.annt} | \
    gunzip -c | \
    awk 'BEGIN{{ OFS="\\t"; }}{{ if ($3 == "gene") {{ print $1, $4, $5 > "{params.tmp_gen}"; }} else if ($3 == "exon") {{ print $1, $4, $5 > "{params.tmp_exn}"; }} }}'

    bedtools sort -i {params.tmp_gen} | \
    bedtools merge -i - | \
    gzip -c > {output.gen}

    bedtools sort -i {params.tmp_exn} | \
    bedtools merge -i - | \
    gzip -c > {output.exn}

    sort -k1,1 -k2,2n {input.ref}.fai | \
    bedtools complement -L -i {output.gen} -g - | \
    gzip -c > {output.inr}

    # Clean up
    rm -f {params.tmp_gen} {params.tmp_exn}
    """       

shell:
    """
    set +o pipefail;

    #cd {params.dir}

    wget -P {params.dir} --no-check-certificate -q https://bismap.hoffmanlab.org/raw/{params.gen}/k100.bismap.bedgraph.gz

    # Create average mappability scores in 10kb windows
    zcat {output.bis} | sort -k1,1 -k2,2n | gzip > {output.sor}

    bedtools makewindows -w 10000 -g {input.ref}.fai | gzip 1> {output.byn} 2> {log}
    bedtools map -a {output.byn} -b {output.sor} -c 4 -o mean | gzip 1> {output.wgt} 2> {log}

    bedtools intersect -a {output.bis} -b {input.cpg} | gzip -c 1> {output.cpg} 2> {log}
    bedtools intersect -a {output.bis} -b {input.cgi} | gzip -c 1> {output.cgi} 2> {log}
    bedtools intersect -a {output.bis} -b {input.exn} | gzip -c 1> {output.exn} 2> {log}
    bedtools intersect -a {output.bis} -b {input.gen} | gzip -c 1> {output.gen} 2> {log}
    bedtools intersect -a {output.bis} -b {input.inr} | gzip -c 1> {output.inr} 2> {log}
    bedtools intersect -a {output.bis} -b {input.msk} | gzip -c 1> {output.msk} 2> {log}
    """

shell:
    """
    bash workflow/scripts/create_context_beds.sh \
        -t {threads} \
        -o {params.dir} \
        {input.ref} \
        {params.cpg}/{input.cpg} 2> {log}
    """

shell:
    """
    set +o pipefail
    mkdir -p {params.outdir}

    # Find genomic coverage
    samtools view -hb -F 0x4 -q 40 {input.bam} | \
    bedtools genomecov -bg -ibam stdin | \
    gzip -c > {output.cov}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.cpg} -b {input.cov} -wo | gzip -c > {output.cpg}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.cgi} -b {input.cov} -wo | gzip -c > {output.cgi}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.intergenic} -b {input.cov} -wo | gzip -c > {output.intergenic}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.exon} -b {input.cov} -wo | gzip -c > {output.exon}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.genic} -b {input.cov} -wo | gzip -c > {output.genic}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.rmsk} -b {input.cov} -wo | gzip -c > {output.rmsk}
    """

shell:
    """
    set +o pipefail

    # Find intersections
    bedtools intersect -a {params.mapped} -b {input.cov} -wo | gzip -c > {output.mapped}
    """

script:
    '../scripts/calculate_feature_sizes.py'

script:
    '../scripts/plot_obs_exp.py'

script:
    '../scripts/rename.R'

shell:
    """
    fastq_screen --bisulfite --conf {params.conf} --outdir {params.output_dir} {input} 2> {log.fastq_screen}
    """

shell:
    """
    set +o pipefail;
    QC.sh \
        -o {params.output_dir} \
        --vcf {input.vcf} \
        {params.assets} \
        {input.ref} \
        {wildcards.sample} \
        {input.bam} \
        2> {log}
    """

script:
    '../scripts/plot_percent_covered.py'

shell:
    """
    mkdir -p {params.lambda_dir}
    mkdir -p {params.puc19_dir}

    # >J02459.1 Escherichia phage Lambda, complete genome - UNMETHYLATED CONTROL
    zcat {input.bed} | {{ grep '^J02459.1' || true; }} > {output.lambda_bed} 2> {log.lambda_log}

    # >M77789.2 Cloning vector pUC19, complete sequence - METHYLATED CONTROL
    zcat {input.bed} | {{ grep '^M77789.2' || true; }}  > {output.puc19_bed} 2> {log.puc19_log}
    """

script:
    '../scripts/control_vector.R'

shell:
    """
    python3 workflow/scripts/region_centered_bin_averages.py {params.args_list} {input.region_bed} | \
    sort -k1,1 -k2,2n > {output.region_centered_bins}
    """

shell:
    """
    # do the intersecting for this sample
    # column 14 has the methylation
    bedtools intersect \
        -a {input.region_bins} \
        -b {input.merged_sample_bed} \
        -sorted -wo | \
    bedtools groupby \
        -g 1-10 \
        -c 14 \
        -o mean | \
    awk '{{print $7,"\\t",$8,"\\t",$9,"\\t",$4,"\\t",$10,"\\t",$11}}' > {output.bed}
    """

set -euo pipefail

# Check for bedtools, awk, and parallel in PATH
function check_path {
    if [[ `which bedtools 2>&1 > /dev/null` ]]; then
        >&2 echo "bedtools does not exist in PATH"
        exit 1
    else
        >&2 echo "Using bedtools found at: `which bedtools`"
    fi
    if [[ `which awk 2>&1 > /dev/null` ]]; then
        >&2 echo "awk does not exist in PATH"
        exit 1
    else
        if awk --version | grep -q GNU; then
            >&2 echo "Using GNU awk found at: `which awk`"
        else
            >&2 echo "It doesn't appear you are using GNU awk"
            >&2 echo "Try adding GNU awk at the front of PATH"
            exit 1
        fi
    fi
    if [[ `which parallel 2>&1 > /dev/null` ]]; then
        >&2 echo "parallel does not exist in PATH"
        >&2 echo "Make sure to add GNU parallel to PATH"
        exit 1
    else
        if parallel --version | grep -q GNU; then
            >&2 echo "Using GNU parallel found at: `which parallel`"
        else
            >&2 echo "It doesn't appear you are using GNU parallel."
            >&2 echo "Try adding GNU parallel at the front of PATH"
            exit 1
        fi
    fi
}

expand_windows() {
    # Load chromosome sizes to help with keeping extracted ranges within bounds of the chromosomes
    chr_sizes=""
    while read line; do
        chr="$(echo ${line} | awk '{ print $1 }')"
        siz="$(echo ${line} | awk '{ print $2 }')"

        if [ -z "$chr_sizes" ]; then
            chr_sizes="${chr},${siz}"
        else
            chr_sizes="${chr_sizes};${chr},${siz}"
        fi
    done < $1

    awk -v sizes="${chr_sizes}" \
    'BEGIN { 
        OFS="\t"

        # Create a key-value dictionary of chromosome sizes
        split(sizes, tmp1, ";")
        for (i in tmp1) {
            split(tmp1[i], tmp2, ",")
            chr_lengths[tmp2[1]] = tmp2[2]
        }
    } {
        chr = $1
        if ($2 - 35 < 0) { beg = 0 }
        else { beg = $2 - 35 }
        if ($3 + 35 > chr_lengths[chr]) { end = chr_lengths[chr] }
        else { end = $3 + 35 }

        print chr, beg, end
    }' $2
}
export -f expand_windows

find_context() {
    awk \
    'BEGIN { OFS = "\t" } {
        # Extract position
        match($1, /([^:]*):([0-9]*)-([0-9]*)/, position)
        chr = position[1]
        beg = position[2] + 35
        end = position[3] - 35

        # Reference sequence
        a1=toupper($2)
        a2=a1
        a3=a1

        # 4-base sequence and context
        rseq_4 = substr(a1,35,4)
        if (rseq_4~/[CG]CG[CG]/) { context = "SCGS" }
        else if (rseq_4~/[AT]CG[AT]/) { context = "WCGW" }
        else { context = "SCGW" }

        n_cpgs = gsub(/CG/,"",a1)                # number of CpGs in string
        gc_per = gsub(/[CG]/,"",a2) / length(a3) # GC-content of string
        rseq_6 = substr(a3,34,6)                 # 6-base sequence
        rseq_2 = substr(a3,36,2)                 # 2-base sequence

        # Need to have an easy way to extract CpGs into specific context and number of neighboring CpGs files
        if (n_cpgs == 1) { tag = "0" }
        else if (n_cpgs == 2) { tag = "1" }
        else if (n_cpgs == 3) { tag = "2" }
        else { tag = "3" }

        if (!(rseq_4~/N/)) print chr, beg, end, n_cpgs-1, gc_per, rseq_4, rseq_6, context, rseq_2, context"_"tag
    }' $1;
}
export -f find_context

create_files() {
    check_path

    if [ ! -d ${outdir} ]; then
        mkdir -p ${outdir}
    fi

    cd ${outdir}

    # Split up CpGs for quicker processing
    zcat ${cpgbed} | split -d -l 1000000 - cpgsplit_

    # Retrieve sequences around CpGs from genome FASTA file
    export genome
    parallel -j ${thread} 'expand_windows ${genome}.fai {} | bedtools getfasta -bed - -fi ${genome} -tab -fo cpgcontext_{}.tab' ::: cpgsplit_*

    # Annotate CpGs
    parallel -j ${thread} 'find_context {} > {.}.context' ::: cpgcontext_*.tab

    # Extract CpGs into specific context and number of neighboring CpG files
    grep --no-filename WCGW_0 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > wcgw_0_neighbors.bed.gz
    grep --no-filename WCGW_1 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > wcgw_1_neighbors.bed.gz
    grep --no-filename WCGW_2 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > wcgw_2_neighbors.bed.gz
    grep --no-filename WCGW_3 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > wcgw_3p_neighbors.bed.gz
    grep --no-filename SCGW_0 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgw_0_neighbors.bed.gz
    grep --no-filename SCGW_1 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgw_1_neighbors.bed.gz
    grep --no-filename SCGW_2 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgw_2_neighbors.bed.gz
    grep --no-filename SCGW_3 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgw_3p_neighbors.bed.gz
    grep --no-filename SCGS_0 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgs_0_neighbors.bed.gz
    grep --no-filename SCGS_1 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgs_1_neighbors.bed.gz
    grep --no-filename SCGS_2 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgs_2_neighbors.bed.gz
    grep --no-filename SCGS_3 *.context | cut -f1-9 | sort -k1,1 -k2,2n -k3,3n | gzip > scgs_3p_neighbors.bed.gz

    rm cpgsplit_* cpgcontext_*
}

usage() {
    >&2 echo -e "\nUsage: create_context_beds.sh [-h,--help] [-o,--outdir] genome cpgbed\n"
    >&2 echo -e "Required inputs:"
    >&2 echo -e "\tgenome       : Path to reference FASTA file used in creating CpG BED"
    >&2 echo -e "\tcpgbed       : CpG BED file containing CpG locations in genome\n"
    >&2 echo -e "Optional inputs:"
    >&2 echo -e "\t-h,--help    : Print help message and exit"
    >&2 echo -e "\t-o,--outdir  : Output directory [DEFAULT: assets]"
    >&2 echo -e "\t-t,--threads : Number of threads to use [DEFAULT: 1]"
}

# Initialize default variable values
outdir="assets"
thread=1

# Process command line arguments
OPTS=$(getopt \
    --options ho:t: \
    --long help,outdir:,threads: \
    --name "$(basename "$0")" \
    -- "$@"
)
eval set -- ${OPTS}

while true; do
    case "$1" in
        -h|--help )
            usage
            exit 0
            ;;
        -o|--outdir )
            outdir="$2"
            shift 2
            ;;
        -t|--threads )
            thread="$2"
            shift 2
            ;;
        -- )
            shift
            break
            ;;
        * )
            >&2 echo "Unknown option: $1"
            usage
            exit 1
            ;;
    esac
done

# Make sure there are the correct number of inputs
if [[ $# -ne 2 ]]; then
    >&2 echo "$0: Missing inputs"
    usage
    exit 1
fi

# Fill required positional arguments
genome=$1
cpgbed=$2

# Input checks
if [[ ! -f "${genome}.fai" ]]; then
    >&2 echo "Cannot locate fai-indexed reference: ${genome}.fai"
    >&2 echo "Please provide a viable path to the reference genome FASTA file."
    exit 1
fi

if [[ ! -f "${cpgbed}" ]]; then
    >&2 echo "Cannot locate CpG BED file: ${cpgbed}"
    >&2 echo "Please provide an existing CpG BED file"
    exit 1
fi

>&2 echo -e "Inputs:"
>&2 echo -e "\tOutput directory : ${outdir}"
>&2 echo -e "\tNumber of threads: ${thread}"
>&2 echo -e "\tGenome Reference : ${genome}"
>&2 echo -e "\tCpG Location BED : ${cpgbed}"

create_files

set -euo pipefail



# Find binned average values
function find_averages {
    tmp=$(echo $RANDOM|md5sum|head -c 20; echo) # get random name for tmp file
    if [[ `file ${infile}` =~ "gzip" ]]; then
        zcat ${infile} |
        awk -v filt=${filter} '{ if ($5 >= filt) { print } }' > $tmp
    else
        awk -v filt=${filter} '{ if ($5 >= filt) { print } }' \
        ${infile} > $tmp
    fi

    # Find average beta value across bins
    bedtools makewindows -w ${window} \
        -g ${reffai} | \
    sort -k1,1 -k2,2n | \
    bedtools map \
        -prec 2 \
        -a - \
        -b $tmp \
        -c 4 -o mean |
    gzip > ${otfile}.gz

    rm -f $tmp
}

# Print helpful usage information
function usage {
    >&2 echo -e "\nUsage: find_binned_average.sh [-h,--help] reference window filter infile outfile"
    >&2 echo -e "Required inputs:"
    >&2 echo -e "\treference    : FAI-index for reference to create windows from"
    >&2 echo -e "\twindow       : Window size for binned average finding"
    >&2 echo -e "\tfilter       : Coverage filter to apply before finding averages"
    >&2 echo -e "\tinfile       : BISCUIT CG BED file to calculate methylation average from"
    >&2 echo -e "\toutfile      : Name of BED file to write output to"
    >&2 echo -e "Optional inputs:"
    >&2 echo -e "\t-h,--help    : Print help message and exit\n"
}

# Process command line arguments
OPTS=$(getopt \
    --options h \
    --long help \
    --name "$(basename "$0")" \
    -- "$@"
)
eval set -- ${OPTS}

while true; do
    case "$1" in
        -h|--help )
            usage
            exit 0
            ;;
        -- )
            shift
            break
            ;;
        * )
            >&2 echo "Unknown option: $1"
            usage
            exit 1
            ;;
    esac
done

# Make sure there are the correct number of inputs
if [[ $# -ne 5 ]]; then
    >&2 echo "$0: Missing inputs"
    usage
    exit 1
fi

# Fill required positional arguments
reffai="${1}"
window="${2}"
filter="${3}"
infile="${4}"
otfile="${5}"

# Do some checks on the given inputs
if [[ ! -f ${reffai} ]]; then
    >&2 echo "Doesn't appear like you provided an existing FAI file: ${reffai}"
    exit 1
fi

if [[ ! "${reffai}" =~ .fai$ ]]; then
    >&2 echo "Doesn't appear like you provided a viable FAI file: ${reffai}"
    exit 1
fi

if ! [[ "${window}" =~ ^[0-9]+$ ]]; then
    >&2 echo "window needs to be an integer: ${window}"
    exit 1
fi

if ! [[ ${filter} =~ ^[0-9]+$ ]]; then
    >&2 echo "filter needs to be an integer: ${filter}"
    exit 1
fi

if [[ ! -f ${infile} ]]; then
    >&2 echo "Doesn't appear that infile exists: ${infile}"
    exit 1
fi

find_averages