Mendelian randomization

Mendelian randomization introduction

Comparison between RCT and MR

Instrumental Variables (IV)

Instrumental variable (IV) can be defined as a variable that is correlated with the exposure X and uncorrelated with the error \(\epsilon\) in the following regression:

\[ Y = X\beta + \epsilon \]

\(Y\) is the outcome

Assumptions

Key Assumptions

Assumptions	Description
Relevance	Instrumental variables are strongly associated with the exposure.
Exclusion restriction	Instrumental variables do not affect the outcome except through the exposure.
Independence	There are no confounders of the instrumental variables and the outcome.
Monotonicity	Variants affect the exposure in the same direction for all individuals
No assortative mating	Assortative mating might cause bias in MR

Two-stage least-squares (2SLS)

\[ X = \mu_1 + \beta_{IV} IV + \epsilon_1 \]

\[ Y = \mu_2 + \beta_{2SLS} \hat{X} + \epsilon_2 \]

Two-sample MR

Two-sample MR refers to the approach that the genetic effects of the instruments on the exposure can be estimated in an independent sample other than that used to estimate effects between instruments on the outcome. As more and more GWAS summary statistics become publicly available, the scope of MR also expands with Two-sample MR methods.

\[ \hat{\beta}_{X,Y} = {{\hat{\beta}_{IV,Y}}\over{\hat{\beta}_{IV,X}}} \]

Caveats

For two-sample MR, there is an additional key assumption:

The two samples used for MR are from the same underlying populations. (The effect size of instruments on exposure should be the same in both samples.)

Therefore, for two-sample MR, we usually use datasets from similar non-overlapping populations in terms of not only ancestry but also contextual factors.

IV selection

One of the first things to do when you plan to perform any type of MR is to check the associations of instrumental variables with the exposure to avoid bias caused by weak IVs.

The most commonly used method here is the F-statistic, which tests the association of instrumental variables with the exposure.

Practice

In this tutorial, we will walk you through how to perform a minimal TwoSampleMR analysis. We will use the R package TwoSampleMR, which provides easy-to-use functions for formatting, clumping and harmonizing GWAS summary statistics.

This package integrates a variety of commonly used MR methods for analysis, including:

> mr_method_list()
                             obj
1                  mr_wald_ratio
2               mr_two_sample_ml
3            mr_egger_regression
4  mr_egger_regression_bootstrap
5               mr_simple_median
6             mr_weighted_median
7   mr_penalised_weighted_median
8                         mr_ivw
9                  mr_ivw_radial
10                    mr_ivw_mre
11                     mr_ivw_fe
12                mr_simple_mode
13              mr_weighted_mode
14         mr_weighted_mode_nome
15           mr_simple_mode_nome
16                       mr_raps
17                       mr_sign
18                        mr_uwr

                                                        name PubmedID
1                                                 Wald ratio
2                                         Maximum likelihood
3                                                   MR Egger 26050253
4                                       MR Egger (bootstrap) 26050253
5                                              Simple median
6                                            Weighted median
7                                  Penalised weighted median
8                                  Inverse variance weighted
9                                                 IVW radial
10 Inverse variance weighted (multiplicative random effects)
11                 Inverse variance weighted (fixed effects)
12                                               Simple mode
13                                             Weighted mode
14                                      Weighted mode (NOME)
15                                        Simple mode (NOME)
16                      Robust adjusted profile score (RAPS)
17                                     Sign concordance test
18                                     Unweighted regression

File Preparation

To perform two-sample MR analysis, we need summary statistics for exposure and outcome generated from independent populations with the same ancestry.

In this tutorial, we will use sumstats from Biobank Japan pheweb and KoGES pheweb.

Type 2 diabetes sumstats from BBJ : wget -O bbj_t2d.zip https://pheweb.jp/download/T2D
BMI sumstats from KoGES : wget -O koges_bmi.txt.gz https://koges.leelabsg.org/download/KoGES_BMI

R package TwoSampleMR

First, to use TwosampleMR, we need R>= 4.1. To install the package, run:

library(remotes)
install_github("MRCIEU/TwoSampleMR")

Loading package

library(TwoSampleMR)

Reading exposure sumstats

#format exposures dataset

exp_raw <- fread("koges_bmi.txt.gz")

Extracting instrumental variables

# select only significant variants
exp_raw <- subset(exp_raw,exp_raw$pval<5e-8)

