#!/usr/bin/env bash
# Build a SIMULATED proband VCF: real GIAB HG002 variants as background, plus one real ClinVar
# pathogenic variant in a disease gene, which becomes the diagnosis the triage should find.
#
# Spiking a known variant into a public genome is how phenotype-driven triage tools are
# benchmarked. Every input here is public and openly licensed, and this is not a real patient.
# Provenance and citations: docs/data.md
set -euo pipefail

IMAGE=${BCFTOOLS_IMAGE:-quay.io/biocontainers/bcftools:1.20--h8b25389_0}
GENE=${GENE:-NF2}
REGION=${REGION:-chr22:20000000-31000000}
BACKGROUND=${BACKGROUND:-12}          # kept small: VEP's database mode is slow per variant
OUT_DIR=${OUT_DIR:-data}
CLINVAR=${CLINVAR:-https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz}
GIAB=${GIAB:-https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NISTv4.2.1/GRCh38/HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz}

mkdir -p "$OUT_DIR"
docker run --rm -v "$PWD/$OUT_DIR:/out" "$IMAGE" bash -eu -c "
  echo '==> background: GIAB HG002 $REGION' >&2
  bcftools view -H -r '$REGION' '$GIAB' \
    | awk '\$5 !~ /,/ && length(\$4) < 20 && length(\$5) < 20' \
    | awk 'NR % 149 == 0' \
    | head -n $BACKGROUND \
    | awk -v OFS='\t' '{ sub(/^chr/, \"\", \$1); print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \
    > /out/_background.tsv

  echo '==> diagnosis: a ClinVar pathogenic variant in $GENE' >&2
  bcftools view -H -r 22 -i 'INFO/CLNSIG ~ \"Pathogenic\" && INFO/GENEINFO ~ \"${GENE}:\" && INFO/CLNREVSTAT ~ \"multiple_submitters\"' '$CLINVAR' \
    | awk 'length(\$4) < 20 && length(\$5) < 20' \
    | head -n 1 \
    | awk -v OFS='\t' '{ print \$1, \$2, \".\", \$4, \$5, \".\", \"PASS\", \".\" }' \
    > /out/_spike.tsv

  test -s /out/_background.tsv || { echo 'no background variants found' >&2; exit 1; }
  test -s /out/_spike.tsv || { echo 'no 2-star pathogenic ClinVar variant found for $GENE' >&2; exit 1; }

  {
    printf '##fileformat=VCFv4.2\n'
    printf '##source=rarelens SIMULATED proband: GIAB HG002 background + one ClinVar pathogenic %s variant\n' '$GENE'
    printf '##contig=<ID=22,length=50818468>\n'
    printf '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n'
    sort -k2,2n /out/_background.tsv /out/_spike.tsv
  } | bgzip > /out/proband-simulated.vcf.gz
  tabix -f -p vcf /out/proband-simulated.vcf.gz
  rm -f /out/_background.tsv
"
echo
echo "wrote $OUT_DIR/proband-simulated.vcf.gz"
echo "the planted diagnosis (chrom pos ref alt):"
awk '{print "  " $1, $2, $4, $5}' "$OUT_DIR/_spike.tsv"
rm -f "$OUT_DIR/_spike.tsv"