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

IN_DIR="/data/ddpanchenko/main_dom/ChiPseq_25122025/alignments"
MACS_OUT="/data/ddpanchenko/main_dom/ChiPseq_25122025/macs_bowtie2_25122025"

CHROM_SIZES="/data/ddpanchenko/main_dom/genomes/hg38/chrom.sizes"
GENOME_SIZE="hs"

mkdir -p "$MACS_OUT"

# CHIPS=(
# CTCF1 CTCF2 H3K27Ac1 H3K27Ac2
# )

CHIPS=(
H3K27Ac2_2BGI
)

find_bam() {
  local d="$1"
  local b=""

  b="$(find "$d" -maxdepth 1 -type f \
      \( -name "*.q30.bam" -o -name "*.mapq30.bam" \) \
      ! -name "*circ19-98*" \
      | sort | head -n1 || true)"

  if [[ -z "$b" ]]; then
    b="$(find "$d" -maxdepth 1 -type f \
        -name "*.dedup.bam" \
        ! -name "*circ19-98*" \
        | sort | head -n1 || true)"
  fi

  if [[ -z "$b" ]]; then
    b="$(find "$d" -maxdepth 1 -type f \
        -name "*.sorted.bam" \
        ! -name "*circ19-98*" \
        | sort | head -n1 || true)"
  fi

  [[ -n "$b" ]] && printf '%s\n' "$b"
}

for chip in "${CHIPS[@]}"; do
  input="${chip}_input"

  chip_dir="${IN_DIR}/${chip}"
  input_dir="${IN_DIR}/${input}"

  if [[ ! -d "$chip_dir" ]]; then
    echo "⚠️  No dir for ChIP: $chip_dir — skip"
    continue
  fi
  if [[ ! -d "$input_dir" ]]; then
    echo "⚠️  No dir for input: $input_dir — skip"
    continue
  fi

  chip_bam="$(find_bam "$chip_dir" || true)"
  input_bam="$(find_bam "$input_dir" || true)"

  if [[ -z "$chip_bam" ]]; then
    echo "⚠️  No BAM found for ChIP $chip in $chip_dir — skip"
    continue
  fi
  if [[ -z "$input_bam" ]]; then
    echo "⚠️  No BAM found for input $input in $input_dir — skip"
    continue
  fi

  # убедимся, что есть индексы (полезно для speed/совместимости)
  [[ -f "${chip_bam}.bai" ]] || samtools index "$chip_bam"
  [[ -f "${input_bam}.bai" ]] || samtools index "$input_bam"

  sample="$chip"
  outdir="${MACS_OUT}/${sample}"
  mkdir -p "$outdir"

  echo "MACS3 callpeak: ChIP=$chip  Input=$input"
  macs3 callpeak \
    -t "$chip_bam" \
    -c "$input_bam" \
    -f BAM \
    -g "$GENOME_SIZE" \
    -n "$sample" \
    --outdir "$outdir" \
    --broad --broad-cutoff 0.1 \
    -B --SPMR

  # --- BigWig из treat pileup (как раньше) ---
  treat_bdg="${outdir}/${sample}_treat_pileup.bdg"
  if [[ -s "$treat_bdg" ]]; then
    sort -k1,1 -k2,2n "$treat_bdg" -o "${treat_bdg%.bdg}.sorted.bdg"
    bedGraphToBigWig "${treat_bdg%.bdg}.sorted.bdg" "$CHROM_SIZES" "${outdir}/${sample}.bw"
    rm -f "${treat_bdg%.bdg}.sorted.bdg"
  else
    echo "⚠️  No treat_pileup.bdg for $sample (expected: $treat_bdg)"
  fi

  # --- ДОБАВЛЕНО: logFE + control_lambda BigWig ---
  control_bdg="${outdir}/${sample}_control_lambda.bdg"
  control_bw="${outdir}/${sample}_control_lambda.bw"

  out_prefix="${sample}_lfc"
  out_bdg="${outdir}/${out_prefix}_FE.bdg"
  out_bw="${outdir}/${out_prefix}_FE.bw"

  if [[ -s "$treat_bdg" && -s "$control_bdg" ]]; then
    echo "MACS3 bdgcmp logFE for $sample"
    macs3 bdgcmp \
      -t "$treat_bdg" \
      -c "$control_bdg" \
      --outdir "$outdir" \
      --o-prefix "$out_prefix" \
      -m FE -p 1

    # macs3 может вернуть unsorted bedGraph — сортируем перед BigWig
    if [[ -s "$out_bdg" ]]; then
      sort -k1,1 -k2,2n "$out_bdg" -o "${out_bdg%.bdg}.sorted.bdg"
      bedGraphToBigWig "${out_bdg%.bdg}.sorted.bdg" "$CHROM_SIZES" "$out_bw"
      rm -f "${out_bdg%.bdg}.sorted.bdg" "$out_bdg"
    else
      echo "⚠️  Expected logFE bdg not found: $out_bdg"
    fi

    # control_lambda -> bw (тоже лучше сортировать)
    sort -k1,1 -k2,2n "$control_bdg" -o "${control_bdg%.bdg}.sorted.bdg"
    bedGraphToBigWig "${control_bdg%.bdg}.sorted.bdg" "$CHROM_SIZES" "$control_bw"
    rm -f "${control_bdg%.bdg}.sorted.bdg"

    echo "✅ Done $sample"
  else
    echo "❌ Skipping bdgcmp for $sample: required .bdg files not found."
    echo "    treat:   $treat_bdg"
    echo "    control: $control_bdg"
  fi

done

echo "DONE"
