#!/usr/bin/env bash
# End-to-end SimurgArray run on the public GSA v4 PGx demo set (36 samples):
#   IDAT → GTC + sample QC → HTML QC report → VCF (per sample and cohort BCF) → PLINK / bedGraph / final report.
#
# Prerequisites (from the repository root):
#   uv venv && uv pip install -e .
#   bash scripts/fetch_data.sh p0      # demo IDATs + manifest + cluster file (~1.4 GB)
#   bash scripts/fetch_data.sh p6      # GRCh38 FASTA for VCF output (~0.9 GB download)
#
# Usage: bash examples/gsa_demo.sh [OUT_DIR] [N_SAMPLES]      (defaults: results/gsa_demo, all 36)
# Takes ~10 min for 36 samples on a Raspberry Pi 5 (2 worker processes); much less on a desktop.
set -euo pipefail
cd "$(dirname "$0")/.."
OUT=${1:-results/gsa_demo}
N=${2:-0}
THREADS=${THREADS:-2}
SIMURG=${SIMURG:-.venv/bin/simurg}

D=data/p0/gsa_v4_pgx
BPM=$D/GSA-PGx-48v4-0_20079540_E2_bpm/GSA-PGx-48v4-0_20079540_E2.bpm
CSV=$D/GSA-PGx-48v4-0_20079540_E2_csv/GSA-PGx-48v4-0_20079540_E2.csv
EGT=$D/GSA-PGx-48v4-0_20079540_E2_ClusterFile/GSA-PGx-48v4-0_20079540_E2_ClusterFile.egt
IDATS=$D/idats_Demo_36/idats_Demo_36
FASTA=data/p6/GRCh38.primary.fa
for f in "$BPM" "$CSV" "$EGT" "$IDATS" "$FASTA"; do
  [[ -e $f ]] || { echo "missing $f — run scripts/fetch_data.sh p0 and p6 first"; exit 1; }
done
mkdir -p "$OUT"

# Optional subset: symlink the first N IDAT pairs into $OUT/idats (same <barcode>/<barcode>_<pos>_<ch>.idat layout).
if [[ $N -gt 0 ]]; then
  rm -rf "$OUT/idats"
  find "$IDATS" -name '*_Grn.idat' | sort | head -n "$N" | while read -r g; do
    b=$(basename "$(dirname "$g")"); mkdir -p "$OUT/idats/$b"
    ln -s "$(realpath "$g")" "$OUT/idats/$b/"; ln -s "$(realpath "${g%_Grn.idat}_Red.idat")" "$OUT/idats/$b/"
  done
  IDATS=$OUT/idats
fi

echo "== 1. inspect the inputs"
$SIMURG inspect "$BPM" "$EGT" | head -20

echo "== 2. genotype: GTC v5 per sample, sample_summary.tsv, control-probe QC, contamination, relatedness"
$SIMURG genotype --bpm "$BPM" --egt "$EGT" --csv "$CSV" --idats "$IDATS" --out "$OUT/gtc" \
  --qc-profile genotyping --threads "$THREADS" --force
column -t -s$'\t' "$OUT/gtc/sample_summary.tsv" | cut -c1-150 | head -8

echo "== 3. HTML QC report + QC table"
$SIMURG qc report --data "$OUT/gtc" --out "$OUT/qc_report" --output-format both

echo "== 4. VCF: one bgzipped, tabix-indexed VCF per sample, then one cohort BCF"
$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/vcf" \
  --threads "$THREADS" --force
$SIMURG gtc-to-vcf --bpm "$BPM" --csv "$CSV" --fasta "$FASTA" --gtcs "$OUT/gtc" --out "$OUT/bcf" --merge \
  --threads "$THREADS" --force

echo "== 5. exports"
$SIMURG export plink --vcfs "$OUT/vcf" --gtcs "$OUT/gtc" --out "$OUT/plink" --format both --force
$SIMURG export bedgraph --bpm "$BPM" --gtcs "$OUT/gtc" --out "$OUT/bedgraph" --chrom-style ucsc --force
$SIMURG export finalreport --bpm "$BPM" --csv "$CSV" --egt "$EGT" --gtcs "$OUT/gtc" --out "$OUT/finalreport" --force

echo "== 6. differential test: re-call with the older normalization and compare field by field"
$SIMURG genotype --bpm "$BPM" --egt "$EGT" --idats "$IDATS" --out "$OUT/gtc_norm112" --norm-version 1.1.2 \
  --no-controls --no-cohort-qc --threads "$THREADS" --quiet --force
$SIMURG compare "$OUT/gtc" "$OUT/gtc_norm112" --bpm "$BPM" --out "$OUT/compare" --title "norm 1.2.0 vs 1.1.2" \
  --threads "$THREADS" --force || true

echo
echo "Done. Open $OUT/qc_report/*.html and $OUT/compare/dashboard.html in a browser."
find "$OUT" -maxdepth 2 -type f \( -name '*.html' -o -name '*.bcf' -o -name '*.pgen' -o -name '*.bed' \) | sort