exp_dat <- format_data( exp_raw,
    type = "exposure",
    snp_col = "rsids",
    beta_col = "beta",
    se_col = "sebeta",
    effect_allele_col = "alt",
    other_allele_col = "ref",
    eaf_col = "af",
    pval_col = "pval",
)

Clumping exposure variables

clumped_exp <- clump_data(exp_dat,clump_r2=0.01,pop="EAS")

outcome

out_raw <- fread("hum0197.v3.BBJ.T2D.v1/GWASsummary_T2D_Japanese_SakaueKanai2020.auto.txt.gz",
                    select=c("SNPID","Allele1","Allele2","BETA","SE","p.value"))
out_dat <- format_data( out_raw,
    type = "outcome",
    snp_col = "SNPID",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "Allele2",
    other_allele_col = "Allele1",
    pval_col = "p.value",
)

Harmonizing data

harmonized_data <- harmonise_data(clumped_exp,out_dat,action=1)

Perform MR analysis

res <- mr(harmonized_data)

id.exposure id.outcome  outcome exposure    method  nsnp    b   se  pval
<chr>   <chr>   <chr>   <chr>   <chr>   <int>   <dbl>   <dbl>   <dbl>
9J8pv4  IyUv6b  outcome exposure    MR Egger    28  1.3337580   0.69485260  6.596064e-02
9J8pv4  IyUv6b  outcome exposure    Weighted median 28  0.6298980   0.09401352  2.083081e-11
9J8pv4  IyUv6b  outcome exposure    Inverse variance weighted   28  0.5598956   0.23225806  1.592361e-02
9J8pv4  IyUv6b  outcome exposure    Simple mode 28  0.6097842   0.15180476  4.232158e-04
9J8pv4  IyUv6b  outcome exposure    Weighted mode   28  0.5946778   0.12820220  8.044488e-05

Sensitivity analysis

Heterogeneity

mr_heterogeneity(harmonized_data)

id.exposure id.outcome  outcome exposure    method  Q   Q_df    Q_pval
<chr>   <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>
9J8pv4  IyUv6b  outcome exposure    MR Egger    670.7022    26  1.000684e-124
9J8pv4  IyUv6b  outcome exposure    Inverse variance weighted   706.6579    27  1.534239e-131

Horizontal Pleiotropy

Intercept in MR-Egger

mr_pleiotropy_test(harmonized_data)

id.exposure id.outcome  outcome exposure    egger_intercept se  pval
<chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>
9J8pv4  IyUv6b  outcome exposure    -0.03603697 0.0305241   0.2484472

Single SNP MR and leave-one-out MR

Single SNP MR

res_single <- mr_singlesnp(harmonized_data)
res_single

