孟德尔随机化：twosampleMR 实战

孟德尔随机化：twosampleMR 实战 by 医学和生信笔记

关注公众号，发送R语言或python，可获取资料

💡专注R语言在🩺生物医学中的使用

设为“星标”，精彩不错过

两样本孟德尔随机化专用R包TwoSampleMR，可以和 IEU open GWAS 数据库无缝连接，可通过在线的方式直接查询数据库资源，并实现分析、下载等功能。

除此之外，还可以使用本地数据进行MR分析，是目前用的最多的MR分析R包，只要你做MR分析，基本上不可能绕过这个包。所以我认真学习了它的官方文档，整理下分享给大家~

本文目录：

安装
overview
基础使用
敏感性分析
Steiger方向性检验
1-to-many forest plot

例1
例2

多变量MR（MVMR）
MR-MoE
整合所有结果

安装

remotes::install_github("MRCIEU/TwoSampleMR")

overview

使用TwoSampleMR进行MR分析主要就是4步：

为暴露因素选择合适的工具变量（必要时进行连锁不平衡的过滤）
对感兴趣的结局提取工具变量
统一暴露因素和结局的数据，使得reference allele保持一致
进行MR分析、敏感性分析、可视化

基础使用

直接使用IEU的API，使用在线的方式查询、分析、提取IEU open GWAS 数据库中的工具变量。这个功能其实是借助ieugwasr实现的。

首先我们要选择合适的暴露因素和结局，暴露因素我这里选的是BMI（body mass index），结局选的是CHD（冠心病）。也就是我想研究BMI对冠心病的因果关系。

BMI在IEU open GWAS 数据库中的ID是ieu-a-2，CHD的ID是ieu-a-7。

选好之后，就是用代码查询暴露和结局相关的SNP。

在提取暴露因素相关的工具变量时，会自动进行去除连锁不平衡的SNP，由于这一步是在线进行的，需要联网，所以这一步对网络有要求，有可能会经常失败，此时我们也可以下载数据到本地进行处理（下次再讲），也可以得到一样的结果。

提取结局相关的SNP这一步也是在线进行的，也有可能会失败。

harmonise_data这一步是必须的，以确保暴露和结局的SNP的效应等位基因是一样的。

library(TwoSampleMR)

# Registered S3 method overwritten by 'data.table':
#   method           from
# print.data.table     
# TwoSampleMR version 0.5.8 
# [>] New: Option to use non-European LD reference panels for clumping etc
# [>] Some studies temporarily quarantined to verify effect allele
# [>] See news(package='TwoSampleMR') and https://gwas.mrcieu.ac.uk for further details

# List available GWASs
#ao <- available_outcomes()

# 提取暴露相关的SNP
exposure_dat <- extract_instruments("ieu-a-2", access_token = NULL)

# API: public: http://gwas-api.mrcieu.ac.uk/
#save(exposure_dat,file = "./ieu/exposure_data_ieu-a-2.rdata")

# 提取结局相关的SNP
outcome_dat <- extract_outcome_data(snps=exposure_dat$SNP, 
                                    outcomes = "ieu-a-7",
                                    access_token = NULL)

# Extracting data for 79 SNP(s) from 1 GWAS(s)
#save(outcome_dat, file = "./ieu/outcome_data_ieu-a-7.rdata")

# 统一数据
dat <- harmonise_data(exposure_dat, outcome_dat)

# Harmonising Body mass index || id:ieu-a-2 (ieu-a-2) 
# and Coronary heart disease || id:ieu-a-7 (ieu-a-7)

# 进行MR分析
res <- mr(dat)

# Analysing 'ieu-a-2' on 'ieu-a-7'

res

# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 3     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 4     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 5     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure                    method nsnp
# 1 Body mass index || id:ieu-a-2                  MR Egger   79
# 2 Body mass index || id:ieu-a-2           Weighted median   79
# 3 Body mass index || id:ieu-a-2 Inverse variance weighted   79
# 4 Body mass index || id:ieu-a-2               Simple mode   79
# 5 Body mass index || id:ieu-a-2             Weighted mode   79
# b         se         pval
# 1 0.5024935 0.14396056 8.012590e-04
# 2 0.3870065 0.07111020 5.258272e-08
# 3 0.4459091 0.05898302 4.032020e-14
# 4 0.3401554 0.15985362 3.650217e-02
# 5 0.3888249 0.10593531 4.412728e-04

结果中给出了：

method:使用的因果效应估计方法
nsnp:SNP的数量
b:暴露因素对结局的效应量
se:效应量的标准差
pval:MR分析的P值，小于0.05说明有因果关系。

还可以计算OR值和可信区间：

generate_odds_ratios(res)

# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 3     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 4     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 5     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure                    method nsnp
# 1 Body mass index || id:ieu-a-2                  MR Egger   79
# 2 Body mass index || id:ieu-a-2           Weighted median   79
# 3 Body mass index || id:ieu-a-2 Inverse variance weighted   79
# 4 Body mass index || id:ieu-a-2               Simple mode   79
# 5 Body mass index || id:ieu-a-2             Weighted mode   79
# b         se         pval      lo_ci     up_ci       or
# 1 0.5024935 0.14396056 8.012590e-04 0.22033081 0.7846562 1.652838
# 2 0.3870065 0.07111020 5.258272e-08 0.24763049 0.5263825 1.472566
# 3 0.4459091 0.05898302 4.032020e-14 0.33030238 0.5615158 1.561909
# 4 0.3401554 0.15985362 3.650217e-02 0.02684233 0.6534685 1.405166
# 5 0.3888249 0.10593531 4.412728e-04 0.18119168 0.5964581 1.475246
# or_lci95 or_uci95
# 1 1.246489 2.191653
# 2 1.280987 1.692797
# 3 1.391389 1.753328
# 4 1.027206 1.922196
# 5 1.198645 1.815676

or就是OR值，如果大于1就是风险因素，小于1就是保护因素，然后再结合pval，如果小于0.05就是有因果关系。

主要是看Inverse variance weighted的结果，其他方法作为辅助即可。

对MR分析结果进行可视化，这个图很明显使用ggplot2画的，所以可以使用相关语法进行修改：

library(ggplot2)
library(ggsci)
p1 <- mr_scatter_plot(res, dat = dat)

p1[[1]]+scale_color_lancet()+coord_fixed(0.5)

横坐标是SNP对暴露的效应，纵坐标是SNP对结局的效应，

这个图中每一个黑点表示一个SNP，横线表示SNP对暴露因素的效应量的可信区间（beta±se），竖线表示SNP对结局因素的效应量的可信区间（beta±se）。

5条直线是5种方法的拟合线（图中有2条线重合了），每种方法其实就是进行一种回归分析，所以有5条拟合线。每条直线的斜率（就是上面结果中的b）表示该方法估计的因果效应大小，如果斜率大于0，就表示是风险因素，如果小于0就表示保护因素。

支持多种因果效应估计方法，其中还给出了部分方法的参考文献，use_by_default会告诉你是不是默认的方法，heterogeneity_test会告诉你该方法能不能用于异质性检验：

# 18种方法
mr_method_list()

# 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
# 省略

如果你要更改方法也是非常简单：

mr(dat, method_list = c("mr_egger_regression", "mr_ivw"))

# Analysing 'ieu-a-2' on 'ieu-a-7'
# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure                    method nsnp
# 1 Body mass index || id:ieu-a-2                  MR Egger   79
# 2 Body mass index || id:ieu-a-2 Inverse variance weighted   79
# b         se        pval
# 1 0.5024935 0.14396056 8.01259e-04
# 2 0.4459091 0.05898302 4.03202e-14

敏感性分析

属于事后检验，目的是看看MR分析的结果是否靠谱。主要包括4个方面：

