Diagnose and resolve SNP/variant-id mismatch in S-PrediXcan and Predict.py

[hakyim]

Summary

SNP/variant-id mismatch between the model DB and the GWAS/genotype input (rsid vs varID, hg19 vs hg38, non-rsID GWAS ids, a covariance keyed differently from the model) is the single largest recurring source of user-reported failures against this repo. Matching is a plain exact-string merge that silently succeeds even at ~0% overlap, producing a near-empty or all-NA results file with only an INFO-level log line as a clue.

This PR does two things: it makes that failure loud, and then it stops most of it from happening.

Diagnostics — check_snp_overlap() compares the model's SNP set against the GWAS's before any per-gene work, and check_covariance_overlap() compares the model's against the covariance's. Each logs a WARNING naming the overlap %, example ids from both sides, and the specific flags to try. Predict.py (individual-level genotype) has the same silent-failure shape; since it streams genotypes rather than materializing a variant-id set, it reuses the id-match count already tracked during the prediction loop and warns post-loop, with genotype-specific advice (--variant_mapping, --on_the_fly_mapping, --skip_palindromic, --liftover).

Resolution — metax/misc/SnpKeyResolution.py samples the ids on all three sides before anything is loaded and lines them up: it picks which of the model db's id columns to match on, translates the GWAS's ids to it when needed, and lifts the loader's rsid-only filter when the GWAS isn't rsID-keyed. It says what it did, at WARNING.

Why the diagnostic keys on id format, not a percentage

The first draft fired on a low overlap percentage alone. Testing against real data showed that was wrong in both directions:

invocation	model↔GWAS	model↔cov	usable genes (before)	after
--snp_column panel_variant_id --model_db_snp_key varID --keep_non_rsid (correct)	100%	100%	15/15	15/15
--snp_column panel_variant_id --keep_non_rsid (classic mistake)	8.11%	8.11%	3/15	15/15
--snp_column panel_variant_id --model_db_snp_key varID	0%	100%	0	15/15
--snp_column variant_id	91.89%	8.11%	0 (all NA)	15/15

1. A 1% floor misses the canonical MASHR mistake (row 2, at 8.11%): a MASHR db's rsid column falls back to the varID string for variants that have no rsID, and those few match.
2. Raising the floor is not the fix — a GWAS covering a single chromosome legitimately overlaps a genome-wide model by a few percent, so any floor high enough to catch row 2 false-alarms on that.

So the primary signal is the variant id format differing between the two sides (metax/misc/SnpOverlapDiagnostics.py), which holds regardless of how much of the genome the GWAS covers. min_pct=1.0 is kept only as a floor for what formats can't see — an hg19/hg38 mismatch, where both sides are chr_pos_ref_alt but no position ever lines up. Two unrecognized formats are never called a mismatch, since a false alarm is worse than a miss here.

Row 4 is a failure the GWAS-overlap check structurally cannot see: model↔GWAS matching is fine, but PredictDB's MASHR covariances are keyed by varID, so an rsID-keyed model finds nothing in the covariance and writes an all-NA file. check_covariance_overlap applies the same format-aware rule one layer over, backed by a new MatrixManager.snps(). StreamedMatrixManager.snps() returns None and the check is skipped, since a streamed covariance never has the full id set in hand.

How the resolution works, and why it constructs nothing

A PredictDB model db carries both rsid and varID for the same weights row — the same variant, with the same alleles. So the db is itself an exact id table, and no id has to be built out of chromosome/position/alleles. That matters: constructing a chr_pos_ref_alt key means inferring allele orientation, and a mistake there silently flips the sign of an association rather than failing loudly.

1. Pick the model's key column. Chosen to agree with the covariance when the covariance's format is known, otherwise with the GWAS — a model the covariance can't be looked up in yields all-NA no matter how well the GWAS matched.
2. Translate the GWAS's ids to that key through the db's own weights table, when the GWAS names its variants by a different db column (row 4). An id naming more than one variant in the target column is dropped rather than resolved arbitrarily.
3. Lift the rsid-only load filter when the GWAS's ids aren't rsIDs and would otherwise all be dropped at read time (row 3).

Downstream is untouched: the same exact-string merge, and align_data_to_alleles still doing the allele matching and sign correction on unchanged input.

Overrides: --model_db_snp_key is never overridden, only diagnosed. New --gwas_snp_key names the GWAS's source column by hand. --snp_map_file turns resolution off entirely, since ids are then already being mapped. Resolution runs on a copy of args, because MetaMany drives SPrediXcan.run in a loop over model dbs with one args and a key resolved for one db need not exist in the next.

Known limit: a variant the db has no rsID for can't be named on the GWAS's side of an rsID→varID lookup, so it's dropped and the count is logged. On the fixtures that's 3 of 37 — row 4 produces 15/15 genes from 34 variants rather than 37.

Test plan

• 150 tests pass (python -m unittest discover -s . -p "test_*.py" -t . from software/), except one pre-existing failure noted below.
• Real-data fixtures: software/tests/_td/qgt_chr1_subset/ (36 KB) holds a public slice of the QGT course data — a 15-gene chr1 subset of GTEx v8 MASHR Whole Blood plus its covariance, and the 37 matching variants of the imputed CARDIoGRAM_C4D_CAD GWAS. The GWAS carries both variant_id (rsID) and panel_variant_id (varID), so every row of the table is reproducible by changing only --snp_column.
• Allele-flip coverage, since a sign error is the failure mode that would be silent rather than empty: the fixture GWAS with every variant reported against the other allele (effect/non-effect swapped, zscore negated) must produce bit-identical results, both through the key switch and through the translation.
• Rows 2 and 3 now produce byte-identical output to the explicitly-correct row 1. For row 4, every gene that keeps all its snps matches row 1 exactly, and the rest differ only by having had fewer snps.
• Unit tests for the pure helpers: id classification and dominant format (test_snp_overlap_diagnostics.py), overlap checks (test_M04_zscores.py), key choice and translation (test_snp_key_resolution.py).
• Predict.py verified to import/compile cleanly; no genotype fixtures exist in tests/_td for a full Predict.run() integration test.
• Maintainer: run the full suite in the imlabtools conda env.

Two notes for whoever runs the suite, both pre-existing and unrelated to this PR:

• tests/test_M01_covariances_correlations.py::testProcessWeightDBRun fails on master. Commit a6278e5 ("Replace covariance column delimiter by tab instead of single space") changed M01_covariances_correlations.py without updating the test, which still expects space-delimited output. The covariance reader uses sep="\s+" so both parse, but any downstream consumer expecting single spaces would break — worth deciding whether the test or the format change is the thing to fix.
• The suite must be run as a package from software/; the -s tests form in CLAUDE.md fails on the relative imports. sqlalchemy is needed for the model-db tests and is missing from conda_env.yaml.

🤖 Generated with Claude Code