exposure    outcome id.exposure id.outcome  samplesize  SNP b   se  p
<chr>   <chr>   <chr>   <chr>   <lgl>   <chr>   <dbl>   <dbl>   <dbl>
1   exposure    outcome 9J8pv4  IyUv6b  NA  rs10198356  0.6323140   0.2082837   2.398742e-03
2   exposure    outcome 9J8pv4  IyUv6b  NA  rs10209994  0.9477808   0.3225814   3.302164e-03
3   exposure    outcome 9J8pv4  IyUv6b  NA  rs10824329  0.6281765   0.3246214   5.297739e-02
4   exposure    outcome 9J8pv4  IyUv6b  NA  rs10938397  1.2376316   0.2775854   8.251150e-06
5   exposure    outcome 9J8pv4  IyUv6b  NA  rs11066132  0.6024303   0.2232401   6.963693e-03
6   exposure    outcome 9J8pv4  IyUv6b  NA  rs12522139  0.2905201   0.2890240   3.148119e-01
7   exposure    outcome 9J8pv4  IyUv6b  NA  rs12591730  0.8930490   0.3076687   3.700413e-03
8   exposure    outcome 9J8pv4  IyUv6b  NA  rs13013021  1.4867889   0.2207777   1.646925e-11
9   exposure    outcome 9J8pv4  IyUv6b  NA  rs1955337   0.5442640   0.2994146   6.910079e-02
10  exposure    outcome 9J8pv4  IyUv6b  NA  rs2076308   1.1176226   0.2657969   2.613132e-05
11  exposure    outcome 9J8pv4  IyUv6b  NA  rs2278557   0.6238587   0.2968184   3.556906e-02
12  exposure    outcome 9J8pv4  IyUv6b  NA  rs2304608   1.5054682   0.2968905   3.961740e-07
13  exposure    outcome 9J8pv4  IyUv6b  NA  rs2531995   1.3972908   0.3130157   8.045689e-06
14  exposure    outcome 9J8pv4  IyUv6b  NA  rs261967    1.5303384   0.2921192   1.616714e-07
15  exposure    outcome 9J8pv4  IyUv6b  NA  rs35332469  -0.2307314  0.3479219   5.072217e-01
16  exposure    outcome 9J8pv4  IyUv6b  NA  rs35560038  -1.5730870  0.2018968   6.619637e-15
17  exposure    outcome 9J8pv4  IyUv6b  NA  rs3755804   0.5314915   0.2325073   2.225933e-02
18  exposure    outcome 9J8pv4  IyUv6b  NA  rs4470425   0.6948046   0.3079944   2.407689e-02
19  exposure    outcome 9J8pv4  IyUv6b  NA  rs476828    1.1739083   0.1568550   7.207355e-14
20  exposure    outcome 9J8pv4  IyUv6b  NA  rs4883723   0.5479721   0.2855004   5.494141e-02
21  exposure    outcome 9J8pv4  IyUv6b  NA  rs509325    0.5491040   0.1598196   5.908641e-04
22  exposure    outcome 9J8pv4  IyUv6b  NA  rs55872725  1.3501891   0.1259791   8.419325e-27
23  exposure    outcome 9J8pv4  IyUv6b  NA  rs6089309   0.5657525   0.3347009   9.096620e-02
24  exposure    outcome 9J8pv4  IyUv6b  NA  rs6265  0.6457693   0.1901871   6.851804e-04
25  exposure    outcome 9J8pv4  IyUv6b  NA  rs6736712   0.5606962   0.3448784   1.039966e-01
26  exposure    outcome 9J8pv4  IyUv6b  NA  rs7560832   0.6032080   0.2904972   3.785077e-02
27  exposure    outcome 9J8pv4  IyUv6b  NA  rs825486    -0.6152759  0.3500334   7.878772e-02
28  exposure    outcome 9J8pv4  IyUv6b  NA  rs9348441   -4.9786332  0.2572782   1.992909e-83
29  exposure    outcome 9J8pv4  IyUv6b  NA  All - Inverse variance weighted 0.5598956   0.2322581   1.592361e-02
30  exposure    outcome 9J8pv4  IyUv6b  NA  All - MR Egger  1.3337580   0.6948526   6.596064e-02

leave-one-out MR

res_loo <- mr_leaveoneout(harmonized_data)
res_loo

