#!/usr/bin/env bash
set -euo pipefail

DATA_DIR=/data/ddpanchenko/main_dom/ChiPseq_25122025/data
OUT_DIR=/data/ddpanchenko/main_dom/ChiPseq_25122025/alignments
THREADS=8

GENOME="/data/ddpanchenko/main_dom/genomes/hg38/hg38.fna"
GENOME_INDEX="/data/ddpanchenko/main_dom/genomes/hg38/bowtie_index"

# GENOME="/data/ddpanchenko/main_dom/hybrid_anopheles_human/Repeats/TandemRepeats/circ19-98/ref/prefix_circ19-98_monomer.fa"
# GENOME_INDEX="/data/ddpanchenko/main_dom/hybrid_anopheles_human/Repeats/TandemRepeats/circ19-98/ref/bowtie_index"

# WINDOW="/home/ddpanchenko/hybrid_anopheles_human/ChiPseq/aligned_bam_bowtie/windows.bed"
CHROM_SIZES="/data/ddpanchenko/main_dom/genomes/hg38/chrom.sizes"

# MAP=(H3K27Ac1 H3K27Ac2 CTCF1_input CTCF2_input H3K27Ac1_input H3K27Ac2_input) # CTCF1 CTCF2
MAP=(H3K27Ac2_2BGI) # CTCF1 CTCF2

mkdir -p "$OUT_DIR"

for sample_name in "${MAP[@]}"; do
    sample_dir="$DATA_DIR/${sample_name}/"
    sample=$(basename "$sample_dir")

    echo "== Sample: ${sample}"

    SDIR="${OUT_DIR}/${sample}"
    mkdir -p "${SDIR}"
    cd "${SDIR}"

    mkdir -p RawReads FilteredReads logs
    mapfile -t files < <(
    find "$sample_dir" -maxdepth 1 -type f \
        \( -name "${sample}_*_1.f*q.gz" -o -name "${sample}_*_2.f*q.gz" \
        -o -name "${sample}_*R1*.f*q.gz" -o -name "${sample}_*R2*.f*q.gz" \) \
    | sort
    )

    if [ ${#files[@]} -eq 0 ]; then
        echo "No FASTQ files found for sample ${sample} in ${DATA_DIR}/${sample}/"
    else
        for f in "${files[@]}"; do
        if [ -f "$f" ]; then
            ln -sf "$f" RawReads/
            echo "🔗 Linked: $(basename "$f")"
        else
            echo "Missing file: $f"
        fi
        done
    fi

    shopt -s nullglob
    R1=$(ls RawReads/${sample}_*_1.f*q.gz RawReads/${sample}_*R1*.f*q.gz 2>/dev/null | head -n1)
    R2=$(ls RawReads/${sample}_*_2.f*q.gz RawReads/${sample}_*R2*.f*q.gz 2>/dev/null | head -n1)
    shopt -u nullglob
    echo "    R1='$R2'"

    # Проверка, что оба файла найдены
    if [[ -z "$R1" || -z "$R2" ]]; then
        echo "Could not find both R1 and R2 for ${sample} in RawReads/"
        echo "    R1='$R1'"
        echo "    R2='$R2'"
    fi

    # 1) trim 3' адаптеров
    echo "Triming adapters for ${sample}"
    cutadapt \
        -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
        -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
        -m 70 \
        -o FilteredReads/R1.trim.fq.gz -p FilteredReads/R2.trim.fq.gz \
        "${R1}" "${R2}" \
        --report=full > logs/cutadapt_trim3_full_report.txt

    cutadapt --cores 4 \
      -g ^TGACTTatcatgtctgctcgaagcGTCGAGCTAAGTCGTATGGTATGCGCAATTACATTCAGCT \
      -g ^ATGGTATGCGCAATTACATTCAGCT \
      -a AGCTGAATGTAATTGCGCATACCAT \
      -A AGCTGAATGTAATTGCGCATACCAT \
      -e 0.2 --overlap 10 \
      -o FilteredReads/R1.cln.fq.gz -p FilteredReads/R2.cln.fq.gz \
      FilteredReads/R1.trim.fq.gz FilteredReads/R2.trim.fq.gz \
      --report=full > logs/cutadapt_filter5_report.txt

    # 2) bowtie "${R1_CLEAN}" "${R2_CLEAN}" FilteredReads/R1.cln.fq.gz FilteredReads/R2.cln.fq.gz
    bowtie2 -p $THREADS -x "$GENOME_INDEX" \
        -1 FilteredReads/R1.cln.fq.gz \
        -2 FilteredReads/R2.cln.fq.gz \
        | samtools view -b \
        | samtools sort -@ "${THREADS}" -o "${sample}.sorted.bam"
    samtools index "${sample}.sorted.bam"

    ### проверка bam файла
    bam_file="${sample}.sorted.bam"

    if [[ ! -s "$bam_file" ]]; then
        echo "⚠️  No BAM file generated for ${sample}, skipping..."
        continue
    fi

    # Проверяем, что BAM содержит хотя бы одно выравнивание
    if [[ $(samtools view -c "$bam_file") -eq 0 ]]; then
        echo "⚠️  BAM file for ${sample} is empty, skipping..."
        continue
    fi

    # 5) dedup по UMI
    umi_tools dedup \
      --paired --method=unique --ignore-umi \
      -I "${sample}.sorted.bam" \
      -S "${sample}.dedup.bam" \
      --log="logs/${sample}.dedup.log"
    samtools index "${sample}.dedup.bam"

    # cp "${sample}.circ19-98.sorted.bam" "${sample}.circ19-98.dedup.bam"

    # Проверяем, что BAM содержит хотя бы одно выравнивание
    if [[ $(samtools view -c "${sample}.dedup.bam") -eq 0 ]]; then
        echo "⚠️  BAM file after umi_tools dedup for ${sample} is empty, skipping..."
        continue
    fi

    # 6) MAPQ≥30 → BigWig 
    samtools view -@ "${THREADS}" -h -q 30 -f 0x2 "${sample}.dedup.bam" \
    | samtools sort -@ "${THREADS}" -o "${sample}.q30.bam"
    samtools index -@ "${THREADS}" "${sample}.q30.bam"


    bedtools genomecov -ibam "${sample}.q30.bam" -bg \
        > "${sample}.q30.bedGraph" 

    sort -k1,1 -k2,2n "${sample}.q30.bedGraph" -o "${sample}.q30.sorted.bedGraph"
    bedGraphToBigWig "${sample}.q30.sorted.bedGraph" "${CHROM_SIZES}" "${sample}.q30.bw"
    
    rm "${sample}.q30.sorted.bedGraph"
done

echo "DONE"
