How to do Mendelian randomization: from GWAS data to sensitivity analyses

Workflow

Seven steps of two-sample MR: functions, key parameters and what to report

Functions and defaults were checked against TwoSampleMR 0.7.12 (released 2026-10-08) and ieugwasr 1.2.0.

Core assumptions

The three core assumptions: how to test each, and how to argue the ones you cannot test

Only relevance can be tested directly. Independence and exclusion restriction can be checked only in part, and the paper has to argue them with sensitivity analyses and biological knowledge. STROBE-MR item 5 asks for all three to be stated explicitly in the methods.

Of the 20 STROBE-MR items (JAMA 2021 statement and BMJ explanation), reviewers most often check: item 5, stating the three core assumptions; item 7, how the assumptions were assessed; item 9, software, package versions and settings, and whether the protocol was preregistered; item 10d, for two-sample MR, the similarity of the exposure and outcome samples and the number of overlapping participants; item 11, results on an interpretable scale (for example OR per 1 SD); and item 13, sensitivity analyses and tests of direction. Writing item by item avoids most comments about incomplete methods.

The 2023 update of the Burgess et al. MR guidelines also recommends a positive control outcome: an outcome for which the causal effect is established, used to check that the instruments reproduce it. The LDL-C to coronary disease analysis on this page is itself a recognised positive control and is a good way to check that a pipeline works.

Data sources

Where to get GWAS summary statistics: status of each source in October 2026

The table reflects access checks on 2026-10-10. In the methods, state for each source: dataset ID or file name, version, ancestry, sample size, genome build and effect allele column.

Direction of bias from sample overlap: without overlap, weak instrument bias pulls the estimate towards zero and does not create false positives; with overlap, the bias points towards the observational association and grows linearly with the overlap fraction (Burgess et al. 2023 guidelines). With F around 30, even complete overlap gives a relative bias of only about 1/F, roughly 3%. The usual way to avoid it is to take the exposure from European consortia or UKB and the outcome from FinnGen; when overlap cannot be avoided, report its extent and correct with methods such as MRlap.

OpenGWAS

Setting up the OpenGWAS token and fixing 401 errors

Since 0.6.0 (2024-04), TwoSampleMR depends on the ieugwasr release with the new authentication system, and since 0.6.20 it requires ieugwasr ≥ 1.1.0. Older versions cannot connect even with a token.

Measured here (ieugwasr 1.2.0): without a token, extract_instruments returns the 401 message in the first row; ld_clump without bfile goes through the API and also needs a token. The variable name IUEUGWAS_TOKEN that circulates in Chinese posts is wrong; ieugwasr only reads OPENGWAS_JWT.

# 1. Sign in at https://api.opengwas.io/profile and generate a token (valid for 14 days)
# 2. Add it to ~/.Renviron (or .Renviron in the project directory); the variable must be OPENGWAS_JWT
usethis::edit_r_environ()
# Add one line to the opened file: OPENGWAS_JWT=eyJhbGciOi... (no quotes, end with a newline)

# 3. Restart R and check
ieugwasr::get_opengwas_jwt()   # a long string means the token was read
ieugwasr::user()               # account info means the token works; 401 means expired or mis-copied

# 4. The current argument name is opengwas_jwt; access_token from older tutorials no longer exists
exp <- TwoSampleMR::extract_instruments("ieu-a-300", p1 = 5e-8, r2 = 0.001, kb = 10000)

Offline data

No token needed: download summary statistics and read them locally

The offline route suits batch analyses and papers that must be reproducible: data files, reference panel and scripts stay local, and results do not change when the OpenGWAS database is updated.

Found in this page's test: the GLGC 2013 LDL-C file contains three P values below the smallest double (1.24e-652, 6.56e-397, 3.85e-326), so fread reads the whole P column as text. pval < 5e-8 then compares strings: 2,435,913 of 2,437,751 rows "pass", clumping yields 1,837 instruments, and nothing raises an error. Check the P column with class() after reading; as.numeric turns those three values into 0, and format_data then truncates them at min_pval = 1e-200.

