bio-causal-genomics-mendelian-randomization

## Version Compatibility Reference examples tested with: TwoSampleMR 0.6.0+, MendelianRandomization 0.10+, MR-PRESSO 1.0+, cause 1.2+, MVMR 0.4+, ieugwasr 1.0+, MRlap 0.0.3.2+, coloc 5.2+, mrclust 0.1+, lhcMR 0.0.1+, R 4.4+. Both TwoSampleMR 0.6.0 and ieugwasr 1.0 are the JWT-transition versions; older versions still expect deprecated OAuth. Before using code patterns, verify installed versions match. If versions differ: - R: `packageVersion('<pkg>')` then `?function_name` to verify parameters - CLI (plink, GCTA-GSMR): `<tool> --version` then `<tool> --help` If code throws an error referencing a function that has moved (e.g. `ieugwasr::ld_clump` vs `TwoSampleMR::clump_data`) or an OAuth token failure, introspect the installed API and adapt the example rather than retrying. # Mendelian Randomization **"Test whether trait X causally affects trait Y from GWAS summary statistics"** -> Use genetic variants as instrumental variables (IVs) that satisfy three assumptions (relevance, independence, exclusion restriction) to estimate `beta_causal = beta_outcome / beta_exposure` under the IV framework (Davey Smith & Ebrahim 2003 IJE 32:1; Burgess & Thompson 2021 Chapman & Hall/CRC, 2nd ed.). Tool choice is a decision about the **regime** (one-sample vs two-sample, sparse vs polygenic, drug-target vs polygenic exposure) and the **pleiotropy model** (balanced, directional InSIDE, correlated horizontal). Wrong tool inflates Type-I error or attenuates true effects in a direction predictable from the bias structure. - R: `TwoSampleMR::mr()` orchestrates IVW + Egger + weighted median + weighted mode in one call - R: `MendelianRandomization::mr_ivw / mr_egger / mr_median / mr_mbe / mr_conmix` per-method API (S4 objects; MR-RAPS is NOT in this package -- use `TwoSampleMR::mr_raps()` which wraps the GitHub `mr.raps`) - R: `MRPRESSO::mr_presso()` global / outlier / distortion tests - R: `cause::cause()` correlated horizontal pleiotropy mixture - R: `MVMR::strength_mvmr() + MVMR::ivw_mvmr()` multivariable conditional-F + IVW ## Statistical Model Taxonomy | Method | Pleiotropy assumption | Min instruments | Strength | Fails when | |--------|------------------------|-----------------|----------|------------| | IVW (fixed) | All IVs valid | 2 | Most efficient under no pleiotropy | Any directional or balanced pleiotropy inflates Type-I | | IVW (random effects) | Balanced + InSIDE | 3 | Standard primary; absorbs heterogeneity into wider SE | Directional pleiotropy biases the point estimate | | MR-Egger | Directional pleiotropy + InSIDE | 10+ for power | Detects + corrects directional pleiotropy via intercept (Bowden 2015 IJE 44:512) | NOME violated (`I^2_GX < 0.9`); SIMEX correction required; underpowered <10 SNPs | | Weighted median | Up to 50% invalid IVs | 3 | Robust to a minority of bad instruments (Bowden 2016 Genet Epidemiol 40:304) | >50% invalid IVs | | Weighted mode | Zero modal pleiotropy (ZEMPA) | 3 | Robust if the modal estimate is unbiased (Hartwig 2017 IJE 46:1985) | Bimodal pleiotropy; small numbers | | MR-RAPS | Balanced pleiotropy + weak instruments | 10+ | Profile-score robust to weak-IV + balanced horizontal pleiotropy (Zhao 2020 Ann Stat 48:1742) | Strong directional pleiotropy; CRAN-archived 2025-03-01 | | CAUSE | Correlated horizontal pleiotropy (CHP) | 100+ sig SNPs | Explicit shared-factor mixture; protects against CHP-driven false positives (Morrison 2020 Nat Genet 52:740) | Sparse polygenic exposures; <100 sig SNPs | | GSMR + HEIDI-outlier | Outlier removal under InSIDE | 10+ | Alternative outlier detection; integrates with LD reference (Zhu 2018 Nat Commun 9:224) | Requires individual-level LD; HEIDI conservative | | MR-PRESSO | Outlier-driven horizontal pleiotropy | 4+ | Global / outlier / distortion three-step (Verbanck 2018 Nat Genet 50:693) | Blind to CHP; computationally heavy at large NbDistribution | | MVMR (IVW) | Conditional independence after measured pleiotropy | 1+ per exposure | Accounts for measured horizontal pleiotropy via multivariable regression (Sanderson 2019 IJE 48:713) | Conditional F < 10 on any exposure | | MR-Clust | Heterogeneous causal effects (multiple mechanisms) | 30+ | Clusters SNPs by their causal-effect estimate (Foley 2021 Bioinformatics 37:531) | Single causal mechanism; small instrument sets | | Contamination mixture | Mixture of valid + invalid IVs | 10+ | Profile-likelihood mixture (Burgess 2020 Nat Commun 11:376) | Sparse signal | | LCV | Genome-wide; distinguishes causation vs genetic correlation | All SNPs | Tests `gcp` parameter using LDSC-style block jackknife (O'Connor & Price 2018 Nat Genet 50:1728) | Two-trait covariance dominated by a third confounder | | LHC-MR | Bidirectional + heritable confounder | All SNPs | Joint likelihood over genome-wide markers; estimates both directions + confounder (Darrous 2021 Nat Commun 12:7274) | Computationally heavy; rare-variant trait | | MRlap | Sample overlap + winner's curse + weak-IV jointly | Genome-wide sumstats | LDSC-scaffolded joint correction (Mounier & Kutalik 2023 Genet Epidemiol 47:314) | LDSC intercept poorly estimated (h^2 < 0.05); non-EUR without matched LD scores | | Doubly-Ranked MR (DRMR) | Non-linear, non-parametric | 5+ strata | Non-parametric stratification (Tian 2023 PLoS Genet 19:e1010823); replaces residual stratification when linearity fails | Continuous exposures only; needs individual-level data; Hamilton 2023 medRxiv 23293658 shows stratum-specific bias from age/sex | Methodology evolves; benchmark consensus shifts every 2-3 years. Verify against the current Slob & Burgess 2020 *Genet Epidemiol*, Burgess 2023 *Wellcome Open Res* "Guidelines for performing Mendelian randomization" (v3+), and STROBE-MR 2021 reporting standards before locking a method as primary. ## Decision Tree by Experimental Scenario | Scenario | Primary method | Sensitivity battery | Why | |----------|----------------|----------------------|-----| | Standard two-sample, independent cohorts, polygenic exposure | IVW (random) | Egger + weighted median + MR-PRESSO + MR-RAPS + Steiger | Default; covers balanced, directional, outlier, weak-IV regimes | | One-sample (e.g. UK Biobank both ends) | IVW with weak-IV-aware (MR-RAPS) | Egger + LCV + jackknife SE | One-sample F-stat floor shifts to F >= 20; jackknife SE preferred over analytic at one-sample scale; do NOT run exposure GWAS and outcome GWAS on the same individuals then claim two-sample (Barry 2021 PLoS Genet 17:e1009703 collider bias); within-stratum MR (e.g. "MR among smokers") risks collider bias from the stratification variable | | Partial sample overlap (UKB exposure + UKB outcome) | MR-RAPS with overlap correction | Sample-overlap-adjusted IVW (Burgess 2016 Genet Epidemiol 40:597) | Bias is intermediate, proportional to z-score correlation | | Drug-target / cis-MR (cis-pQTL, cis-eQTL) | IVW restricted to cis window | Colocalization PP.H4 + LD-prune within window | Exclusion restriction relaxed because the protein/transcript directly mediates effect (Schmidt 2020 Nat Commun 11:3255) | | MVMR for measured pleiotropy (e.g. LDL adjusted for HDL/TG) | `MVMR::ivw_mvmr` | Conditional F + Q_A heterogeneity | Required when exposures correlate via shared SNPs | | Mediation MR (X -> M -> Y) | MVMR difference of total vs direct | Two-step MR + product-of-coefficients (Carter 2021 Eur J Epidemiol 36:465) | Network MR; quantifies indirect effect | | Polygenic exposure with potential CHP (e.g. BMI -> CHD) | CAUSE (primary) + IVW (secondary) | Egger + MR-PRESSO + LCV | CAUSE explicitly models CHP via shared-factor; needs >=100 sig SNPs | | Binary outcome (e.g. T2D) on linear scale | IVW on log-OR with log-additive coding | All sensitivity on log-OR; report exp(beta) | Linearity of MR estimating equation holds on log-OR not OR | | Time-to-event (Cox) outcome | IVW on log-HR | Non-collapsibility-aware log-HR reporting | Non-collapsibility caveats apply | | Non-linear MR (e.g. alcohol J-curve) | DRMR (Tian 2023) + residual stratification side-by-side | Negative-control outcomes (genotype-vs-sex within strata); Hamilton 2023 limitation cited | Both methods produce stratum-specific bias from age/sex effects (Hamilton 2023 medRxiv 23293658); pre-specify the non-linear hypothesis, do not data-snoop the J-curve, report negative-control sanity checks | | Single-patient rare disease | Not MR -- use FRASER/DROP outlier framework | See alternative-splicing/outlier-splicing-detection | MR requires summary stats; n=1 is wrong regime | ## One-Sample vs Two-Sample Bias Direction | Design | Weak-IV bias direction | Reason | |--------|------------------------|--------| | One-sample, F<10 | Toward confounded observational estimate (overestimates causal effect if confounding is in same direction) | Sample correlation between IV-X and IV-Y residuals | | Two-sample non-overlapping, F<10 | Toward null | Independent samples decouple residuals (Burgess 2011 IJE 40:755) | | Two-sample with partial overlap | Intermediate; proportional to overlap fraction and z-score correlation | Burgess 2016 Genet Epidemiol 40:597; correction available | **Operational rule:** Whenever both GWAS came from UK Biobank (or any single biobank), treat the analysis as one-sample-equivalent and prefer MR-RAPS as primary. Treating it as "two-sample because separate GWAS files" is a common error and produces overestimates. ### MRlap: unified correction for sample overlap + winner's curse + weak instruments MRlap (Mounier & Kutalik 2023 Genet Epidemiol 47:314) jointly corrects three biases that previously required three separate tools: sample overlap, winner's curse, and weak-instrument bias. It builds on an LDSC scaffold (cross-trait LD-score regression intercept estimates the overlap-induced covariance) and reweights the IVW estimate against the analytical bias-correction formula. ```r remotes::install_github('n-mounier/MRlap') # never on CRAN; bioconductor unsuitable library(MRlap) fit <- MRlap( exposure = gwas_X_df, exposure_name = 'BMI', outcome = gwas_Y_df, outcome_name = 'T2D', ld = 'eur_w_ld_chr/', hm3 = 'w_hm3.snplist', # LDSC reference files MR_threshold = 5e-8, MR_pruning_dist = 500, MR_pruning_LD = 0.05 ) fit$MRcorrection$corrected_effect # overlap + winner's curse + weak-IV corrected fit$MRcorrection$corrected_effect_se fit$LDSC$h2_exp # exposure heritability sanity check fit$LDSC$int_crosstrait # cross-trait LDSC intercept; ~0 means no sample overlap ``` **Decision rule -- prefer MRlap when:** (a) any sample overlap is suspected, (b) only sumstats are available (no individual-level data for re-running GWAS on disjoint samples), (c) exposure discovery and outcome were both run inside the same biobank (UKB-on-UKB, FinnGen-on-FinnGen). MRlap returns NA / unstable estimates when h^2 < 0.05; in that regime, fall back to Burgess 2016 overlap-corrected IVW plus MR-RAPS for the weak-IV component. ## Drug-Target / cis-MR Framework cis-MR restricts instruments to the cis-regulatory window of the gene encoding the protein/transcript exposure (Schmidt 2020 Nat Commun 11:3255), relaxing the exclusion-restriction assumption because the protein product directly mediates the SNP's effect on the outcome. Operational core: extract cis-pQTL/cis-eQTL within +/-500 kb of the gene; clump at r2 < 0.1 (looser than polygenic MR to retain power within a narrow window); require colocalization PP.H4 >= 0.7; flag protein-altering variants (PAV) which can break SomaScan/Olink aptamer/antibody binding rather than reflect biology. Full drug-target cis-MR workflow including UKB-PPP / deCODE / Fenland pQTL panels, PAV flagging, Olink vs SomaScan replication (~15-30% cross-platform disagreement), and the operational claim ladder lives in causal-genomics/proteome-mr-drug-target. Use that skill for any drug-target nomination. ### Binary outcomes and non-collapsibility MR with logistic-GWAS sumstats returns per-allele log-OR on the **population-averaged** (marginal) scale, NOT the conditional log-OR (Burgess 2017 Stat Methods Med Res 26:2333). For rare disease (prevalence < 10%), OR ~= RR ~= HR and the distinction is harmless. For common disease, OR diverges from RR/HR and the MR estimate cannot be back-converted to a conditional effect without strong assumptions; report as "per 1-SD increase in genetically-predicted X, OR for Y = ..." rather than implying an individual-level intervention effect. Collider bias when conditioning on a collider variable (Coscia 2022 Eur J Epidemiol 37:671 formalizes this for stratified MR): case-only or disease-progression designs condition the sample on disease status, opening a collider path between any cause of disease and any cause of progression. MR within affected subsets without explicit adjustment for selection probability is fragile; weight by inverse probability of selection or restrict claims to the unconditioned population. ## Per-Method Failure Modes ### IVW under directional pleiotropy **Trigger:** Several SNPs affect the outcome through pathways not via the exposure, in a consistent direction. **Mechanism:** IVW is a weighted regression through origin; non-zero mean pleiotropy shifts the slope. **Symptom:** Egger intercept p < 0.05 with non-zero estimate; IVW differs from weighted median; MR-PRESSO global test p < 0.05. **Fix:** Use Egger (if `I^2_GX >= 0.9` -- otherwise SIMEX-correct via the `simex` package applied to the Egger fit, treating `se.exposure` as measurement error in `beta.exposure`); cross-check with weighted median, MR-PRESSO, and CAUSE; report IVW only as one of a panel, never alone. The `MendelianRandomization::mr_egger()` function accepts `distribution='normal'` and reports the `I.sq` (I^2_GX) NOME diagnostic but applies no NOME/SIMEX correction to the estimate itself and does NOT expose a SIMEX wrapper. ### Weak-instrument bias direction **Trigger:** Mean per-instrument F-statistic < 10, or several individual F < 10. **Mechanism:** Weak IVs amplify finite-sample correlation between IV-X and IV-Y errors; bias direction depends on overlap regime (see table above). **Symptom:** Estimates shift markedly when removing the weakest instruments; one-sample MR estimates much larger than two-sample. **Fix:** Compute F per instrument from the EXPOSURE GWAS, not the outcome; exclude F < 10; use MR-RAPS (handles weak IVs by design); for two-sample, also report unweighted IVW (less weak-IV-bias-inflated than weighted in some regimes). ### Winner's curse at P~5e-8 **Trigger:** Discovery GWAS is the source of both instrument selection and effect-size estimates. **Mechanism:** SNPs that just cross 5e-8 in discovery have over-estimated effect sizes (regression toward the mean in independent replication); MR uses inflated `beta_X`, biasing causal estimate. **Symptom:** MR effect shrinks substantially when using effect sizes from an independent replication GWAS. **Fix:** (1) Three-sample design (discovery / replication-for-instrument-effect / outcome) where feasible. (2) When sumstats-only: MRlap (Mounier 2023), MR-SimSS (sample-splitting from sumstats), or RIVW (Ma 2023 Ann Statist 51:211 -- rerandomized IVW) jointly correct winner's curse + weak IVs + overlap. (3) Jiang 2023 IJE 52:1209 empirical magnitude: variant-level inflation ~50-400% near the genome-wide-significance threshold, dropping to <25% when the minimum P <= 1e-13. ### NOME violation invalidating Egger **Trigger:** Running MR-Egger with `I^2_GX < 0.9`. **Mechanism:** Egger assumes NO Measurement Error in exposure effect sizes (NOME); when violated, Egger slope is attenuated toward null with reciprocal bias on the intercept. **Symptom:** `mr_pleiotropy_test()` Egger estimate disagrees with weighted median in magnitude but agrees in direction; `Isq()` function returns <0.9. **Fix:** Compute `Isq(beta_X, se_X)` (Bowden 2016 IJE 45:1961); if <0.9, apply SIMEX correction via `simex` package or report Egger as exploratory only. The MendelianRandomization package's `mr_egger()` reports the `I.sq` (I^2_GX) NOME diagnostic but does not itself apply a NOME/SIMEX correction; SIMEX must be run separately. ### Steiger filter false flag under unmeasured confounding **Trigger:** Applying `steiger_filtering()` on traits with unmeasured shared confounders (e.g. SES). **Mechanism:** Steiger compares variance explained in exposure vs outcome per SNP; an unmeasured confounder upstream of both produces SNPs that explain more variance in the outcome than the exposure, falsely flagging "reverse causation" (Lutz 2022 Genet Epidemiol 46:139). **Symptom:** Many SNPs flagged as wrong-direction yet biology and prior MR support forward causation. **Fix:** Treat Steiger as a heuristic, not gospel; cross-validate direction with bidirectional MR (forward + reverse with independent instrument sets); for known-confounder-rich domains (psychiatric traits, SES proxies) use LCV or LHC-MR instead, which jointly model confounders. ### Palindromic SNP harmonization **Trigger:** SNPs with alleles A/T or C/G near MAF 0.5. **Mechanism:** Strand orientation is ambiguous for palindromic SNPs when allele frequencies are intermediate; flipping introduces sign errors that look like pleiotropy. **Symptom:** `harmonise_data()` reports many palindromic SNPs dropped; remaining SNPs show heterogeneity from a handful. **Fix:** Default `action = 2` (infer from allele frequencies) drops MAF~0.5 palindromes; `action = 3` drops ALL palindromes (most conservative); never use `action = 1` (assumes forward strand) unless both GWAS are guaranteed to use the same strand convention. Document choice in methods. ## Quantitative Thresholds | Threshold | Source | Rationale | |-----------|--------|-----------| | F-statistic > 10 per instrument | Staiger & Stock 1997 (linear IV) | Heuristic; debated (Burgess 2011 IJE; Zhao 2020 argues 10 is too low for one-sample) | | Conditional F > 10 per exposure (MVMR) | Sanderson 2019 IJE 48:713 | Total F can be high while conditional F low; per-exposure F is what matters | | `I^2_GX >= 0.9` for Egger | Bowden 2016 IJE 45:1961 | NOME assumption; below this, SIMEX correction required | | CAUSE >= 100 significant SNPs | Morrison 2020 Nat Genet 52:740 | Mixture model needs signal density for shared-factor estimation | | Egger >= 10 instruments | Bowden 2015 IJE 44:512 | Power for slope test in weighted regression | | Sample-overlap z-score correlation | Burgess 2016 Genet Epidemiol 40:597 | Use LDSC bivariate intercept as proxy; correct IVW SE accordingly | | Clumping r2 < 0.001, 10 Mb window | TwoSampleMR default; matches GWAS LD norms | Polygenic MR; cis-MR uses r2 < 0.1 within window | | Steiger p < 0.05 | Hemani 2017 PLoS Genet 13:e1007081 | Heuristic; subject to confounder caveat (Lutz 2022) | | MR-PRESSO NbDistribution | 1000 exploratory; >= 5000 publication; >= 10000 stringent | Verbanck 2018 Nat Genet 50:693; precision of global p-value scales with NbDistribution |