孟德尔随机化代码实现

安装并载入软件包：

install.packages("devtools")
devtools::install_github("MRCIEU/TwoSampleMR")
install_github("MRCIEU/TwoSampleMR")

library("devtools")
library(TwoSampleMR)

运行TSMR

提取暴露因素的工具变量

options(ieugwasr_api = 'gwas-api.mrcieu.ac.uk/')#让IEU网站即时出结果
exposure_dat <- extract_instruments(outcomes = 'ebi-a-GCST90029019')
head(exposure_dat)

如果想要调整clump参数,可以把R方调整成0.01

exposure_dat <- extract_instruments(outcomes = 'prot-c-2992_59_2',
                                         clump = TRUE, r2 = 0.01,
                                         kb = 5000)

如果想要调整P值，r2<0.001,kb=10000用于去除连锁不平衡

exposure_dat <- extract_instruments(outcomes = 'ebi-a-GCST006956',
                                               p1 = 5e-08,
                                               clump = TRUE,r2 = 0.001,kb = 10000, 
                                               )

将工具SNP作为一个新的csv文件输出

write.csv(exposure_dat,"exposure_dat.csv", row.names=F)

在胃癌ebi-a-GCST90018849这个数据中找到与上面SNP匹配的数据，这样就生成了暴露数据

outcome_dat<-extract_outcome_data(snps=exposure_dat$SNP,
                                  outcomes='ebi-a-GCST90018849', 
                                  proxies = FALSE,
                                maf_threshold = 0.01)

效应等位与效应量保持统一

dat <- harmonise_data(
  exposure_dat = exposure_dat,
  outcome_dat = outcome_dat
)

进行MR分析，b是效应值，se是标准误，pval是P值，mr默认使用五种方法（ MR Egger，Weighted median，Inverse variance weighted，Simple mode ，Weighted mode ）最重要的是Inverse variance weighted

res <- mr(dat)
res

异质性检验 mr_heterogeneity

mr_heterogeneity(dat)

水平多效性检验，如果变量工具不通过暴露影响结果（p值大于0.05是正确的），就违反了孟德尔的假设，就是存在多水平效应。

mr_pleiotropy_test(dat)

Leave-one-out analysis,留一图分析是指逐步剔除SNP后观察剩余的稳定性，理想的是剔除后变化不大，这和我们的meta分析剔除法很相似。

res_loo <- mr_leaveoneout(dat)
mr_leaveoneout_plot(res_loo)

可视化，做出散点图，看相关趋势。

p1 <- mr_scatter_plot(res, dat)
p1

森林图，森林图主要是看总效应

res_single <- mr_singlesnp(dat)
mr_forest_plot(res_single)