bio-causal-genomics-transcriptome-wide-association

Version Compatibility

Reference examples tested with: FUSION (head of gusevlab/fusion_twas, scripts dated 2023+), MetaXcan / S-PrediXcan / S-MultiXcan 0.7.5+ (hakyimlab/MetaXcan), PrediXcan model files from PredictDB (GTEx v8 elastic-net + MASHR), UTMOST (head of Joker-Jerome/UTMOST), pyfocus 0.8+ (bogdanlab/focus), MA-FOCUS (head of mancusolab/ma-focus), TIGAR-V2 (head of yanglab-emory/TIGAR), PLINK 1.9 + PLINK 2.0, R 4.3+, Python 3.9-3.11.

Before using code patterns, verify installed versions match. If versions differ:

• R: Rscript --version; for FUSION scripts inspect --help flags directly in the source
• Python: pip show pyfocus (MetaXcan is git-cloned, not on PyPI) then SPrediXcan.py --help, SMulTiXcan.py --help, focus finemap --help
• CLI: plink2 --version; FUSION ships as R scripts not a binary

If a script throws an error about an argument that has moved (e.g. --gwas_file vs --gwas-file) or a model database schema change, introspect the installed script with --help and adapt rather than retrying. PredictDB model file paths change with GTEx version; pin the version explicitly in scripts.

Transcriptome-Wide Association

"Find genes whose predicted tissue expression is associated with my GWAS trait" -> Train SNP -> expression prediction models on a reference eQTL panel, apply the per-gene SNP weights to GWAS summary statistics or genotypes, and produce a gene-level Z-score equivalent to a weighted sum of SNP Z-scores. The output is a gene-by-tissue association, but TWAS is NOT direct evidence of causal mediation: an LD-tagged eQTL signal produces the same statistical association as a truly causal one, and the dominant failure modes are LD-induced false positives at gene-dense loci, tissue mis-specification, and ancestry mismatch between GWAS and prediction weights.

• CLI (sumstat TWAS, R): FUSION.assoc_test.R --sumstats g.sumstats --weights weights.pos --weights_dir wgt/ --ref_ld_chr 1KG/EUR. --chr 22 --out chr22.dat
• CLI (S-PrediXcan, Python): SPrediXcan.py --model_db_path gtex_v8.db --covariance gtex_v8.cov --gwas_file g.txt --output_file out.csv
• CLI (S-MultiXcan joint): SMulTiXcan.py --models_folder mashr_models/ --gwas_folder gwas/ --metaxcan_folder spredixcan_per_tissue/ --output joint.csv
• CLI (UTMOST cross-tissue): joint test across tissues via UTMOST's per-tissue GBJ / GBJ2 step
• CLI (FOCUS fine-mapping): focus finemap gwas.sumstats 1KG_EUR focus.db --chr 22 --p-threshold 5e-8 --out chr22.focus
• CLI (MA-FOCUS multi-ancestry): focus finemap with colon-separated per-ancestry sumstats / LD / weights and hyphen-joined ancestry codes in --locations (e.g. 38:EUR-EAS-AFR)

TWAS, cis-eQTL MR, and coloc operate on overlapping evidence: TWAS asks "is the gene's predicted expression associated with the trait?"; cis-eQTL MR asks "does the eQTL effect on expression mediate the trait effect under IV assumptions?"; coloc asks "do the GWAS and eQTL share a causal variant?". Strong causal claims require triangulation, not single-method significance.

Algorithmic Taxonomy

