#!/bin/bash
# Phase P8 scientific validation on GEO GSE165845 (GSA-MD v1, 120 1000 Genomes individuals × 3 replicates):
#   1. fetch IDATs (GEO), GSA v1 C2 build38 BPM/CSV + C1 EGT, GRCh38 FASTA, NYGC 30x per-chromosome VCFs, GIAB HG001;
#   2. `simurg genotype` → `simurg gtc-to-vcf --merge`;
#   3. truth subsets at the array positions (bcftools, split + left-normalised like the array records);
#   4. scripts/p8/validate_p8.py (plan §16.5 metrics).
# Args: <unused data dir>  <results dir>. Run with OMNI=0 (the Omni2.5 archive is not needed).
set -uo pipefail
OUT=$2
W=/work/p8; mkdir -p $W $OUT/p8
R=$OUT/p8
N=$(nproc)
export DEBIAN_FRONTEND=noninteractive
apt-get install -y -q bcftools tabix > /dev/null
bcftools --version | head -1

GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
NYGC=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV
GIAB=https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/NISTv4.2.1/GRCh38/HG001_GRCh38_1_22_v4.2.1_benchmark
CHROMS="$(seq 1 22) X"

echo "=== downloads $(date -u +%T) ==="
mkdir -p $W/nygc $W/gsa $W/giab
( for e in vcf.gz vcf.gz.tbi bed; do curl -sSfL --retry 5 -o $W/giab/hg001.$e $GIAB.$e; done ) &
( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 &
FASTA_PID=$!
cd $W/gsa
for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
         infinium-global-screening-array-24-v1-0-c2-manifest-file-csv-build38.zip \
         infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
  curl -sSfL --retry 5 -O "$GSA/$z"
done
sha256sum -c - <<'EOF' || { echo "!! GSA product file checksum mismatch"; exit 1; }
10ba8e692f6f14f9c50b3bed00d92a96af7ef8aa873df115abc52b5ae7c209b8  infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip
fed9af6de6a8f000b42f44ff1aeecc9fde2a4ae29c834c542bedb2f28879a815  infinium-global-screening-array-24-v1-0-c2-manifest-file-csv-build38.zip
8bd33c286feded0e77b93cacc8533c27dd879a06e97ef9a575dff2037348cb67  infinium-global-screening-array-v1-0-c1-cluster-file.zip
EOF
for z in *.zip; do unzip -oq $z; done
cd - > /dev/null
BPM=$W/gsa/GSA-24v1-0_C2.bpm; CSV=$W/gsa/GSA-24v1-0_C2.csv; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
ls -la $W/gsa

# NYGC per-chromosome VCFs (~40 GB, EBI serves ~15 MB/s) are streamed and cut down on the fly to the manifest
# positions (±100 bp, so left-normalised indels stay inside) and the GSE165845 individuals: no copy on disk, and
# the CPU stays busy (a download-only phase would trip the idle alarm, which terminates the instance).
.venv/bin/python - $BPM > $W/targets_chr.tsv <<'PY'
import sys
from simurg.formats.bpm import read_bpm
m = read_bpm(sys.argv[1])
iv = {}
for c, p in zip(m.chroms, m.map_infos.tolist()):
    if p > 0 and c in {str(i) for i in range(1, 23)} | {"X"}:
        iv.setdefault(c, []).append((max(1, p - 100), p + 100))
for c in [str(i) for i in range(1, 23)] + ["X"]:
    cur = None
    for s, e in sorted(iv.get(c, [])):
        if cur and s <= cur[1]:
            cur[1] = max(cur[1], e)
            continue
        if cur:
            print(f"chr{c}\t{cur[0]}\t{cur[1]}")
        cur = [s, e]
    if cur:
        print(f"chr{c}\t{cur[0]}\t{cur[1]}")
PY
wc -l < $W/targets_chr.tsv
for c in $(seq 1 22) X Y M; do echo -e "chr$c\t$c"; done > $W/chrmap.txt
cut -f2 scripts/p8/gse165845_samples.tsv | tail -n +2 | sort -u | paste -sd, > $W/inds.txt
export W NYGC
( for c in $CHROMS; do echo $c; done | xargs -P 23 -I{} sh -c \
    's=""; [ {} = X ] && s=".v2"
     curl -sSfL --retry 5 "$NYGC/1kGP_high_coverage_Illumina.chr{}.filtered.SNV_INDEL_SV_phased_panel$s.vcf.gz" \
     | bcftools view -T $W/targets_chr.tsv -s $(cat $W/inds.txt) --force-samples -Ou \
     | bcftools annotate --rename-chrs $W/chrmap.txt -Ob -o $W/nygc/raw.{}.bcf 2> $W/nygc/{}.log \
     && echo "chr{} done $(date -u +%T)" || echo "!! chr{} failed"'
  echo "NYGC streamed at $(date -u +%T)" ) > $W/nygc.log 2>&1 &
NYGC_PID=$!

/usr/bin/time -v .venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo \
  --out $W/idats --workers 24 2> $R/prepare_time.txt || { tail -30 $R/prepare_time.txt; exit 1; }
cp $W/idats/samples_map.tsv $R/

echo "=== genotype $(date -u +%T) ==="
GT=$W/gt
/usr/bin/time -v .venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $GT \
  --threads "$N" --quiet 2> $R/genotype_time.txt || echo "genotype exit $?"
tail -22 $R/genotype_time.txt
cp $GT/sample_summary.tsv $GT/relatedness.tsv $GT/controls_qc.tsv $GT/run_manifest.json $R/ 2>/dev/null
cp -r $GT/logs $R/genotype_logs 2>/dev/null

wait $FASTA_PID; tail -2 $W/fasta.log
FA=data/p6/GRCh38.primary.fa
echo "=== gtc-to-vcf $(date -u +%T) ==="
/usr/bin/time -v .venv/bin/simurg gtc-to-vcf --bpm $BPM --csv $CSV --fasta $FA --gtcs $GT --out $W/vcf --merge \
  --threads "$N" --quiet 2> $R/vcf_time.txt || echo "gtc-to-vcf exit $?"
tail -22 $R/vcf_time.txt
cp $W/vcf/run_manifest.json $R/vcf_run_manifest.json; cp -r $W/vcf/logs $R/vcf_logs 2>/dev/null
ARR=$W/vcf/cohort.snv.bcf
[[ -f $ARR ]] || { echo "!! no cohort BCF"; exit 1; }

echo "=== truth subsets $(date -u +%T) ==="
bcftools query -f '%CHROM\t%POS\n' $ARR | awk '{print "chr"$1"\t"$2}' | uniq > $W/sites_chr.tsv
wc -l < $W/sites_chr.tsv
bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/array.norm.bcf $ARR 2> $R/norm_array.txt
wait $NYGC_PID; cat $W/nygc.log; df -h /work | tail -1
export FA
echo $CHROMS | tr ' ' '\n' | xargs -P "$N" -I{} sh -c \
  'bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/truth.{}.bcf $W/nygc/raw.{}.bcf 2> $W/truth.{}.log \
   || echo "!! truth chr{} failed"'
bcftools concat -Ob -o $W/truth.bcf $(for c in $CHROMS; do echo $W/truth.$c.bcf; done)
bcftools view -T $W/sites_chr.tsv -Ou $W/giab/hg001.vcf.gz | bcftools annotate --rename-chrs $W/chrmap.txt -Ou \
  | bcftools norm -f $FA -m-any --check-ref w -Ob -o $W/giab.bcf
echo "array $(bcftools view -H $W/array.norm.bcf | wc -l) records, truth $(bcftools view -H $W/truth.bcf | wc -l)," \
  "giab $(bcftools view -H $W/giab.bcf | wc -l)"

echo "=== validate $(date -u +%T) ==="
/usr/bin/time -v .venv/bin/python scripts/p8/validate_p8.py --array $W/array.norm.bcf --truth $W/truth.bcf \
  --giab $W/giab.bcf --giab-bed $W/giab/hg001.bed --map $W/idats/samples_map.tsv --gt-dir $GT --out $R \
  2> $R/validate_time.txt || echo "validation exit $?"
tail -22 $R/validate_time.txt
gzip -f $R/replicate_pairs.tsv $R/systematic_sites.tsv 2>/dev/null
echo "=== task done $(date -u +%T) ==="