#!/usr/bin/env bash
# 01_align_call_peaks.sh
# Per-sample fastq -> alignment -> dedup -> ATAC-seq peak calling
#
# Usage:
#   ./01_align_call_peaks.sh <SRR_ID> <LINE_NAME> <DISEASE>
# Example:
#   ./01_align_call_peaks.sh SRR11735719 HGMDFN090 control
#
# Pipeline:
#   1. bowtie2 paired-end alignment to GRCh37/hg19 (matches GSE150136 Eq2018 protocol)
#   2. samtools sort + remove duplicates (markdup, then -F 1024)
#   3. Filter for properly paired, MAPQ >= 30, primary alignments
#   4. Remove mitochondrial reads (chrM is over-amplified in ATAC-seq)
#   5. macs3 callpeak --nomodel --shift -75 --extsize 150 -f BAMPE (canonical
#      ATAC-seq narrow-peak calling per ENCODE recipe)
#
# Output:
#   <line>.dedup.bam   -- final deduplicated, MAPQ-filtered BAM
#   <line>_peaks.narrowPeak  -- input for chromperiod
#
# All steps stream output incrementally to keep memory low.

set -euo pipefail

SRR=$1
LINE=$2
DISEASE=$3

ROOT="/Users/tonesonvox/Projects/Wavelet/hgps_atac_pipeline_2026-05-02"
INDEX="${ROOT}/data/reference/GRCh37/GRCh37"   # bowtie2 index prefix
FASTQ_DIR="${ROOT}/data/fastq"
BAM_DIR="${ROOT}/data/bam"
PEAK_DIR="${ROOT}/data/peaks"
LOG_DIR="${ROOT}/logs"
THREADS=${THREADS:-8}

mkdir -p "$BAM_DIR" "$PEAK_DIR" "$LOG_DIR"

R1="${FASTQ_DIR}/${SRR}_1.fastq.gz"
R2="${FASTQ_DIR}/${SRR}_2.fastq.gz"
SAM="${BAM_DIR}/${LINE}.sam"
RAW="${BAM_DIR}/${LINE}.raw.bam"
SORTED="${BAM_DIR}/${LINE}.sorted.bam"
DEDUP="${BAM_DIR}/${LINE}.dedup.bam"
PEAKS_PREFIX="${PEAK_DIR}/${LINE}"
LOG="${LOG_DIR}/${LINE}.alignment.log"

if [[ ! -f "$R1" || ! -f "$R2" ]]; then
  echo "ERROR: fastq files missing: $R1 / $R2" >&2
  exit 1
fi

echo "=== ${LINE} (${SRR}, ${DISEASE}) ==="

# 1. Alignment
if [[ ! -f "$SORTED" ]]; then
  echo "[$(date +%H:%M:%S)] bowtie2 alignment..."
  bowtie2 --very-sensitive -X 2000 -p "$THREADS" \
    -x "$INDEX" -1 "$R1" -2 "$R2" -S "$SAM" 2>"$LOG"
  echo "[$(date +%H:%M:%S)] sort..."
  samtools view -@ "$THREADS" -bS -q 30 -f 2 -F 256 "$SAM" \
    | samtools sort -@ "$THREADS" -o "$SORTED" -
  rm -f "$SAM"
  samtools index "$SORTED"
fi

# 2. Remove duplicates + chrM
if [[ ! -f "$DEDUP" ]]; then
  echo "[$(date +%H:%M:%S)] markdup + dedup + chrM filter..."
  samtools sort -n -@ "$THREADS" "$SORTED" \
    | samtools fixmate -m -@ "$THREADS" - - \
    | samtools sort -@ "$THREADS" - \
    | samtools markdup -@ "$THREADS" - - \
    | samtools view -@ "$THREADS" -b -F 1024 - \
    | samtools view -@ "$THREADS" -b -e 'rname != "chrM" && rname != "MT"' - \
      > "$DEDUP"
  samtools index "$DEDUP"
fi

# 3. Peak call (ATAC-seq canonical params)
if [[ ! -f "${PEAKS_PREFIX}_peaks.narrowPeak" ]]; then
  echo "[$(date +%H:%M:%S)] macs3 callpeak..."
  macs3 callpeak \
    -t "$DEDUP" \
    -f BAMPE \
    -g hs \
    -n "$LINE" \
    --outdir "$PEAK_DIR" \
    --nomodel --shift -75 --extsize 150 \
    --keep-dup all \
    --call-summits \
    -q 0.01 \
    2> "${LOG_DIR}/${LINE}.macs3.log"
fi

# 4. Quick QC
NPEAKS=$(wc -l < "${PEAKS_PREFIX}_peaks.narrowPeak")
DEDUPED=$(samtools view -c "$DEDUP")
echo "[$(date +%H:%M:%S)] DONE  reads=${DEDUPED}  peaks=${NPEAKS}"

# 5. Run chromperiod immediately on the freshly called peaks
echo "[$(date +%H:%M:%S)] chromperiod CWT..."
mkdir -p "${ROOT}/results/per_chrom" "${ROOT}/results/phase_maps"
python3 "${ROOT}/scripts/02_run_chromperiod_per_sample.py" \
  --bed "${PEAKS_PREFIX}_peaks.narrowPeak" \
  --line "$LINE" \
  --disease "$DISEASE" \
  --out-summary "${ROOT}/results/per_chrom/${LINE}_per_chrom_summary.csv" \
  --out-phase-map "${ROOT}/results/phase_maps/${LINE}_phase_map_500kb.csv" \
  > "${LOG_DIR}/${LINE}.chromperiod.log" 2>&1
tail -20 "${LOG_DIR}/${LINE}.chromperiod.log"

echo "=== ${LINE} complete ==="