When format_data meets duplicated rsIDs it keeps the first row and warns. Multi-allelic sites in files such as FinnGen produce duplicated rsIDs (rs7534572 in this test), so match ref/alt to the exposure alleles before calling format_data.

library(TwoSampleMR); library(data.table)   # fread on .gz files also needs R.utils

## A. Exposure: text summary statistics from the authors' site or GWAS Catalog (GLGC 2013 LDL-C)
ldl <- fread("jointGwasMc_LDL.txt.gz")
setnames(ldl, c("P-value", "Freq.A1.1000G.EUR"), c("pval", "eaf"))
class(ldl$pval)              # measured here: "character", because the file contains values such as 1.24e-652
ldl[, pval := as.numeric(pval)]                 # without conversion, pval < 5e-8 is a string comparison
ldl[, `:=`(A1 = toupper(A1), A2 = toupper(A2))] # GLGC alleles are lower case; A1 is the effect allele
exp_dat <- format_data(as.data.frame(ldl[pval < 5e-8]), type = "exposure",
  snp_col = "rsid", beta_col = "beta", se_col = "se", eaf_col = "eaf",
  effect_allele_col = "A1", other_allele_col = "A2", pval_col = "pval", samplesize_col = "N")

## B. Outcome: GWAS Catalog harmonised file (*.h.tsv.gz, GRCh38, hm_ columns share one orientation)
# Keep only instrument rows with awk first: about 80 s for the 370 MB file, so R never reads 8.6 million rows
#   gzip -dc GCST003116.h.tsv.gz | awk -F'\t' 'NR==FNR{a[$1];next} FNR==1 || ($2 in a)' snps.txt - > cad_subset.tsv
cad <- fread("cad_subset.tsv")
out_dat <- format_data(as.data.frame(cad), type = "outcome", snps = exp_dat$SNP,
  snp_col = "hm_rsid", beta_col = "hm_beta", se_col = "standard_error",
  eaf_col = "hm_effect_allele_frequency", effect_allele_col = "hm_effect_allele",
  other_allele_col = "hm_other_allele", pval_col = "p_value", chr_col = "hm_chrom", pos_col = "hm_pos")

## C. IEU OpenGWAS VCF: FORMAT is ES:SE:LP:AF:SS:ID, LP is -log10(P), the effect allele is ALT
#   bcftools query -f '%ID\t%CHROM\t%POS\t%ALT\t%REF\t[%ES]\t[%SE]\t[%LP]\t[%AF]\t[%SS]\n' ieu-a-300.vcf.gz > ieu-a-300.tsv
v <- fread("ieu-a-300.tsv", col.names = c("SNP","chr","pos","effect_allele","other_allele","beta","se","lp","eaf","samplesize"))
exp_vcf <- format_data(as.data.frame(v), type = "exposure", pval_col = "lp", log_pval = TRUE)

## D. Pan-UKB: the P column is neglog10_pval_EUR, so also use log_pval = TRUE; alt is the effect allele, GRCh37

# FinnGen files are tabix-indexed: fetch only the instrument positions instead of the full 810 MB file
# regions_hg38.tsv has two columns: chromosome (1-23, no "chr") and GRCh38 position,
# e.g. from hm_chrom and hm_pos in the GWAS Catalog harmonised file
tabix -h -R regions_hg38.tsv \
  https://storage.googleapis.com/finngen-public-data-r13/summary_stats/finngen_R13_I9_CHD.gz \
  > finngen_chd_subset.tsv
# Columns: #chrom pos ref alt rsids nearest_genes pval mlogp beta sebeta af_alt ... (alt is the effect allele)
# Measured here: 77 positions, 4 min 26 s, 79 rows returned (including multi-allelic sites); a .tbi file is left in the working directory

Instruments

Instrument selection: P threshold, clumping parameters and the F statistic

Fix the threshold, clumping parameters and reference panel before looking at results, and state them in the methods.

