也是一个计算孟德尔随机化的常用包
CRAN - Package MendelianRandomization (r-project.org)
https://academic.oup.com/ije/article/46/6/1734/3112150
MendelianRandomization: Mendelian Randomization Package (r-project.org)
里面还有多种方法，注意带有mr_m为多变量的，不带有的为单变量

image.png

path.noproxy <- system.file("extdata", "vitD_snps_PhenoScanner.csv",
                            package = "MendelianRandomization")
path.proxies <- system.file("extdata", "vitD_snps_PhenoScanner_proxies.csv",
                            package = "MendelianRandomization")
# these two files from PhenoScanner are provided
# as part of the MendelianRandomization package
extract.pheno.csv(
  exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
  outcome = "Tanner stage", pmidO = 24770850, ancestryO = "European",
  file = path.noproxy)
extract.pheno.csv(
  exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
  outcome = "Tanner stage", pmidO = 24770850, ancestryO = "European",
  rsq.proxy = 0.6, file = path.proxies)
extract.pheno.csv(
  exposure = "log(eGFR creatinine)", pmidE = 26831199, ancestryE = "European",
  outcome = "Asthma", pmidO = 20860503, ancestryO = "European",
  rsq.proxy = 0.6, file = path.proxies)

单变量方法

IVW方法

library(MendelianRandomization) #加载R包
MRInputObject <- mr_input(bx = ldlc,bxse= ldlcse,by = chdlodds,byse = chdloddsse) #指定输入文件
MRInputObject 
### IVW方法估算结果
IVWObject1 <- mr_ivw(MRInputObject,model= "default",robust = FALSE,penalized = FALSE,correl = FALSE,weights ="simple", psi = 0,distribution = "normal",alpha = 0.05)
IVWObject1 
## 使用稳健回归，并且对SNP中的outliers进行“惩罚”
IVWObject2 <- mr_ivw(MRInputObject,model= "default",robust = TRUE,penalized = TRUE,correl = FALSE,weights ="simple", psi = 0,distribution = "normal",alpha = 0.05)
IVWObject2 # 查看结果

median-based方法

## median-based方法，可以作为对IVW方法的补充
# MRInputObject <- mr_input(bx = ldlc,bxse= ldlcse,by = chdlodds,byse = chdloddsse) #指定输入文件
WeightedMedianObject1 <-mr_median(MRInputObject,weighting = "weighted",distribution ="normal",alpha = 0.05,iterations = 10000,seed = 314159265)

WeightedMedianObject1 #结果显著，并且LDL的升高可以增加CHD的发病风险

### penalized“加权法
WeightedMedianObject2 <-mr_median(MRInputObject,weighting = "penalized",distribution ="normal",alpha = 0.05,iterations = 10000,seed = 314159265)
WeightedMedianObject2
### 不采用加权法来计算一下结果
WeightedMedianObject3 <-mr_median(MRInputObject,weighting = "simple",distribution ="normal",alpha = 0.05,iterations = 10000,seed = 314159265)
WeightedMedianObject3

### 加权模型下增加迭代次数（iterations）
WeightedMedianObject4 <-mr_median(MRInputObject,weighting = "weighted",distribution ="normal",alpha = 0.05,iterations = 100000,seed = 314159265) #修改iterations参数为100000
WeightedMedianObject4
# 在同一种模型（比如加权模型）之下，增加bootstrap的迭代次数，可以减少误差，使得结果更加准确，但是增加迭代次数之后，计算量会显著增大，计算时间会相应延长

MR-Egger法

### MR-Egger法
library(MendelianRandomization) #加载R包
MRInputObject <- mr_input(bx = ldlc,bxse= ldlcse,by = chdlodds,byse = chdloddsse) #指定输入文件
EggerObject1 <-mr_egger(MRInputObject,robust = FALSE,penalized = FALSE,correl =FALSE,distribution = "normal",alpha = 0.05)
EggerObject1 

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

#### 使用稳健回归后，MR估计值的误差变小了，但显著性并未改变。这也充分说明，LDL的升高可以增加CHD的发病风险。

Maximum likelihood方法（极大似然估计法)

#Maximum likelihood方法（极大似然估计法)

# 相对于IVW方法，极大似然估计法也有着自身的两点优势：（1）它充分考虑了SNP-exposure关联的不确定性，这一点在IVW的简单加权中是被忽略的；（2）它考虑了双样本MR研究中样本重叠的情况


### 通过参数psi来控制，在样本完全不重叠的情况下（传统的独立双样本MR研究），psi为0；
###如果样本完全重叠的话，psi等于观察性研究中暴露和结局的相关性（相关系数）
MaxLikObject1 <-mr_maxlik(MRInputObject,model = "default",correl = FALSE,psi =0,distribution = "normal",alpha = 0.05)

MaxLikObject1 
MaxLikObject2 <-mr_maxlik(MRInputObject,model = "default",correl = FALSE,psi = 0.3,distribution= "normal",alpha = 0.05)
MaxLikObject2
MaxLikObject3 <-mr_maxlik(MRInputObject,model = "default",correl = FALSE,psi = 0.6,distribution= "normal",alpha = 0.05)
MaxLikObject3

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

多样本

MRMVInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig),
                              bxse = cbind(ldlcse, hdlcse, trigse),
                              by = chdlodds, byse = chdloddsse) #MVMR的input格式会和单变量的有所不同
MRMVInputObject

MRMVObject1 <- mr_mvivw(MRMVInputObject,
                        model = "default",
                        correl = FALSE,
                        distribution = "normal",
                        alpha = 0.05)
MRMVObject2 <-  mr_mvegger(MRMVInputObject, 
                           orientate = 1,
                           correl = FALSE,
                           distribution = "normal", 
                           alpha = 0.05)

绘图

结果取决于传递给 mr_plot 的对象类型。当对象是 MRInput对象，该函数使用 plot 命令（如果 interactive 设置为 FALSE）或 plotly语法（如果 interactive 设置为 TRUE）来绘制关联估计值。当。。。的时候object 是一个 MRMVInput 对象，功能类似，除了我们绘制估计的关联结果在 y 轴上，关联的拟合值与来自x轴上的逆方差加权方法。如果 interactive 设置为 FALSE，则为静态图被生产。通过将标签设置为 TRUE，基因变体的名称出现在点上方。这会产生一个视觉上不太吸引人的图表，但更容易识别个人遗传变异。如果 interactive 设置为 TRUE，则绘图是交互式的，用户可以悬停在各个点上查看相关遗传变异的名称及其关联估计。当对象是 MRAll 对象时，该函数生成一个 ggplot 来比较因果估计通过不同的方法提出。

mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
        line="egger", orientate = TRUE)

image.png

mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
        line="ivw", interactive=FALSE) # produces a static graph

image.png

mr_plot(mr_allmethods(mr_input(bx = ldlc, bxse = ldlcse,
                               by = chdlodds, byse = chdloddsse), method="all", iterations = 50))

image.png

https://zhuanlan.zhihu.com/p/439009634