异质性检验
水平多效性
单个SNP分析
留一法分析

异质性检验就是检验不同的工具变量（SNP）之间是不是有非常大的差异，如果差异很大，可能会对结果的可靠性产生影响。

1行代码即可实现，某些方法既能进行因果效应分析，又可以进行异质性检验，比如MR-Egger和IVM：

mr_heterogeneity(dat)

# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure                    method        Q
# 1 Body mass index || id:ieu-a-2                  MR Egger 143.3046
# 2 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508
# Q_df       Q_pval
# 1   77 6.841585e-06
# 2   78 8.728420e-06

主要就是看结果中的Q_pval，也就是 Cochran's Q，小于0.05就是没有异质性。

也支持更改方法（支持的方法参考mr_method_list()）：

mr_heterogeneity(dat, method_list = c("mr_egger_regression", "mr_ivw"))

# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# 2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure                    method        Q
# 1 Body mass index || id:ieu-a-2                  MR Egger 143.3046
# 2 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508
# Q_df       Q_pval
# 1   77 6.841585e-06
# 2   78 8.728420e-06

下面是进行水平多效性检验，工具变量的多效性是指一个工具变量（也就是遗传变异，或者SNP）对多个表型（也就是暴露因素）都会产生影响的现象，这是违反独立性假设和排他性假设的。

这个水平多效性的检验是通过MR-Egger法的截距来判断的，如果截距（egger_intercept）在0附近（可以结合上面的图），且pval大于0.05说明没有水平多效性。

mr_pleiotropy_test(dat)

# id.exposure id.outcome                              outcome
# 1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
# exposure egger_intercept          se
# 1 Body mass index || id:ieu-a-2    -0.001719304 0.003985962
# pval
# 1 0.6674266

单个SNP分析，也就是对每一个SNP单独进行一次MR分析，目的就是看看是不是有某个SNP明显和其他的SNP不一样，识别离群值：

res_single <- mr_singlesnp(dat)
head(res_single)

# exposure
# 1 Body mass index || id:ieu-a-2
# 2 Body mass index || id:ieu-a-2
# 3 Body mass index || id:ieu-a-2
# 4 Body mass index || id:ieu-a-2
# 5 Body mass index || id:ieu-a-2
# 6 Body mass index || id:ieu-a-2
# outcome id.exposure id.outcome
# 1 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 2 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 3 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 4 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 5 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 6 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# samplesize        SNP          b        se          p
# 1     184305  rs1000940 -0.7901087 0.5347554 0.13953788
# 2     184305 rs10132280  0.5506335 0.4825339 0.25381623
# 3     184305  rs1016287  0.6178509 0.4637632 0.18277634
# 4     184305 rs10182181  0.5920712 0.3002913 0.04864884
# 5     184305 rs10733682  0.2415426 0.5066223 0.63352560
# 6     184305 rs10840100  0.7086893 0.4623010 0.12528548

这个结果也是可以获取OR值和可信区间的：

head(generate_odds_ratios(res_single))

# exposure
# 1 Body mass index || id:ieu-a-2
# 2 Body mass index || id:ieu-a-2
# 3 Body mass index || id:ieu-a-2
# 4 Body mass index || id:ieu-a-2
# 5 Body mass index || id:ieu-a-2
# 6 Body mass index || id:ieu-a-2
# outcome id.exposure id.outcome
# 1 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 2 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 3 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 4 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 5 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# 6 Coronary heart disease || id:ieu-a-7     ieu-a-2    ieu-a-7
# samplesize        SNP          b        se          p
# 1     184305  rs1000940 -0.7901087 0.5347554 0.13953788
# 2     184305 rs10132280  0.5506335 0.4825339 0.25381623
# 3     184305  rs1016287  0.6178509 0.4637632 0.18277634
# 4     184305 rs10182181  0.5920712 0.3002913 0.04864884
# 5     184305 rs10733682  0.2415426 0.5066223 0.63352560
# 6     184305 rs10840100  0.7086893 0.4623010 0.12528548
# lo_ci    up_ci        or  or_lci95 or_uci95
# 1 -1.838229348 0.258012 0.4537955 0.1590989 1.294354
# 2 -0.395133032 1.496400 1.7343514 0.6735904 4.465584
# 3 -0.291124912 1.526827 1.8549373 0.7474223 4.603545
# 4  0.003500324 1.180642 1.8077287 1.0035065 3.256464
# 5 -0.751437234 1.234522 1.2732116 0.4716881 3.436737
# 6 -0.197420583 1.614799 2.0313271 0.8208453 5.026879

