#!/usr/bin/env bash
# Build a training table from ClinVar, shaped like VEP --tab output so rarelens_ml.train reads it
# unchanged.
#
# ClinVar already carries the molecular consequence (MC), the gene (GENEINFO) and an allele
# frequency (AF_EXAC), which is the feature set serving sends. That avoids running VEP over
# hundreds of thousands of variants to produce a training set. CADD and AlphaMissense are left
# missing, exactly as they are when the pipeline runs without plugin data.
#
# Only 2-star-and-above records are kept: "criteria provided, multiple submitters, no conflicts"
# or better. Labels come from CLNSIG (see docs/data.md on evaluating this honestly).
set -eu

IMAGE=${BCFTOOLS_IMAGE:-quay.io/biocontainers/bcftools:1.20--h8b25389_0}
CLINVAR=${CLINVAR:-https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz}
REGION=${REGION:-}                      # empty means the whole genome
OUT=${OUT:-data/clinvar-training.vep.tsv}

mkdir -p "$(dirname "$OUT")"
docker run --rm -v "$PWD/$(dirname "$OUT"):/out" "$IMAGE" bash -eu -c "
  echo '==> streaming ClinVar ${REGION:-(whole genome)}' >&2
  bcftools query ${REGION:+-r '$REGION'} \
    -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/MC\t%INFO/GENEINFO\t%INFO/AF_EXAC\t%INFO/CLNSIG\t%INFO/CLNREVSTAT\n' \
    '$CLINVAR' \
  | awk -F'\t' -v OFS='\t' '
      BEGIN {
        split(\"transcript_ablation splice_acceptor_variant splice_donor_variant stop_gained frameshift_variant stop_lost start_lost transcript_amplification\", h, \" \");
        for (i in h) impact[h[i]] = \"HIGH\";
        split(\"inframe_insertion inframe_deletion missense_variant protein_altering_variant\", m, \" \");
        for (i in m) impact[m[i]] = \"MODERATE\";
        split(\"splice_region_variant synonymous_variant start_retained_variant stop_retained_variant\", l, \" \");
        for (i in l) impact[l[i]] = \"LOW\";
        print \"#Uploaded_variation\", \"Location\", \"Allele\", \"Consequence\", \"IMPACT\", \"SYMBOL\", \"gnomADe_AF\", \"CLIN_SIG\", \"CADD_PHRED\", \"am_pathogenicity\";
      }
      \$9 !~ /multiple_submitters|expert_panel|practice_guideline/ { next }   # 2 stars and up
      \$5 == \".\" || \$6 == \".\" { next }
      {
        split(\$5, mc, \",\"); split(mc[1], so, \"|\"); csq = so[2];
        if (csq == \"\") next;
        split(\$6, gi, \"|\"); split(gi[1], g, \":\"); gene = g[1];
        af = (\$7 == \".\" ? \"-\" : \$7);
        imp = (csq in impact ? impact[csq] : \"MODIFIER\");
        print \$1 \"_\" \$2 \"_\" \$3 \"_\" \$4, \$1 \":\" \$2, \$4, csq, imp, gene, af, \$8, \"-\", \"-\";
      }' > /out/$(basename "$OUT")
"
echo "wrote $OUT: $(($(wc -l < "$OUT") - 1)) labelled variants"
awk -F'\t' 'NR>1 {print $8}' "$OUT" | sort | uniq -c | sort -rn | head -6