Tool	Model	Input	Output	Strength	Fails when
FUSION (Gusev 2016 Nat Genet 48:245)	Weighted sum of SNP Z-scores; per-gene multi-model (BLUP, lasso, elnet, top1, bslmm) selected by cross-validation R^2	GWAS sumstats + pre-computed weights (.pos + per-gene RData)	TWAS Z, p; conditional joint analysis	Mature, ENCODE-style, pre-trained weights for many tissues (GTEx, CMC, METSIM, YFS, NTR)	LD-induced false positives at gene-dense loci; needs ancestry-matched LD reference; weights are heritability-thresholded so low-h2 genes drop
PrediXcan (Gamazon 2015 Nat Genet 47:1091)	Elastic-net SNP -> expression prediction; individual-level genotypes	PLINK genotypes + GWAS phenotype + GTEx prediction DB	Per-gene association test with full regression machinery	Most flexible (allows covariates, interactions, binary outcomes); same TWAS interpretation	Requires individual-level data; biobank-scale compute
S-PrediXcan (Barbeira 2018 Nat Commun 9:1825)	Summary-statistic equivalent of PrediXcan; Z-score weighted sum analogous to FUSION	GWAS sumstats + PredictDB model + covariance	Per-gene Z, p	Public PredictDB models (GTEx v8 elastic-net + MASHR-EUR); minimal compute	Pre-trained weights ancestry-specific (EUR primarily); covariance file must match model DB
S-MultiXcan (Barbeira 2019 PLoS Genet 15:e1007889)	Joint multi-tissue test combining per-tissue S-PrediXcan via PCA-regularised regression	Folder of per-tissue S-PrediXcan outputs + model folder	Single joint p per gene + per-tissue significance	Boosts power when causal tissue is unknown; standard for transcriptome-wide screens	Joint test cannot pinpoint causal tissue; correlated tissues produce ill-conditioned regression
UTMOST (Hu 2019 Nat Genet 51:568)	Cross-tissue elastic-net (group lasso) for joint tissue prediction + GBJ test	Per-tissue eQTL data + GWAS sumstats	Cross-tissue joint statistic	Often better-powered than S-MultiXcan at cross-tissue genes	Computationally heavier than MetaXcan; tissue weights are less interpretable
MOSTWAS (Bhattacharya 2021 PLoS Genet 17:e1009398)	Mediator-aware TWAS adding distal trans-mediating SNPs to cis-only models	GWAS sumstats + MOSTWAS weights	Per-gene Z, p (TWAS + distal mediator extension)	Recovers signal at genes with non-cis genetic regulation	Trans-mediation models need large reference panel; weights less broadly available
kTWAS (Cao 2021 Brief Bioinform 22:bbaa270)	Kernel-based TWAS using SKAT-style aggregation	GWAS sumstats or genotypes + per-gene SNP set	Per-gene p	Robust to non-linear and rare-variant contributions	Loses the eQTL-weighting interpretability; less standardised
EpiXcan (Zhang 2019 Nat Commun 10:3834)	Adds epigenome-derived per-SNP prior to weight training	eQTL + epigenome + GWAS sumstats	Per-gene Z, p (epigenome-informed)	Higher prediction R^2 in epigenome-rich tissues	Requires matched epigenome data for the prediction tissue
TIGAR-V2 (Parrish 2022 HGG Adv 3:100068)	Dirichlet process regression (non-parametric Bayes) for SNP -> expression	Reference eQTL + GWAS sumstats	Per-gene Z, p	Captures non-elastic-net effect structures; more flexible weight learning	Slower training; benefits depend on locus genetics
FOCUS (Mancuso 2019 Nat Genet 51:675)	Probabilistic gene-level fine-mapping over TWAS Z-scores using gene-by-gene predicted-expression correlation as analog of LD	FUSION/S-PrediXcan TWAS Z + ancestry-matched LD reference	Per-gene PIP + credible gene set	Resolves co-regulated gene clusters into a probabilistic causal gene; standard add-on after TWAS	Requires the same prediction-weight panel that produced TWAS Z; PIPs depend on prior
MA-FOCUS (Lu 2022 AJHG 109:1388-1404)	Multi-ancestry FOCUS; joint gene fine-mapping across ancestries with shared causal-gene assumption	Per-ancestry TWAS sumstats + per-ancestry weights + per-ancestry LD	Cross-ancestry PIP	Smaller credible gene sets when AFR/EAS contribute non-EUR LD information	Trans-ethnic gene-effect heterogeneity violated; weights must be ancestry-matched
JEPEG / JEPEG-Mix (Lee 2015 / 2016)	Gene-based test combining eQTL and functional weights	GWAS sumstats + JEPEG annotation database	Per-gene p	Lightweight gene-burden alternative to TWAS	Less granular than full PrediXcan/FUSION machinery; minimally updated

