#!/bin/bash
# v2 CN/LOH validation on GSE165845 (360 GSA arrays, 120 1000 Genomes individuals x 3): simurg cn call vs the NYGC
# 30x structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between replicates.
# Args: <unused data dir>  <results dir>. Run with OMNI=0.
set -uo pipefail
OUT=$2
W=/work/cnv; R=$OUT/cnv; mkdir -p $W $R
N=$(nproc)
export DEBIAN_FRONTEND=noninteractive
apt-get install -y -q bcftools tabix gcc > /dev/null
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal > /dev/null 2>&1
export PATH="$HOME/.cargo/bin:$PATH"
VIRTUAL_ENV=$PWD/.venv uv pip install -q ./rust/simurg-core || echo "!! kernels not built (pure Python)"

SV=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20210124.SV_Illumina_Integration/1KGP_3202.gatksv_svtools_novelins.freeze_V3.wAF.vcf.gz
( curl -sSfL --retry 5 -o $W/sv.vcf.gz $SV && curl -sSfL --retry 5 -o $W/sv.vcf.gz.tbi $SV.tbi \
  && echo "SV callset $(du -h $W/sv.vcf.gz | cut -f1) at $(date -u +%T)" ) > $W/sv.log 2>&1 &
SV_PID=$!
( bash scripts/fetch_data.sh p6 ) > $W/fasta.log 2>&1 &
FASTA_PID=$!

echo "=== data $(date -u +%T) ==="
GSA=https://webdata.illumina.com/downloads/productfiles/global-screening-array/v1-0
mkdir -p $W/gsa && cd $W/gsa
for z in infinium-global-screening-array-24-v1-0-c2-manifest-file-bpm-build38.zip \
         infinium-global-screening-array-v1-0-c1-cluster-file.zip; do
  curl -sSfL --retry 5 -O "$GSA/$z" && unzip -oq $z
done
cd - > /dev/null
BPM=$W/gsa/GSA-24v1-0_C2.bpm; EGT=$W/gsa/GSA-24v1-0_C1_ClusterFile.egt
.venv/bin/python scripts/p8/prepare.py --tsv scripts/p8/gse165845_samples.tsv --dl $W/geo --out $W/idats --workers 24 \
  | tail -2 || exit 1

echo "=== genotype $(date -u +%T) ==="
.venv/bin/simurg genotype --bpm $BPM --egt $EGT --idats $W/idats --out $W/gt --threads $N --quiet
cp $W/idats/samples_map.tsv $W/gt/sample_summary.tsv $R/

echo "=== model $(date -u +%T) ==="
wait $FASTA_PID; tail -1 $W/fasta.log
/usr/bin/time -v .venv/bin/simurg cn build-model --bpm $BPM --fasta data/p6/GRCh38.primary.fa --egt $EGT \
  --reference-gtcs $W/gt --max-reference 120 --out $W/gsa1.simurg-cyto 2> $R/model_time.txt | grep -v "GC windows"
grep -E "Elapsed|Maximum resident" $R/model_time.txt

echo "=== cn call $(date -u +%T) ==="
/usr/bin/time -v .venv/bin/simurg cn call --bpm $BPM --model $W/gsa1.simurg-cyto --gtcs $W/gt \
  --fasta data/p6/GRCh38.primary.fa --out $W/cn --threads $N --quiet 2> $R/call_time.txt
grep -E "Elapsed|Maximum resident" $R/call_time.txt
.venv/bin/simurg cn annotate --vcfs $W/cn --quiet || echo "annotate failed"
cp $W/cn/cn_sample_summary.tsv $W/cn/cn_call_manifest.json $R/

echo "=== truth $(date -u +%T) ==="
wait $SV_PID; cat $W/sv.log
cut -f3 $W/idats/samples_map.tsv | tail -n +2 | sort -u > $W/inds.txt
bcftools view -S $W/inds.txt --force-samples -i 'INFO/SVTYPE="DEL" || INFO/SVTYPE="DUP"' -Ou $W/sv.vcf.gz \
  | bcftools query -i 'GT="alt"' -f '[%CHROM\t%POS\t%INFO/END\t%INFO/SVTYPE\t%SAMPLE\n]' > $W/truth.tsv
wc -l < $W/truth.tsv
.venv/bin/python scripts/v2/cn_validate_truth.py --calls $W/cn --map $W/idats/samples_map.tsv --truth $W/truth.tsv \
  --bpm $BPM --out $R/validation
gzip -c $W/truth.tsv > $R/truth.tsv.gz
echo "=== task done $(date -u +%T) ==="