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

DATA_DIR=/home/ddpanchenko/MPRA/TOMSK_HEK293/data
OUT_DIR=/home/ddpanchenko/MPRA/TOMSK_HEK293/alignments
REF_DIR=/home/ddpanchenko/MPRA/TOMSK_HEK293/data/references
THREADS=8

SAMPLES=(HEK293  HEK293T  Viral_production_cell)
REFERENCES=(hg38  hg38_add_seq)

mkdir -p "$OUT_DIR"

for sample in "${SAMPLES[@]}"; do
  echo "== Sample: ${sample}"

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

  mkdir -p RawReads FilteredReads logs
  mapfile -t files < <(find "${DATA_DIR}/${sample}" -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 адаптеров
  echo "Trimming adapters"
  cutadapt \
    -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \
    -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
    -m 50 \
    -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

  for ref in "${REFERENCES[@]}"; do
    reference_dir="${REF_DIR}/${ref}"
    SIZES="${reference_dir}/chrom.sizes"
    REF="${reference_dir}/${ref}.fna"

    if [[ ! -d "$reference_dir" ]]; then
      echo "❌ Reference directory not found: ${reference_dir}" >&2
      continue
    fi
    
    echo "   → Reference: ${REF}"
    echo "   → Chrom sizes: ${SIZES}"

    # 4) bwa mem "${R1_CLEAN}" "${R2_CLEAN}" FilteredReads/R1.cln.fq.gz FilteredReads/R2.cln.fq.gz
    bwa mem -t ${THREADS} "${REF}" FilteredReads/R1.cln.fq.gz FilteredReads/R2.cln.fq.gz \
      | samtools view -b \
      | samtools sort -n -@ ${THREADS} -o "${sample}.${ref}.name_sorted.bam"
    # samtools index "${sample}.${ref}.sorted.bam"

    ### проверка bam файла
    bam_file="${sample}.${ref}.name_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

    samtools fixmate -m "${sample}.${ref}.name_sorted.bam" "${sample}.${ref}.fixmate.bam"
    samtools sort -@ ${THREADS} -o "${sample}.${ref}.coord_sorted.bam" "${sample}.${ref}.fixmate.bam"
    samtools markdup -r "${sample}.${ref}.coord_sorted.bam" "${sample}.${ref}.dedup.bam"
    samtools index "${sample}.${ref}.dedup.bam"

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

    # 6) MAPQ≥30 → BigWig (-f 0x2 filter only properly paird) and IS <= 1000
    samtools view -@ "${THREADS}" -h -q 30 -f 0x2 "${sample}.${ref}.dedup.bam" \
    | awk 'BEGIN{FS=OFS="\t"} /^@/{print; next} {t=$9; if(t<0)t=-t; if(t<=10000000) print}' \
    | samtools sort -@ "${THREADS}" -o "${sample}.${ref}.q30.bam"
    samtools index -@ "${THREADS}" "${sample}.${ref}.q30.bam"

    if [[ $(samtools view -c "${sample}.${ref}.q30.bam") -eq 0 ]]; then
      echo "⚠️  BAM file after MQ for ${sample} is empty, skipping..."
      continue
    fi

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

    bedGraphToBigWig "${sample}.${ref}.q30.bedGraph" "${SIZES}" "${sample}.${ref}.q30.bw"
  done
done

echo "DONE"