Methodology evolves; the FUSION-vs-PrediXcan landscape has been stable but probabilistic fine-mapping (FOCUS, MA-FOCUS) and integrative methods (MOSTWAS, EpiXcan, OmicsXcan) continue to advance. Verify against the current PredictDB release notes (predictdb.org) and the latest FUSION weight panels (gusevlab.org/projects/fusion) before locking on a tissue or model.

S-PrediXcan and FUSION are mathematically near-identical: both compute a weighted sum of GWAS SNP Z-scores using per-gene SNP-expression weights and an LD-aware variance correction. The practical differences reduce to (a) the weight panel (FUSION elastic-net vs PredictDB MASHR), (b) the LD reference, and (c) the per-gene heritability threshold; the supplement of Barbeira 2018 works through the algebra. Method choice should therefore be driven by panel availability and ancestry-match, not by the underlying algorithm.

PredictDB Model Choice and GTEx Versioning

Release	Cohort	Status (2026)	When to use
GTEx v8 (838 donors, 49 tissues, 2020)	EUR-dominant (~85%)	Current PredictDB standard	Default; pre-trained MASHR + elastic-net DBs available
GTEx v9 (2023)	Expanded harmonization	Not migrated into PredictDB	Do not use until PredictDB rebuilds
GTEx v10 (2024 AnVIL release)	Re-aligned to GRCh38 v44	Limited harmonization; not PredictDB-default	Wait for community-validated weight panels

Operational rule: Use GTEx v8 unless there is an explicit biological reason to deviate (tissue not in v8, ancestry-specific panel preferred). In methods, pin exactly: "GTEx v8 MASHR-EUR, PredictDB release 2022-01".

MASHR vs Elastic-Net Models

PredictDB ships two cross-validated model families per tissue (Barbeira 2021 Genome Biol 22:49):

Model	Construction	Per-gene SNP count	When to use
MASHR	Cross-tissue posterior mean from DAP-G fine-mapped SNPs	~10x sparser	Primary discovery; higher per-gene R^2 in most genes; standard for S-MultiXcan
Elastic-net	Per-tissue lasso/ridge mix (alpha = 0.5)	Denser	Tissues where MASHR's cross-tissue prior is mis-specified (ovary, testis, isolated-organ traits)

Operational rule: Never mix MASHR and elastic-net within a single S-MultiXcan run; the inter-tissue covariance and condition number assumptions break. Pick one family and apply consistently across all tissues.

Decision Tree by Experimental Scenario

Scenario	Recommended workflow	Why
GWAS summary stats only, EUR, single hypothesis tissue (e.g. liver for LDL)	S-PrediXcan with GTEx v8 MASHR-EUR weights, or FUSION with GTEx liver	Standard pre-trained pipeline; minimal compute
GWAS summary stats only, tissue unknown a priori	S-MultiXcan (standard GTEx v8) OR UTMOST (custom panel) -- see S-MultiXcan vs UTMOST table	Joint multi-tissue inflates power; tissue prioritisation requires LDSC-SEG separately
Multiple TWAS hits at one locus (gene-dense region)	Run TWAS then FOCUS for probabilistic fine-mapping	LD ties co-regulated genes; FOCUS PIP distinguishes likely causal gene
Multi-ancestry GWAS (EUR + EAS + AFR)	Per-ancestry S-PrediXcan with matched weights, then MA-FOCUS to combine	Single-ancestry weights miscalibrated in other ancestries; joint fine-mapping shrinks credible set
Individual-level genotypes available (UKB)	PrediXcan (full regression)	Allows covariates, interactions, binary outcomes natively
Low-N tissue (GTEx N < 100)	Substitute eQTLGen (whole blood, N ~ 31k) OR skip the tissue	Per-gene CV R^2 unstable below N ~ 100; weights overfit
Drug-target prioritisation (TWAS as causal-gene evidence)	TWAS + cis-eQTL MR + coloc + FOCUS triangulation	TWAS alone is associational; triangulation strengthens causal claim
Trans-acting / mediator-aware analysis	MOSTWAS	Adds distal trans-mediating SNPs to cis-only models
Rare-variant or non-linear gene effects	kTWAS or TIGAR-V2	Kernel / non-parametric flexibility
HLA region (chr6:25-35 Mb hg38)	Exclude or use HLA-specific tools	Long-range LD breaks every gene-by-gene method; standard TWAS PIPs not interpretable
Splicing-mediated trait (e.g. neuropsych for sQTL)	sTWAS (sQTL-weighted TWAS) using GTEx splicing models	Splicing mediates many GWAS effects; cis-sQTL panels available in PredictDB
Cell-type-specific trait	sc-eQTL-based TWAS (e.g. OneK1K, Yazar 2022 Science 376:eabf3041)	Bulk-tissue TWAS averages over cell types; single-cell eQTL recovers cell-type specificity