结果可视化:

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

这个图是一个森林图，展示了单个SNP对CHD（冠心病）的影响，如果黑色水平线条完全在0的左边说明这个SNP是保护因素，完全在右边说明是危险因素，如果穿过了0说明没有意义。但是通常都是多个SNP一起作用的，所以单看一个没有太大说服力，直接看最低下的红色线条即可，这是汇总的效应，上图展示了两种方法的汇总效应，都是完全在0的右侧的。

单个SNP分析默认是Wald ratio检验，可以更改，比如使用固定效应meta分析法：

res_single <- mr_singlesnp(dat, single_method = "mr_meta_fixed")

总的效应默认是IVW法和MR-Egger法，也可以更改：

res_single <- mr_singlesnp(dat, all_method = "mr_two_sample_ml",
                           single_method = "mr_meta_fixed")

单个SNP分析还可以使用漏斗图进行可视化，主要是为了看看有没有明显的离群点，正常的话应该是左右基本对称分布的，如果某个点离得很远，可能需要把该SNP去掉：

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

留一法分析，即每次去掉一个SNP，对剩下的SNP进行MR分析，和单个SNP分析刚好相反。这一步也是看看有没有强影响点，正常来说去除某一个SNP之后，总体的效应应该变化不大（如下图所示），如果去掉某个SNP之后产生了很大的变化，说明这个SNP有问题，需要去掉。

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

这个进行分析的方法同样也是可以更改的哈。

Steiger方向性检验

MR-Steiger方向性检验，主要是看结局和暴露之间是不是有反向因果关系：

directionality_test(dat)

# r.exposure and/or r.outcome not present.
# Calculating approximate SNP-exposure and/or SNP-outcome correlations, 
# assuming all are quantitative traits. Please pre-calculate r.exposure 
# and/or r.outcome using get_r_from_lor() for any binary traits
# id.exposure id.outcome                      exposure
# 1     ieu-a-2    ieu-a-7 Body mass index || id:ieu-a-2
# outcome snp_r2.exposure
# 1 Coronary heart disease || id:ieu-a-7      0.01580819
# snp_r2.outcome correct_causal_direction  steiger_pval
# 1    0.001350485                     TRUE 1.748577e-207

该检验可以计算SNP对结局和暴露的 r^2^ ，也就是对结局的解释程度和对暴露的解释程度，理论上肯定是对暴露的解释程度更高（r^2^更大），因为SNP是和暴露因素强相关的。P值很小说明没有反向因果关系。

1-to-many forest plot

1-to-many森林图可以用来可视化多个暴露因素对一个结局的分析结果或者一个暴露因素对多个结局的分析结果。

对于小于50个结果的展示比较好，如果有多于50个结果建议分成两张图然后再拼起来。

例1

结局还是使用CHD，下面我们使用多个暴露因素：

# 6个暴露因素
exp_dat <- extract_instruments(
  outcomes = c("ieu-a-2", "ieu-a-100", "ieu-a-104", 
               "ieu-a-1", "ieu-a-72", "ieu-a-999"))
# 查看每个暴露因素有几个SNP
table(exp_dat$exposure)

# 提取结局相关的SNP
chd_out_dat <- extract_outcome_data(
  snps = exp_dat$SNP,
  outcomes = "ieu-a-7"
)