What F > 10 means: in one-sample settings or with sample overlap, F around 10 corresponds to a 2SLS bias of about 10% of the observational bias; it is a rule of thumb and does not guarantee no bias. The per-SNP F is approximately z², so at P<5e-8 F automatically exceeds 29.7, at P<5e-6 it exceeds 20.8, and at P<1e-5 it exceeds 19.5. After selecting at these thresholds, "all SNPs have F > 10" adds no information; total R², overall F and a comparison of weak SNPs when the threshold is relaxed are more useful.

The formula R² = 2·EAF·(1−EAF)·β² holds only for β in standard deviation units, and β must be squared. Some Chinese tutorials drop the square or plug in β in mg/dL, which puts R² off by orders of magnitude. The code above includes SE and N and does not depend on the units.

# Reference panel: http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz (1.57 GB, EUR/EAS/AFR/AMR/SAS)
# Extract only the population you need: tar -xzf 1kg.v3.tgz EUR.bed EUR.bim EUR.fam
# plink 1.9 from conda (bioconda::plink) or genetics.binaRies::get_plink_binary()
library(ieugwasr)
# First count significant SNPs missing from the panel: ld_clump drops them silently
bim <- data.table::fread("ref/EUR.bim", select = 2)
miss <- exp_dat[!exp_dat$SNP %in% bim$V2, ]
nrow(miss); min(miss$pval.exposure)   # measured here: 132 SNPs, smallest P 4.3e-127

clumped <- ld_clump(
  dplyr::tibble(rsid = exp_dat$SNP, pval = exp_dat$pval.exposure, id = exp_dat$id.exposure),
  clump_kb = 10000, clump_r2 = 0.001, clump_p = 1,
  bfile = "ref/EUR",                  # file prefix, without .bed
  plink_bin = Sys.which("plink"))
exp_dat <- exp_dat[exp_dat$SNP %in% clumped$rsid, ]

# Per-SNP F statistic (approximation)
exp_dat$F <- (exp_dat$beta.exposure / exp_dat$se.exposure)^2

# Variance explained R2 (form that includes SE and N; for a continuous trait with beta in SD units)
b <- exp_dat$beta.exposure; s <- exp_dat$se.exposure
f <- exp_dat$eaf.exposure;  N <- exp_dat$samplesize.exposure
exp_dat$R2 <- 2*b^2*f*(1-f) / (2*b^2*f*(1-f) + 2*s^2*N*f*(1-f))

# Overall F: k SNPs, total R2
k <- nrow(exp_dat); R2 <- sum(exp_dat$R2, na.rm = TRUE)
F_all <- R2 * (median(N) - k - 1) / (k * (1 - R2))
# Measured here (79 SNPs): per-SNP F min 27.8, median 58.4; total R2 0.085; overall F 201

Allele harmonisation

harmonise_data: the action argument, palindromic SNPs and proxies

Harmonisation makes the exposure and outcome β refer to the same effect allele. In this test, 36 of 79 SNPs needed their alleles switched; without harmonisation nearly half the SNPs would point the wrong way.

Palindromic SNPs have A/T or C/G alleles; flipping strand leaves the alleles unchanged, so only frequencies reveal the orientation. The tolerance in the harmonise_data source is 0.08, which gives the 0.42–0.58 window. Removed SNPs are not deleted from the data frame, only flagged with mr_keep = FALSE; mr() ignores them, but if you export supplementary tables without filtering, the SNP count will not match nsnp in the results.

Proxies: extract_outcome_data defaults to proxies = TRUE and rsq = 0.8 and aligns proxy alleles through LD automatically, which relies on the OpenGWAS API. Offline, find proxies in the reference panel with plink and use the in-phase output to map the alleles.

# When an instrument is missing from the outcome, find proxies with r2 >= 0.8 in the local panel
plink --bfile ref/EUR --r2 in-phase with-freqs \
  --ld-snp rs646776 --ld-window-kb 500 --ld-window 99999 --ld-window-r2 0.8 \
  --out proxy_rs646776
# The PHASE column, e.g. CG/TA, means allele C of the target SNP sits on the same haplotype as allele G of the proxy
# When using the proxy's outcome effect, map its effect allele back to the target SNP's allele with this pairing
# Measured here: about 5 s per SNP; 10 proxies for rs646776 (including rs12740374 at r2 = 1)