Tissue Selection Protocol

Tissue choice drives TWAS power and false-positive rate; selecting tissues by inspecting TWAS hit count is circular. Run all three of the following on the GWAS sumstats (independent of any TWAS run) and pick the primary TWAS tissue from the intersection:

1. Stratified LDSC tissue prioritization (Finucane 2018 Nat Genet 50:621): ldsc.py --h2-cts <sumstats> --ref-ld-chr-cts <annot> --w-ld-chr <weights> against 200+ tissue-specific gene expression annotations
2. CELLEX (Timshel 2020 eLife 9:e55851): single-cell tissue / cell-type prioritization on the same GWAS
3. MAGMA gene-property analysis (de Leeuw 2015 PLoS Comput Biol 11:e1004219): cheaper substitute when LDSC unavailable

Operational rule: Primary TWAS tissue = the tissue with FDR-significant enrichment in at least two of the three methods. Run secondary tissues in S-MultiXcan for cross-tissue replication. Bonferroni for tissue selection alone: 0.05 / 200 annotations = 2.5e-4.

S-MultiXcan vs UTMOST

Method	Use case	Rationale
S-MultiXcan (Barbeira 2019)	Standard GTEx v8 analysis	Pre-computed MASHR weights; lower compute barrier; PCA-regularised inter-tissue regression
UTMOST (Hu 2019 Nat Genet 51:568)	Custom eQTL panel with cross-tissue retraining	Higher power at genes with shared cross-tissue eQTL architecture; group-lasso enforces sparsity across tissues

Benchmarks: Hu 2019, Barbeira 2019. Choose by panel availability first; the methods recover overlapping but non-identical gene sets.

Single-Cell and Cell-Type-Resolved TWAS

Bulk-tissue TWAS averages over cell composition; sc-eQTL TWAS recovers cell-type-specific regulation but at lower per-cell-type power.

Panel	Reference	Cells / tissue
OneK1K	Yazar 2022 Science 376:eabf3041	PBMC, ~982 donors, 14 cell types
HipSci iPSC-eQTL	Kilpinen 2017 Nature 546:370	iPSC, ~317 donors
BLUEPRINT	Chen 2016 Cell 167:1398	Monocytes, neutrophils, T cells

Tooling state (2026): No fully pre-built scPrediXcan equivalent to MASHR; sc-eQTL weights are panel-specific. Train custom PredictDB or use TIGAR-V2's Bayesian DPR on the sc-eQTL matrix.

Operational rule: Run standard bulk-tissue TWAS first; run sc-eQTL TWAS in the prioritized cell type as a secondary analysis; require concordance between bulk and sc results before nominating a cell-type-specific gene. Upstream sc preprocessing: cross-reference single-cell/preprocessing.

Per-Tool Failure Modes

LD-induced TWAS false positives (most common pitfall)