dat2 <- harmonise_data(
  exposure_dat = exp_dat, 
  outcome_dat = chd_out_dat
)
res2 <- mr(dat2)

subset_on_method函数对每一个暴露-结局组合只保留一个结果，如果工具变量大于1个就保留IVW法的结果，如果只有1个工具变量就保留Wald-ratio的结果。

# 只保留1个结果
res2 <- subset_on_method(res2) 
# 排序
res2 <- sort_1_to_many(res2, b = "b", sort_action = 4) # 按照效应递减排个序

# 把 Adiponectin || id:ieu-a-1 这种，去掉后面的id，好看一点
res2 <- split_exposure(res2) 

# 计算下权重，以便森林图的黑色方块可以有不同的大小
res2$weight <- 1/res2$se

# 计算下95%可信区间的上下限，方便画图时确定上下端点的值
min(exp(res2$b - 1.96*res2$se)) # 0.3233985，'lo' 
max(exp(res2$b + 1.96*res2$se)) # 2.489726，'up'

# 画图
forest_plot_1_to_many(
  res2,
  b = "b",
  se = "se",
  exponentiate = TRUE,
  ao_slc = FALSE,
  lo = 0.3,
  up = 2.5,
  TraitM = "exposure",
  col1_width = 2,
  by = NULL,
  trans = "log2",
  xlab = "OR for CHD per SD increase in risk factor (95% CI)",
  weight = "weight"
)

这张图就同时展示了6个暴露因素和结局，分别进行MR分析的结果，可以方便的看出每个暴露因素对结局的效应大小，最上面的效应最大，还展示了可信区间等。

还可以多添加几列信息，比如把P值和SNP的数量加到森林图中一起展示：

# 改下P值的显示方法
res2$pval<-formatC(res2$pval, format = "e", digits = 2)

forest_plot_1_to_many(
  res2,
  b = "b",
  se = "se",
  exponentiate = TRUE,
  ao_slc = FALSE,
  lo = 0.3,
  up = 2.5,
  TraitM = "exposure",
  by = NULL,
  trans = "log2",
  xlab = "OR for CHD per SD increase in risk factor (95% CI)",
  weight = "weight",
  subheading_size = 11,
  col1_title = "Risk factor",
  col1_width = 2.5,
  col_text_size = 4,
  addcols = c("nsnp", "pval"), # 把这两列加进去
  addcol_widths = c(1.0, 1.0) # 这两列的宽度
)

例2

上面是合并展示不同的暴露因素，每个MR分析只选了一种分析方法的结果进行展示。也可以根据暴露因素分层，每个暴露因素和结局的MR分析展示多个方法的结果。

res2 <- mr(dat2)
res2 <- split_exposure(res2) 

res2 <-
  sort_1_to_many(
    res2,
    group = "exposure", # 分层变量
    sort_action = 3, # 相同暴露因素的MR分析结果放在一起
    priority = "Inverse variance weighted", # 哪种方法放在最上面
    trait_m = "method"
  )

forest_plot_1_to_many(
  res2,
  b = "b",
  se = "se",
  exponentiate = TRUE,
  trans = "log2",
  ao_slc = FALSE,
  lo = 0.03,
  up = 22,
  col1_width = 2,
  by = "exposure", # 分层变量
  TraitM = "method",
  xlab = "OR for CHD per SD increase in risk factor (95% CI)",
  subheading_size = 12,
  col_text_size = 4
)

多变量MR（MVMR）

当SNP和多个暴露因素都有关系时，可以进行多变量MR分析（Multivariable MR）。大概的过程是这样的：

提取每个暴露因素的SNP
合并这些SNP
去除连锁不平衡的SNP
重新从以上数据中提取SNP
统一数据，确保effect allele一样
使用多变量MR方法进行分析

实际用代码实现是非常简单的过程。

比如我们有以下3个暴露因素，它们彼此之间是相关的：