Main and sensitivity analyses

IVW, MR-Egger, weighted median, MR-PRESSO, Q test, leave-one-out and Steiger

The methods make different assumptions about pleiotropy. IVW is the main analysis; the others check whether its conclusion holds under different assumptions.

MR-PRESSO run time grows with the number of SNPs and NbDistribution. With 77 SNPs here, NbDistribution = 1000 took 66 seconds and warned "Outlier test unstable … The current precision is <0.077"; 2000 took 138 seconds and the warning disappeared. The outlier P values are Bonferroni-corrected for the number of SNPs, so set NbDistribution to at least 20 times the SNP count.

The Burgess et al. guidelines note that outlier-removal methods such as MR-PRESSO have very high false positive rates when several instruments are invalid. Treat the outlier-corrected result as a sensitivity analysis and report it next to the uncorrected IVW.

dat <- harmonise_data(exp_dat, out_dat, action = 2)
attr(dat, "log")                      # counts of switched and removed SNPs, for the supplement
dat <- dat[dat$mr_keep, ]             # removed rows stay in the data frame; filter before exporting tables

res <- mr(dat, method_list = c("mr_ivw", "mr_ivw_fe", "mr_egger_regression",
                               "mr_weighted_median", "mr_weighted_mode"))
generate_odds_ratios(res)
mr_heterogeneity(dat)                 # Cochran's Q
mr_pleiotropy_test(dat)               # MR-Egger intercept
Isq(dat$beta.exposure, dat$se.exposure)   # I2GX; below 0.9 MR-Egger suffers regression dilution
loo <- mr_leaveoneout(dat); ss <- mr_singlesnp(dat)

# Steiger: for a binary outcome compute r.outcome first, otherwise it is approximated as continuous
dat$r.outcome  <- get_r_from_lor(dat$beta.outcome, dat$eaf.outcome,
                                 ncase = 60801, ncontrol = 123504, prevalence = 0.05)
dat$r.exposure <- get_r_from_bsen(dat$beta.exposure, dat$se.exposure, dat$samplesize.exposure)
directionality_test(dat)

# MR-PRESSO: NbDistribution at least 20 times the number of SNPs, otherwise the outlier test lacks precision
library(MRPRESSO)
nb <- max(1000, ceiling(nrow(dat) * 20 / 1000) * 1000)
pr <- mr_presso(BetaOutcome = "beta.outcome", BetaExposure = "beta.exposure",
  SdOutcome = "se.outcome", SdExposure = "se.exposure", data = dat,
  OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = nb, seed = 2026)

# Plots
mr_scatter_plot(res, dat); mr_forest_plot(ss); mr_leaveoneout_plot(loo); mr_funnel_plot(ss)

Interpretation

How to interpret and report results when methods disagree

MR-Egger and the mode estimator have less power than IVW and wider confidence intervals, so a non-significant result from them alone is common. Look at direction and point estimates first, then P values.

Tested on this page

Tested on this page: two-sample MR of LDL-C on coronary disease

Conditions: macOS arm64, 8 cores, 16 GB RAM; R 4.5.3, TwoSampleMR 0.7.12, ieugwasr 1.2.0, MRPRESSO 1.0, MendelianRandomization 0.10.0, meta 8.5.0, PLINK 1.9, htslib 1.24. Exposure: GLGC 2013 LDL-C (2,437,751 SNPs, N up to about 173,000). Outcome: CARDIoGRAMplusC4D 2015 (GWAS Catalog GCST003116 harmonised file, 60,801 cases and 123,504 controls), replicated with FinnGen R13 I9_CHD (90,714 cases and 409,472 controls). No OpenGWAS token was used.

These results show three things. First, with a highly significant Q test, the fixed-effect IVW SE is only 46% of the random-effects SE and the P values differ by 54 orders of magnitude, so reporting fixed effects overstates precision. Second, the MR-PRESSO outliers include known pleiotropic loci such as SH2B3 (rs3184504) and ABO (rs579459); removing them moved the estimate from 0.412 to 0.437, the distortion test was not significant, and the conclusion stands. Third, the same instruments gave ORs of 1.51 and 1.32 in the two outcome sources, with a heterogeneity P of 0.022 between them; FinnGen's I9_CHD is a register-based "major coronary heart disease event" endpoint in a Finnish population, and both the outcome definition and the population change the effect size.

