Kemal Yaylali
197975cc42
"Variants are unscored" was accurate: nothing was ever trained, so a quarter of every rank was dead weight and the UI leaked a connection error at the reader. - scripts/make-training-set.sh derives a training table from ClinVar directly. ClinVar already carries the molecular consequence, the gene and an allele frequency, which is the feature set serving sends, so this avoids running VEP over hundreds of thousands of variants. 2-star records only. - train.py now holds out whole genes (GroupShuffleSplit). docs/data.md had said to do this since the data pass; the code was still doing a random split, which is the leak Grimm 2015 describes. - evaluate() reports missense on its own. On the last run: AUROC 0.986 over 74,239 held-out variants, but 0.872 over the 13,553 missense ones, and the docs say plainly why even that is flattered — within missense the only live feature is allele frequency, and ClinVar's benign calls often use allele frequency as evidence (ACMG BA1/BS1), so the feature partly caused the label. - the 503 now names what is missing (model@alias via tracking URI) and leaves the exception in the server log instead of the UI. - make training-set / make train; the 58 MB table is gitignored. Verified end to end: model registered as v2, the simulated NF2 case scores 0.999 on the planted variant, and it now ranks 1.00 with all four components live. Tests: api 77, ml 22, loader 16, web 32; ruff, mypy, svelte-check clean.
48 lines
2.7 KiB
Bash
Executable File
48 lines
2.7 KiB
Bash
Executable File
#!/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
|