#!/bin/bash
#SBATCH -p short
#SBATCH -t 0-8:00
#SBATCH --mem 10G
#SBATCH -e mergebam.%j.err

set -eo pipefail
shopt -s nullglob


module load conda/miniforge3/24.11.3-0
conda activate /n/groups/neuroduo/Shon/venv/bioenv

module load ucsc-tools/475

############################################################
# Check whether a BigWig exists and is structurally valid
############################################################

valid_bigwig () {

    file="$1"

    # File must exist and be non-empty
    if [[ ! -s "$file" ]]; then
        return 1
    fi

    # bigWigInfo must successfully parse it
    if bigWigInfo "$file" > /dev/null 2>&1; then
        return 0
    else
        return 1
    fi
}


############################################################
# Average a group of replicate BigWigs
############################################################

average_group () {

    output="$1"
    shift

    files=("$@")

    echo
    echo "============================================================"
    echo "Checking:"
    echo "  $output"

    ########################################################
    # If output is complete, skip it
    ########################################################

    if valid_bigwig "$output"; then
        echo "COMPLETE: valid BigWig already exists."
        echo "Skipping."
        return
    fi


    ########################################################
    # If an incomplete/corrupt output exists, remove it
    ########################################################

    if [[ -e "$output" ]]; then
        echo "INCOMPLETE/INVALID output found:"
        echo "  $output"
        echo "Removing it before rerunning."
        rm -f "$output"
    else
        echo "Output does not exist yet."
    fi


    ########################################################
    # Check input replicate count
    ########################################################

    echo "Found ${#files[@]} replicate files:"
    printf '  %s\n' "${files[@]}"

    if [[ ${#files[@]} -lt 2 ]]; then
        echo "WARNING: fewer than 2 replicates found."
        echo "Skipping."
        return
    fi


    ########################################################
    # Run bigwigAverage
    ########################################################

    echo
    echo "Running bigwigAverage..."
    echo "Output:"
    echo "  $output"

    bigwigAverage \
        -b "${files[@]}" \
        -o "$output"


    ########################################################
    # Validate newly generated output
    ########################################################

    if valid_bigwig "$output"; then
        echo "SUCCESS: valid merged BigWig created:"
        echo "  $output"
    else
        echo "ERROR: merged output failed validation:"
        echo "  $output"

        rm -f "$output"
        return 1
    fi
}


############################################################
# H1C1b F10-1
############################################################

for species in Human Chimp; do

    files=(
        ${species}_piN_H1C1b_F10-1_2hrKCl_IgG_R*.hg38.bw
    )

    average_group \
        "${species}_piN_H1C1b_F10-1_2hrKCl_IgG_merged.hg38.bw" \
        "${files[@]}"


    files=(
        ${species}_piN_H1C1b_F10-1_minusKCl_IgG_R*.hg38.bw
    )

    average_group \
        "${species}_piN_H1C1b_F10-1_minusKCl_IgG_merged.hg38.bw" \
        "${files[@]}"

done


############################################################
# H2C2a G6-2
############################################################

for species in Human Chimp; do

    files=(
        ${species}_piN_H2C2a_G6-2_2hrKCl_IgG_R*.hg38.bw
    )

    average_group \
        "${species}_piN_H2C2a_G6-2_2hrKCl_IgG_merged.hg38.bw" \
        "${files[@]}"


    files=(
        ${species}_piN_H2C2a_G6-2_minusKCl_IgG_R*.hg38.bw
    )

    average_group \
        "${species}_piN_H2C2a_G6-2_minusKCl_IgG_merged.hg38.bw" \
        "${files[@]}"

done


############################################################
# H2C2a G11-1
#
# Some files use G11-1 and others use G11_1
############################################################

for species in Human Chimp; do

    files=(
        ${species}_piN_H2C2a_G11-1_2hrKCl_IgG_R*.hg38.bw
        ${species}_piN_H2C2a_G11_1_2hrKCl_IgG_R*.hg38.bam.bw
    )

    average_group \
        "${species}_piN_H2C2a_G11-1_2hrKCl_IgG_merged.hg38.bw" \
        "${files[@]}"


    files=(
        ${species}_piN_H2C2a_G11-1_minusKCl_IgG_R*.hg38.bw
        ${species}_piN_H2C2a_G11_1_minusKCl_IgG_R*.hg38.bam.bw
    )

    average_group \
        "${species}_piN_H2C2a_G11-1_minusKCl_IgG_merged.hg38.bw" \
        "${files[@]}"

done


echo
echo "============================================================"
echo "Finished checking and processing all IgG groups."