# Combine IVW estimates for one exposure from two independent outcome sources (meta 8.5)
library(meta)
m <- metagen(TE = c(0.4123, 0.2741), seTE = c(0.0511, 0.0318),
             studlab = c("CARDIoGRAMplusC4D 2015", "FinnGen R13 I9_CHD"),
             sm = "OR", common = TRUE, random = TRUE)
forest(m)
# Measured here: common effect OR 1.37 (1.30-1.44); random effects 1.40 (1.22-1.60); Q = 5.28, P = 0.022, I2 = 81%

1. 01

   Selection and clumping

   P<5e-8 gave 3,078 SNPs, 3,060 after de-duplication; local plink clumping (r²=0.001, 10,000 kb, 1000G EUR) took 3 seconds and kept 79. Reading and clumping took 22 seconds in total, with peak memory around 2 GB.

2. 02

   Instrument strength

   Per-SNP F minimum 27.8, median 58.4; total R² 0.085; overall F 201.

3. 03

   Harmonisation

   harmonise_data(action = 2) switched alleles for 36 SNPs and removed 2 palindromic SNPs with intermediate frequency, leaving 77.

4. 04

   Main and sensitivity analyses

   IVW gave an OR for coronary disease of 1.51 per 1 SD higher LDL-C; Egger intercept P = 0.13; Q test highly significant; MR-PRESSO found 8 outliers and the estimate barely changed after removing them. Leave-one-out estimates ranged from 0.386 to 0.455, with a largest P of 1.4×10⁻¹³. Steiger: SNP R² 0.083 for the exposure and 0.0058 for the outcome, direction correct.

5. 05

   Replication in an independent outcome

   Remote tabix on FinnGen R13 fetched the 77 positions in 4 min 26 s; IVW OR 1.32, same direction, smaller effect.

Common mistakes

Common mistakes and their consequences

Reviewer comments

Common reviewer concerns and how to respond

Chinese community

Experience from the Chinese community that holds up

The items below come from CSDN. Only posts with original error messages or the author's own tests, and consistent with official documentation or this page's tests, are included. Many 2026 "pitfall guide" posts are templated text; some give the token variable as IUEUGWAS_TOKEN or claim that action = 2 "keeps all SNPs", and both contradict the source code.

Hand it to an agent

Hand this workflow to an agent

You can describe the analysis in one sentence and let Scientify's science agent run it on a cloud computer.

Example instruction: "Run a two-sample MR with GLGC 2013 LDL-C summary statistics as the exposure and coronary artery disease from GWAS Catalog GCST003116 as the outcome: P<5e-8, local clumping with 1000G EUR (r²=0.001, 10,000 kb), compute F and R², harmonise with action=2, run IVW, Egger, weighted median, weighted mode, MR-PRESSO, Q test, leave-one-out and Steiger, replicate with FinnGen R13 I9_CHD, and write the methods paragraph and result tables following STROBE-MR."

The agent installs TwoSampleMR, MRPRESSO and plink on the cloud computer, downloads the summary statistics and reference panel, checks the P column type, effect allele columns and genome build, runs selection, clumping, harmonisation and all analyses, and produces result tables, scatter, forest, leave-one-out and funnel plots, and a draft methods paragraph. The workspace keeps scripts, parameters, logs and data files for reproduction; the agent reviews the results adversarially, for example checking the number of switched SNPs after harmonisation, whether the choice of fixed or random effects matches the Q test, and whether outlier SNPs are known pleiotropic loci. The task keeps running after you shut down your machine.

You still need to check: the ancestry and sample overlap of the exposure and outcome GWAS, whether thresholds and methods were fixed before seeing results, whether the effect scale of a binary exposure was converted correctly, and whether the wording of the conclusion matches the sensitivity analyses.