全代码分享｜孟德尔随机化怎么做？MR一套标准流程

原创

生信初学者

发布于 2024-03-29 09:54:33

1.3K0

发布于 2024-03-29 09:54:33

文章被收录于专栏：备份笔记

1.前言

孟德尔随机化(Mendelian randomization，MR)是一种利用基因变异作为工具变量来评估暴露与结果之间因果关系的统计方法。

它基于这样的原理：基因变异是在出生前就随机分配给个体的，类似于在随机对照试验中随机分配治疗，因此可以帮助区分因果关系和简单相关性。孟德尔随机化通常用于观察性数据，以确定一个特定的**生物标志物、行为**或**其他暴露**是否真正地影响了健康结果，而不是仅仅与之相关。通过这种方法，研究者可以减少混杂因素的影响，避免了传统观察性研究中常见的一些偏差。

可以看一下前一篇的文献分享：

【多线程核糖体】公众号后台回复MR1领取示例文件与代码

2.demo

这里复现的是教育年限对EC的影响

2.1 加载R包

library(devtools)

library(tidyr)

# devtools::install\_github("MRCIEU/MRInstruments")

library(MRInstruments)

library(MendelianRandomization) 

library(TwoSampleMR) 

#install.packages("simex")

library(simex) 

#install\_github("rondolab/MR-PRESSO")

library(MRPRESSO)

2.2 主要MR分析

mydata <- read.table("01mydata\_edu\_cancer.txt", header=T, sep="\t", check.names = F, stringsAsFactors = F)

heterogeneity <- mr\_heterogeneity(mydata)

heterogeneity

# id.exposure id.outcome outcome          exposure                    method        Q Q\_df    Q\_pval

# 1  ieu-a-1001     lAx3ZW outcome  || id:ieu-a-1001                  MR Egger 347.7875  362 0.6951200

# 2  ieu-a-1001     lAx3ZW outcome  || id:ieu-a-1001 Inverse variance weighted 347.8263  363 0.7075905



### 二选一

# 2.1 如果I2<25%或Q\_P>0.05，使用固定效应的逆方差加权法（IVW）

mr\_results <- mr(mydata, method\_list=c('mr\_ivw\_fe')) 

mr\_results

generate\_odds\_ratios(mr\_results)

# 2.2 如果I2>25%且Q\_P<0.05，使用随机效应的逆方差加权法（IVW）

mr\_results <- mr(mydata, method\_list=c('mr\_ivw\_mre')) 

mr\_results

generate\_odds\_ratios(mr\_results)

2.3 MR补充分析、多态性、验证

选择加权还是不加权取决于数据特性和分析目的，当所有基因-暴露估计的可变性（标准误）大致相同时，不加权的方法可能更为适用，反之则使用加权。

教大家一个更简单的方法：**加权和不加权都做，选结果好的（P值更显著的），如果都显著则选择加权的，因为它还照顾到每个点的不确定性，更为全面**

mr\_results1 <- mr(mydata, method\_list=c('mr\_two\_sample\_ml', 'mr\_egger\_regression',

                                        "mr\_simple\_median", "mr\_weighted\_median", "mr\_penalised\_weighted\_median",

                                        "mr\_simple\_mode", "mr\_weighted\_mode")) 

mr\_results1



# 4.1 MR-Egger截距测试

mr\_pleiotropy\_test(mydata)

mr\_egger(mr\_input(bx = mydata$beta.exposure, bxse = mydata$se.exposure, by = mydata$beta.outcome, byse = mydata$se.outcome)) 



# 4.2 MR多态性残差和异常值（MR-PRESSO）测试

# 运行MR-PRESSO全局方法

mr\_presso(BetaOutcome = "beta.outcome", BetaExposure = "beta.exposure", SdOutcome = "se.outcome", SdExposure = "se.exposure", 

          OUTLIERtest = TRUE, DISTORTIONtest = TRUE, data = mydata, NbDistribution = 1000, SignifThreshold = 0.05, seed = 123456) 



# 5.1 IVW radial: mr\_ivw\_radial

mr(mydata, method\_list=c("mr\_ivw\_radial")) 

# 5.2 Egger-SIMEXS

mr\_egger(mr\_input(bx = mydata$beta.exposure, bxse = mydata$se.exposure, by = mydata$beta.outcome, byse = mydata$se.outcome))

BetaXG <- mydata$beta.exposure

BetaYG <- mydata$beta.outcome

seBetaXG <- mydata$se.exposure

seBetaYG <- mydata$se.outcome

BYG <- BetaYG \* sign(BetaXG) # 预处理步骤，确保所有基因-暴露估计都是正的

BXG <- abs(BetaXG)  



# MR-Egger回归（未加权）

Fit1 <- lm(BYG ~ BXG, x=TRUE, y=TRUE)

mod.sim1 <- simex(Fit1, B=1000, measurement.error = seBetaXG, SIMEXvariable="BXG", fitting.method ="quad", asymptotic="FALSE") 

summary(mod.sim1)

# MR-Egger回归（加权）

Fit2 <- lm(BYG ~ BXG, weights=1/seBetaYG^2, x=TRUE, y=TRUE) 

mod.sim2 <- simex(Fit2, B=1000, measurement.error = seBetaXG, SIMEXvariable="BXG", fitting.method ="quad", asymptotic="FALSE") 

summary(mod.sim2)



# IVW不加权

Fit3 <- lm(BYG ~ -1 + BXG, x=TRUE, y=TRUE) 

mod.sim3 <- simex(Fit3, B=1000, measurement.error = seBetaXG, SIMEXvariable="BXG", fitting.method="quad", asymptotic="FALSE") 

summary(mod.sim3)

# IVW加权

Fit4 <- lm(BYG ~ -1 + BXG, weights=1/seBetaYG^2, x=TRUE, y=TRUE) 

mod.sim4 <- simex(Fit4, B=1000, measurement.error = seBetaXG, SIMEXvariable="BXG", fitting.method="quad", asymptotic="FALSE") 

summary(mod.sim4)

2.4 结果可视化

library(ggplot2)

# 散点图

plot1 <- mr\_scatter\_plot(mr\_results1, mydata) 

plot1



# 单个SNP森林图

res\_single <- mr\_singlesnp(

  mydata,

  parameters = default\_parameters(),

  single\_method = "mr\_wald\_ratio",

  all\_method = c("mr\_ivw", 'mr\_two\_sample\_ml', "mr\_weighted\_median")

) 

significant\_snps <- res\_single[res\_single$p < 0.05, ] # 只保留P<0.05

plot2 <- mr\_forest\_plot(significant\_snps)

plot2



# 留一法敏感性测试

single <- mr\_leaveoneout(mydata)

top\_30 <- single[order(single$p), ][1:30, ] # 只保留P最小的前30个

plot3 <- mr\_leaveoneout\_plot(top\_30)

plot3



# 漏斗图

plot4 <- mr\_funnel\_plot(res\_single)

plot4