Trigger: Two or more genes at the same locus have correlated cis-eQTLs (shared causal eQTL SNP or LD-linked eQTL SNPs).

Mechanism: TWAS Z-scores are linear combinations of SNP Z-scores weighted by per-gene SNP effects. When two genes share many high-weight SNPs (e.g. nearby genes regulated by the same enhancer or LD-tagged independent eQTLs), their TWAS Z-scores are positively correlated. A single causal GWAS variant therefore produces significant Z at multiple co-regulated genes (Wainberg 2019 Nat Genet 51:592; Mancuso 2019 Nat Genet 51:675).

Symptom: A GWAS lead locus shows 3-10 genes all passing genome-wide TWAS significance (p < 2.3e-6 ~ 0.05/22k); per-gene LocusZoom-style plots look near-identical; conditional analysis (FUSION.post_process.R runs conditional/joint analysis by default; FUSION.assoc_test.R's --coloc_P adds single-SNP coloc) reveals only 1-2 independent gene signals; the genes lie within 1 Mb of each other.

Fix: Always run FOCUS after TWAS to obtain per-gene PIPs. Report only genes with PIP >= 0.8 as candidate causal; report co-significant genes with PIP < 0.5 as LD-tagged. Cross-check with cis-eQTL coloc (PP.H4 >= 0.7) for the candidate causal gene. FUSION's FUSION.post_process.R performs conditional/joint analysis by default (no --joint flag) as a lighter-weight alternative.

Tissue mis-specification

Trigger: Running TWAS in a tissue that does not host the causal regulatory effect (e.g. whole blood for a psychiatric trait; pancreas for an LDL trait).

Mechanism: A gene's cis-eQTL effect size varies across tissues; a wrong-tissue model has weaker per-gene prediction R^2 and lower power. Conversely, eQTL effects in the wrong tissue can still tag the GWAS signal via LD and produce spurious associations not present in the causal tissue.

Symptom: Strong TWAS signal in a tissue biologically irrelevant to the trait; null in the expected tissue; tissue-prioritisation methods (LDSC-SEG, Finucane 2018 Nat Genet 50:621; RolyPoly (Calderon 2017 AJHG 101:686)) disagree with the TWAS tissue.

Fix: Run S-MultiXcan to combine tissues if causal tissue is unknown. For prioritisation, use LDSC-SEG / CELLEX / EWCE on the GWAS sumstats independently of TWAS, and report TWAS in the prioritised tissues. Never report a single-tissue TWAS hit as causal without independent tissue evidence (single-cell eQTL, chromatin accessibility in matched cell type).

Ancestry mismatch in prediction weights

Trigger: Running TWAS on a non-EUR GWAS using GTEx (~ 85% EUR) weights, or vice versa.

Mechanism: Cis-eQTL effect sizes and LD structure are ancestry-specific; prediction weights trained in one ancestry transfer with reduced R^2 and biased Z-scores in another. Power is lost preferentially at loci where the causal eQTL is not shared across ancestries (Patel 2022 AJHG 109:1286).

Symptom: Genome-wide TWAS hit count much lower than expected given GWAS power; non-EUR-specific GWAS loci fail to produce TWAS hits; per-gene prediction R^2 substantially reduced.

Fix: Use ancestry-matched prediction panels where available: MESA multi-ethnic eQTL (Mogil 2018 PLoS Genet), eQTLGen-Asian, AFGR (Africa) when published, or MAGE (Taliun-style multi-ancestry eQTL). Move to MA-FOCUS for cross-ancestry joint fine-mapping. Document the ancestry assumption explicitly in methods.

Low-N tissue weights are unstable

Trigger: Using a GTEx tissue with N < 100 donors (e.g. several brain sub-regions, kidney cortex in v7).

Mechanism: Per-gene elastic-net weights are cross-validated with the available donors. Below ~ 100 donors, the cross-validation R^2 has high variance and the heritability filter (FUSION requires hsq_p < 0.01) drops many genes. Surviving weights overfit, inflating per-gene Z under the null.

Symptom: Tissue produces unusually high TWAS hit count or unusually high genomic inflation; per-gene CV R^2 distribution is bimodal with a long heavy tail.

Fix: Skip GTEx tissues with N < 100 unless biologically essential. Substitute eQTLGen for whole blood (N ~ 31k, Vosa 2021 Nat Genet 53:1300) where blood is acceptable. For brain, use PsychENCODE (N ~ 1300, Wang 2018 Science 362:eaat8464) or BrainSeq (N ~ 350+) when available; verify the matching prediction-weight panel exists.

HLA region

Trigger: Any gene within chr6:25-35 Mb (hg38; extended MHC) reported by TWAS.

Mechanism: Long-range LD (r2 > 0.5 over many Mb) and extreme structural variation mean per-gene prediction weights at HLA capture haplotype rather than gene-specific regulation. Standard TWAS gene-level inference is biologically meaningless here.

Fix: Exclude chr6:25-35 Mb from genome-wide TWAS summaries by default. For HLA-driven traits (autoimmune, infection, transplantation), impute classical HLA alleles using one of:

Tool	Reference	Notes
HIBAG	Zheng 2014 Pharmacogenomics J 14:192	R package; pre-trained per-ancestry classifiers
SNP2HLA	Jia 2013 PLoS One 8:e64683	Beagle-based imputation; supports T1DGC reference
HLA-TAPAS	Luo 2021 Nat Genet 53:1504	Current standard; multi-ancestry reference; recommended for new analyses

Then test classical alleles plus amino-acid residues (Raychaudhuri 2012 Nat Genet 44:291 set the gold standard for residue-level association in MHC). Do NOT run SNP-level TWAS inside the MHC.

Correlated-expression confounding (co-regulated genes)

Trigger: Two or more genes are functionally co-regulated by a single TF or enhancer, producing nearly-identical predicted-expression vectors.

Mechanism: Even with separate per-gene cis-eQTL prediction, downstream co-regulation makes predicted expression highly correlated; TWAS cannot distinguish which gene mediates the trait.

Symptom: FOCUS credible gene set contains multiple genes with PIP roughly equal (e.g. three genes at 0.3 each); functional follow-up (MPRA, CRISPRi screens) is needed to break the tie.

Fix: Acknowledge the limit of statistical resolution; report the full credible gene set and prioritise on orthogonal evidence (CRISPRi/CRISPRa effect size in matched cell type, e.g. Open Targets-style, or MPRA at allelic series; protein-level pQTL coloc if available).

TWAS - MR - Coloc Triangulation

A TWAS-significant gene is associational, not causal. The strongest defensible claim that a gene mediates a GWAS effect comes from triangulating three orthogonal lines of evidence:

Line	What it tests	Threshold
TWAS	Predicted-expression association with trait	S-MultiXcan joint p < 2.3e-6, OR S-PrediXcan per-tissue p < 4.6e-8 (49-tissue Bonferroni), OR per-tissue FDR < 0.05
cis-eQTL MR	Causal effect of expression on trait under IV assumptions (cross-reference causal-genomics/mendelian-randomization)	Wald-ratio or IVW p < 0.05/n_genes; instrument F > 10
Colocalization	Shared causal variant between GWAS and eQTL (cross-reference causal-genomics/colocalization-analysis)	coloc.abf or coloc.susie PP.H4 >= 0.7
FOCUS	Probabilistic per-gene PIP under TWAS fine-mapping	PIP >= 0.8

Operational rule: Report a gene as a "strong candidate causal gene" only when 3 of the 4 are concordant (TWAS hit + coloc PP.H4 >= 0.7 + FOCUS PIP >= 0.8, with cis-MR as a supporting fourth). 2-of-4 concordance is "suggestive"; 1-of-4 is "associational only". The combination is more conservative than any single method but matches the standards used in modern GWAS-to-target pipelines (Open Targets Genetics Mountjoy 2021 Nat Genet 53:1527; FinnGen R10 release notes).

Quantitative Thresholds