HDL:ieu-a-299
LDL:ieu-a-300
total cholesterol:ieu-a-302

结局是冠心病：

coronary heart disease (CHD):ieu-a-7

id_exposure <- c("ieu-a-299", "ieu-a-300", "ieu-a-302")
id_outcome <- "ieu-a-7"

首先查找这几个暴露相关的SNP，这个函数会自定进行合并、去除LD等，1行代码就解决了：

# 提取SNP
exposure_datm <- mv_extract_exposures(id_exposure)

# Please look at vignettes for options on running this locally if you need to run many instances of this command.
# Clumping 1, 214 variants, using EUR population reference
# Removing 70 of 214 variants due to LD with other variants or absence from LD reference panel
# Extracting data for 144 SNP(s) from 3 GWAS(s)
# Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and LDL cholesterol || id:ieu-a-300 (ieu-a-300)
# Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and Triglycerides || id:ieu-a-302 (ieu-a-302)


head(exposure_datm[order(exposure_datm$SNP),],12)

#     SNP                        exposure id.exposure
# 1   rs10019888 HDL cholesterol || id:ieu-a-299   ieu-a-299
# 145 rs10019888 LDL cholesterol || id:ieu-a-300   ieu-a-300
# 289 rs10019888   Triglycerides || id:ieu-a-302   ieu-a-302
# 2     rs103294 HDL cholesterol || id:ieu-a-299   ieu-a-299
# 146   rs103294 LDL cholesterol || id:ieu-a-300   ieu-a-300
# 290   rs103294   Triglycerides || id:ieu-a-302   ieu-a-302
# 3   rs10468017 HDL cholesterol || id:ieu-a-299   ieu-a-299
# 147 rs10468017 LDL cholesterol || id:ieu-a-300   ieu-a-300
# 291 rs10468017   Triglycerides || id:ieu-a-302   ieu-a-302
# 4    rs1047891 HDL cholesterol || id:ieu-a-299   ieu-a-299
# 148  rs1047891 LDL cholesterol || id:ieu-a-300   ieu-a-300
# 292  rs1047891   Triglycerides || id:ieu-a-302   ieu-a-302
#     effect_allele.exposure other_allele.exposure eaf.exposure
# 1                        G                     A       0.1636
# 145                      G                     A       0.1636
# 289                      G                     A       0.1636
# 2                        T                     C       0.1860
# 146                      T                     C       0.1860
# 290                      T                     C       0.1860
# 3                        T                     C       0.2757
# 147                      T                     C       0.2757
# 291                      T                     C       0.2757
# 4                        A                     C       0.3021
# 148                      A                     C       0.3021
# 292                      A                     C       0.3021
#     beta.exposure se.exposure pval.exposure
# 1         -0.0270      0.0046   4.90095e-08
# 145        0.0182      0.0050   3.23199e-04
# 289        0.0228      0.0045   2.28302e-06
# 2          0.0523      0.0044   3.99485e-30
# 146        0.0073      0.0047   1.22500e-01
# 290       -0.0021      0.0042   7.52101e-01
# 3          0.1179      0.0038  1.21060e-188
# 147        0.0020      0.0042   7.46799e-01
# 291        0.0379      0.0039   7.55962e-21
# 4         -0.0269      0.0039   8.72991e-10
# 148        0.0079      0.0042   1.42200e-01
# 292        0.0000      0.0038   8.60200e-01

可以看到每个SNP对应3种暴露。

然后提取和结局相关的SNP：

outcome_datm <- extract_outcome_data(exposure_datm$SNP, id_outcome)

# Extracting data for 144 SNP(s) from 1 GWAS(s)
# Finding proxies for 1 SNPs in outcome ieu-a-7
# Extracting data for 1 SNP(s) from 1 GWAS(s)

统一数据：

mvdat <- mv_harmonise_data(exposure_datm, outcome_datm)
# Harmonising HDL cholesterol || id:ieu-a-299 (ieu-a-299) and Coronary 
# heart disease || id:ieu-a-7 (ieu-a-7)

