MendelianRandomization | 孟德尔随机化神包更新啦！~（一）（小试牛刀）-编程知识

1写在前面

今天发现MendelianRandomization包更新v0.9了。😜

其实也算不上更新。🫠

跟大家一起分享一下这个包做MR的用法吧。🤩

还有一个包就是TwoSampleMR，大家有兴趣可以去学一下。😅

2用到的包

rm(list = ls())
# install.packages("MendelianRandomization")
library(MendelianRandomization)
library(tidyverse)

3示例数据

这里需要说明一下哦，不是所有的参数都需要填的，这个需要根据你后面用到的方法。😘

比如说一些方法不需要指定bxse。😜

mr_ivw函数是将BXSE设置为零的情况下运行的。😀

示例一

betaX是来自每个SNP的暴露因素的回归分析的β系数，betaXse是标准误。😘

betaY是每个SNP的结局时间的回归分析的β系数数，betaYse是标准误。🥳

MRInputObject <- mr_input(bx = ldlc,
                          bxse = ldlcse,
                          by = chdlodds,
                          byse = chdloddsse)
MRInputObject

示例二

如果SNP是存在存在相关性，这个时候可以通过corr = 加入哦。😯

MRInputObject.cor <- mr_input(bx = calcium,
                              bxse = calciumse,
                              by = fastgluc,
                              byse = fastglucse,
                              corr = calc.rho)
MRInputObject.cor

这里的两个data是出自下面2个原文，大家可以去读一下,试试复现：👇

1️⃣ Burgess S, Scott RA, Timpson NJ, Davey Smith G, Thompson SG; EPIC- InterAct Consortium. Using published data in Mendelian randomization: a blueprint for efficient identification of causal risk factors. Eur J Epidemiol. 2015 Jul;30(7):543-52. doi: 10.1007/s10654-015-0011-zIF: 13.6 Q1 . Epub 2015 Mar 15. PMID: 25773750 IF: 13.6 Q1 B1; PMCID: PMC4516908 IF: 13.6 Q1 B1.

2️⃣ Waterworth DM, Ricketts SL, Song K, Chen L, Zhao JH, Ripatti S, Aulchenko YS, Zhang W, Yuan X, Lim N, Luan J, Ashford S, Wheeler E, Young EH, Hadley D, Thompson JR, Braund PS, Johnson T, Struchalin M, Surakka I, Luben R, Khaw KT, Rodwell SA, Loos RJ, Boekholdt SM, Inouye M, Deloukas P, Elliott P, Schlessinger D, Sanna S, Scuteri A, Jackson A, Mohlke KL, Tuomilehto J, Roberts R, Stewart A, Kesäniemi YA, Mahley RW, Grundy SM; Wellcome Trust Case Control Consortium; McArdle W, Cardon L, Waeber G, Vollenweider P, Chambers JC, Boehnke M, Abecasis GR, Salomaa V, Järvelin MR, Ruokonen A, Barroso I, Epstein SE, Hakonarson HH, Rader DJ, Reilly MP, Witteman JC, Hall AS, Samani NJ, Strachan DP, Barter P, van Duijn CM, Kooner JS, Peltonen L, Wareham NJ, McPherson R, Mooser V, Sandhu MS. Genetic variants influencing circulating lipid levels and risk of coronary artery disease. Arterioscler Thromb Vasc Biol. 2010 Nov;30(11):2264-76. doi: 10.1161/ATVBAHA.109.201020 IF: 8.7 Q1 B1. Epub 2010 Sep 23. PMID: 20864672 IF: 8.7 Q1 B1; PMCID: PMC3891568 IF: 8.7 Q1 B1.

4方法一 Inverse-variance weighted method

在MR的因果推断中，我们常用四种方法，inverse-variance weighted method, median-based method，MR-Egger method，还有Maximum likelihood method。😘

我们注意逐一介绍下。🔽

IVWObject <- mr_ivw(MRInputObject,
                    model = "default",
                    robust = FALSE,
                    penalized = FALSE,
                    correl = FALSE,
                    weights = "simple",
                    psi = 0,
                    distribution = "normal",
                    alpha = 0.05)

IVWObject <- mr_ivw(mr_input(bx = ldlc, 
                             bxse = ldlcse,
                             by = chdlodds, 
                             byse = chdloddsse)
                    )

IVWObject

IVWObject.correl <- mr_ivw(MRInputObject.cor,
                          model = "default",
                          correl = TRUE,
                          distribution = "normal",
                          alpha = 0.05)
                          IVWObject.correl <- mr_ivw(mr_input(bx = calcium, bxse = calciumse,
                          by = fastgluc, byse = fastglucse, corr = calc.rho)
                          )
IVWObject.correl

5方法二 Median-based method

WeightedMedianObject <- mr_median(MRInputObject,
                                  weighting = "weighted",
                                  distribution = "normal",
                                  alpha = 0.05,
                                  iterations = 10000,
                                  seed = 314159265)

WeightedMedianObject <- mr_median(mr_input(bx = ldlc, 
                                           bxse = ldlcse,
                                           by = chdlodds, 
                                           byse = chdloddsse))

WeightedMedianObject

SimpleMedianObject <- mr_median(mr_input(bx = ldlc, 
                                         bxse = ldlcse,
                                         by = chdlodds, 
                                         byse = chdloddsse), 
                                weighting = "simple")

SimpleMedianObject

6方法三 MR-Egger method

EggerObject <- mr_egger(MRInputObject,
                        robust = FALSE,
                        penalized = FALSE,
                        correl = FALSE,
                        distribution = "normal",
                        alpha = 0.05)

EggerObject <- mr_egger(mr_input(bx = ldlc, 
                                 bxse = ldlcse,
                                 by = chdlodds, 
                                 byse = chdloddsse)
                        )

EggerObject

EggerObject.corr <- mr_egger(MRInputObject.cor,
                            correl = TRUE,
                            distribution = "normal",
                            alpha = 0.05)

EggerObject.corr <- mr_egger(mr_input(bx = calcium, 
                                      bxse = calciumse,
                                      by = fastgluc, 
                                      byse = fastglucse, 
                                      corr = calc.rho))

EggerObject.corr

7方法四 Maximum likelihood method

MaxLikObject <- mr_maxlik(MRInputObject,
                          model = "default",
                          correl = FALSE,
                          psi = 0,
                          distribution = "normal",
                          alpha = 0.05)

MaxLikObject <- mr_maxlik(mr_input(bx = ldlc, bxse = ldlcse,by = chdlodds, byse = chdloddsse))

MaxLikObject

MaxLikObject.corr <- mr_maxlik(mr_input(bx = calcium, bxse = calciumse,by = fastgluc, byse = fastglucse, corr = calc.rho))

MaxLikObject.corr

8简单可视化

这里就以IVW为例哦。🙃

mr_plot(MRInputObject,
        error = T, 
        orientate = F, 
        line = "ivw", # "ivw", "egger"
        interactive = F,
        labels = T
        )