孟德尔随机化(Mendelian randomization
,MR)是一种利用基因变异作为工具变量来评估暴露与结果之间因果关系的统计方法。
它基于这样的原理:基因变异是在出生前就随机分配给个体的,类似于在随机对照试验中随机分配治疗,因此可以帮助区分因果关系和简单相关性。孟德尔随机化通常用于观察性数据,以确定一个特定的**生物标志物、行为**或**其他暴露**是否真正地影响了健康结果,而不是仅仅与之相关。通过这种方法,研究者可以减少混杂因素的影响,避免了传统观察性研究中常见的一些偏差。
可以看一下前一篇的文献分享:
【多线程核糖体】公众号后台回复MR1领取示例文件与代码
这里复现的是教育年限
对EC
的影响
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)
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)
选择加权还是不加权取决于数据特性
和分析目的
,当所有基因-暴露估计的可变性(标准误)大致相同时,不加权的方法可能更为适用,反之则使用加权。
教大家一个更简单的方法:**加权和不加权都做,选结果好的(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)
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
最后,示例数据仅供参考,笔者发现这个分析结果并不如意,权作教程示例数据使用,无实际意义!
有需要的小伙伴可以直接替换数据学起来吧~
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。