进行多变量MR分析：

resm <- mv_multiple(mvdat)
resm

# $result
# id.exposure                        exposure id.outcome
# 1   ieu-a-299 HDL cholesterol || id:ieu-a-299    ieu-a-7
# 2   ieu-a-300 LDL cholesterol || id:ieu-a-300    ieu-a-7
# 3   ieu-a-302   Triglycerides || id:ieu-a-302    ieu-a-7
# outcome nsnp           b         se
# 1 Coronary heart disease || id:ieu-a-7   79 -0.08919724 0.05970552
# 2 Coronary heart disease || id:ieu-a-7   68  0.37853543 0.04976846
# 3 Coronary heart disease || id:ieu-a-7   42  0.13584165 0.06738291
# pval
# 1 1.351879e-01
# 2 2.828614e-14
# 3 4.380354e-02

结果说明ieu-a-299，也就是HDL对冠心病是没有因果关系的。

如果你单独对每个暴露因素做MR分析，结果可能和多变量MR的结果不同。

MR-MoE

使用机器学习方法确定最佳的MR方法。MR-MoE:Using a mixture of experts machine learning approach

可以使用随机森林算法评估11种MR方法的准确性。

相关的参考文献（详细的方法学介绍）可以通过链接查看：https://www.biorxiv.org/content/10.1101/173682v2

使用方法是很简单的，首先需要下载官方制作的训练数据，下载链接：https://www.dropbox.com/s/5la7y38od95swcf/rf.rdata?dl=0

然后通过以下代码即可实现使用随机森林确定最佳的MR方法：

# 前面的步骤都是一样的
exposure_dat <- extract_instruments("ieu-a-2")
outcome_dat <- extract_outcome_data(exposure_dat$SNP, "ieu-a-7")
dat <- harmonise_data(exposure_dat, outcome_dat)

# 加载下载的数据
load("./ieu/rf.rdata")

# 整合所有MR方法的结果
res1 <- mr_wrapper(dat)

# Performing MR analysis of 'ieu-a-2' on 'ieu-a-7'
# Body mass index || id:ieu-a-2 - Coronary heart disease || id:ieu-a-7
# Body mass index || id:ieu-a-2 - Coronary heart disease || id:ieu-a-7
# Body mass index || id:ieu-a-2 - Coronary heart disease || id:ieu-a-7
# Body mass index || id:ieu-a-2 - Coronary heart disease || id:ieu-a-7
 

library(dplyr)
# MR-MoE：预测每种方法的性能
res_moe <- mr_moe(res1, rf)

结果查看：

res_moe$`ieu-a-2.ieu-a-7`$estimates

# A tibble: 44 × 14
   id.exposure id.outcome method     nsnp     b     se ci_low ci_upp
   <chr>       <chr>      <chr>     <int> <dbl>  <dbl>  <dbl>  <dbl>
 1 ieu-a-2     ieu-a-7    FE IVW       69 0.464 0.0415  0.383  0.546
 2 ieu-a-2     ieu-a-7    RE IVW       69 0.464 0.0415  0.383  0.546
 3 ieu-a-2     ieu-a-7    Simple m…    69 0.439 0.0465  0.348  0.530
 4 ieu-a-2     ieu-a-7    FE IVW       79 0.446 0.0439  0.330  0.561
 5 ieu-a-2     ieu-a-7    RE IVW       69 0.464 0.0415  0.383  0.546
 6 ieu-a-2     ieu-a-7    Simple m…    79 0.393 0.0683  0.259  0.527
 7 ieu-a-2     ieu-a-7    Simple m…    69 0.397 0.0687  0.262  0.532
 8 ieu-a-2     ieu-a-7    FE IVW       69 0.464 0.0415  0.383  0.546
 9 ieu-a-2     ieu-a-7    Weighted…    69 0.383 0.0747  0.236  0.529