exposure    outcome id.exposure id.outcome  samplesize  SNP b   se  p
<chr>   <chr>   <chr>   <chr>   <lgl>   <chr>   <dbl>   <dbl>   <dbl>
1   exposure    outcome 9J8pv4  IyUv6b  NA  rs10198356  0.5562834   0.2424917   2.178871e-02
2   exposure    outcome 9J8pv4  IyUv6b  NA  rs10209994  0.5520576   0.2388122   2.079526e-02
3   exposure    outcome 9J8pv4  IyUv6b  NA  rs10824329  0.5585335   0.2390239   1.945341e-02
4   exposure    outcome 9J8pv4  IyUv6b  NA  rs10938397  0.5412688   0.2388709   2.345460e-02
5   exposure    outcome 9J8pv4  IyUv6b  NA  rs11066132  0.5580606   0.2417275   2.096381e-02
6   exposure    outcome 9J8pv4  IyUv6b  NA  rs12522139  0.5667102   0.2395064   1.797373e-02
7   exposure    outcome 9J8pv4  IyUv6b  NA  rs12591730  0.5524802   0.2390990   2.085075e-02
8   exposure    outcome 9J8pv4  IyUv6b  NA  rs13013021  0.5189715   0.2386808   2.968017e-02
9   exposure    outcome 9J8pv4  IyUv6b  NA  rs1955337   0.5602635   0.2394505   1.929468e-02
10  exposure    outcome 9J8pv4  IyUv6b  NA  rs2076308   0.5431355   0.2394403   2.330758e-02
11  exposure    outcome 9J8pv4  IyUv6b  NA  rs2278557   0.5583634   0.2394924   1.972992e-02
12  exposure    outcome 9J8pv4  IyUv6b  NA  rs2304608   0.5372557   0.2377325   2.382639e-02
13  exposure    outcome 9J8pv4  IyUv6b  NA  rs2531995   0.5419016   0.2379712   2.277590e-02
14  exposure    outcome 9J8pv4  IyUv6b  NA  rs261967    0.5358761   0.2376686   2.415093e-02
15  exposure    outcome 9J8pv4  IyUv6b  NA  rs35332469  0.5735907   0.2378345   1.587739e-02
16  exposure    outcome 9J8pv4  IyUv6b  NA  rs35560038  0.6734906   0.2217804   2.391474e-03
17  exposure    outcome 9J8pv4  IyUv6b  NA  rs3755804   0.5610215   0.2413249   2.008503e-02
18  exposure    outcome 9J8pv4  IyUv6b  NA  rs4470425   0.5568993   0.2392632   1.993549e-02
19  exposure    outcome 9J8pv4  IyUv6b  NA  rs476828    0.5037555   0.2443224   3.922224e-02
20  exposure    outcome 9J8pv4  IyUv6b  NA  rs4883723   0.5602050   0.2397325   1.945000e-02
21  exposure    outcome 9J8pv4  IyUv6b  NA  rs509325    0.5608429   0.2468506   2.308693e-02
22  exposure    outcome 9J8pv4  IyUv6b  NA  rs55872725  0.4419446   0.2454771   7.180543e-02
23  exposure    outcome 9J8pv4  IyUv6b  NA  rs6089309   0.5597859   0.2388902   1.911519e-02
24  exposure    outcome 9J8pv4  IyUv6b  NA  rs6265  0.5547068   0.2436910   2.282978e-02
25  exposure    outcome 9J8pv4  IyUv6b  NA  rs6736712   0.5598815   0.2387602   1.902944e-02
26  exposure    outcome 9J8pv4  IyUv6b  NA  rs7560832   0.5588113   0.2396229   1.969836e-02
27  exposure    outcome 9J8pv4  IyUv6b  NA  rs825486    0.5800026   0.2367545   1.429330e-02
28  exposure    outcome 9J8pv4  IyUv6b  NA  rs9348441   0.7378967   0.1366838   6.717515e-08
29  exposure    outcome 9J8pv4  IyUv6b  NA  All 0.5598956   0.2322581   1.592361e-02

Visualization

Scatter plot

res <- mr(harmonized_data)
p1 <- mr_scatter_plot(res, harmonized_data)
p1[[1]]

Single SNP

res_single <- mr_singlesnp(harmonized_data)
p2 <- mr_forest_plot(res_single)
p2[[1]]

Leave-one-out

res_loo <- mr_leaveoneout(harmonized_data)
p3 <- mr_leaveoneout_plot(res_loo)
p3[[1]]

Funnel plot

res_single <- mr_singlesnp(harmonized_data)
p4 <- mr_funnel_plot(res_single)
p4[[1]]

MR Steiger directionality test

MR Steiger directionality test is a method to test the causal direction.

harmonized_data$"r.outcome" <- get_r_from_lor(
  harmonized_data$"beta.outcome",
  harmonized_data$"eaf.outcome",
  45383,
  132032,
  0.26,
  model = "logit",
  correction = FALSE
)

out <- directionality_test(harmonized_data)
out

id.exposure id.outcome  exposure    outcome snp_r2.exposure snp_r2.outcome  correct_causal_direction    steiger_pval
<chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <lgl>   <dbl>
rvi6Om  ETcv15  BMI T2D 0.02125453  0.005496427 TRUE    NA

Reference: Hemani, G., Tilling, K., & Davey Smith, G. (2017). Orienting the causal relationship between imprecisely measured traits using GWAS summary data. PLoS genetics, 13(11), e1007081.

MR-Base (web app)

MR-Base web app

Reference

Sanderson, E., Glymour, M. M., Holmes, M. V., Kang, H., Morrison, J., Munafò, M. R., ... & Davey Smith, G. (2022). Mendelian randomization. Nature Reviews Methods Primers, 2(1), 1-21.
Hemani, G., Zheng, J., Elsworth, B., Wade, K. H., Haberland, V., Baird, D., ... & Haycock, P. C. (2018). The MR-Base platform supports systematic causal inference across the human phenome. elife, 7, e34408.

Mendelian randomization