10 ieu-a-2     ieu-a-7    FE IVW       79 0.446 0.0439  0.330  0.561
# ℹ 34 more rows
# ℹ 6 more variables: pval <dbl>, steiger_filtered <lgl>,
#   outlier_filtered <lgl>, selection <chr>, method2 <chr>,
#   MOE <dbl>
# ℹ Use `print(n = ...)` to see more rows

这个结果res_moe是一个列表，主要包含以下几个元素：

estimates：每种MR方法的结果
heterogeneity：异质性检验的结果
directional_pleiotropy：egger intercepts的结果
info：用于产生MOE的指标

主要是看estimates中的MOE列，这是每个方法的AUC值，越大说明这个方法越好。

整合所有结果

就是把所有结果放到一个数据框中，感觉用处不大...

res <- mr(dat)
het <- mr_heterogeneity(dat)
plt <- mr_pleiotropy_test(dat)
sin <- mr_singlesnp(dat)

# 合并所有结果
all_res <-
  combine_all_mrresults(
    res,
    het,
    plt,
    sin,
    ao_slc = TRUE,
    Exp = TRUE,
    split.exposure = FALSE,
    split.outcome = TRUE
  )
head(all_res[, c(
  "Method",
  "outcome",
  "exposure",
  "nsnp",
  "b",
  "se",
  "pval",
  "intercept",
  "intercept_se",
  "intercept_pval",
  "Q",
  "Q_df",
  "Q_pval",
  "consortium",
  "ncase",
  "ncontrol",
  "pmid",
  "population"
)])

# 1 Inverse variance weighted Coronary heart disease
# 2                  MR Egger Coronary heart disease
# 3               Simple mode Coronary heart disease
# 4                Wald ratio Coronary heart disease
# 5                Wald ratio Coronary heart disease
# 6                Wald ratio Coronary heart disease
#                       exposure nsnp          b         se
# 1 Body mass index || id:ieu-a-2   79  0.4459091 0.05898302
# 2 Body mass index || id:ieu-a-2   79  0.5024935 0.14396056
# 3 Body mass index || id:ieu-a-2   79  0.3401554 0.15878926
# 4 Body mass index || id:ieu-a-2    1  1.1494985 0.37045723
# 5 Body mass index || id:ieu-a-2    1 -0.7901087 0.53475543
# 6 Body mass index || id:ieu-a-2    1  0.5506335 0.48253394
#           pval    intercept intercept_se intercept_pval        Q
# 1 4.032020e-14           NA           NA             NA 143.6508
# 2 8.012590e-04 -0.001719304  0.003985962      0.6674266 143.3046
# 3 3.529868e-02           NA           NA             NA       NA
# 4 1.916225e-03           NA           NA             NA       NA
# 5 1.395379e-01           NA           NA             NA       NA
# 6 2.538162e-01           NA           NA             NA       NA
#              Q_df       Q_pval        consortium ncase ncontrol     pmid
# 1   78 8.728420e-06 CARDIoGRAMplusC4D 60801   123504 26343387
# 2   77 6.841585e-06 CARDIoGRAMplusC4D 60801   123504 26343387
# 3   NA           NA CARDIoGRAMplusC4D 60801   123504 26343387
# 4   NA           NA CARDIoGRAMplusC4D 60801   123504 26343387
# 5   NA           NA CARDIoGRAMplusC4D 60801   123504 26343387
# 6   NA           NA CARDIoGRAMplusC4D 60801   123504 26343387
#    population
# 1      Mixed
# 2      Mixed
# 3      Mixed
# 4      Mixed
# 5      Mixed
# 6      Mixed

OVER！

后面的推文就是介绍如何从各大数据库下载数据到本地进行MR分析了。

联系我们，关注我们
免费QQ交流群1：613637742
免费QQ交流群2：608720452
公众号消息界面关于作者获取联系方式
知乎、CSDN、简书同名账号
哔哩哔哩：阿越就是我

ixxmu / mp_duty