孟de尔R代码

setwd("D:\\MDR database")

library(TwoSampleMR)

CRP_exp_dat <- extract_instruments(outcomes = 'met-d-Remnant_C')

#参数默认为P=5*10-8，r方=0.001，kb=10000

#查看暴露SNP基本情况

dim(CRP_exp_dat)  #查看筛选出多少个SNP

head(CRP_exp_dat) #查看前五行信息

write.csv(CRP_exp_dat, file="D:\\MDR database\\CRP.csv")

##结局指标：在暴露因素的SNP里面找结局指标的（从56--52个）

CHD_out <- extract_outcome_data(

  

    snps=CRP_exp_dat$SNP, #根据暴露SNP编号在结局中查找对用SNP

    outcomes="ieu-a-7", #查找数据库ID号为ebi-a-GCST006906的SNP

    proxies = TRUE,#在结局中找不到的SNP是否用1000组的代替，默认TRUE，R方0.8

    maf_threshold = 0.01,#最小等位基因频率

    access_token = NULL#是否通过谷歌访问，国内选NULL

  )

#查看暴露SNP基本情况，发现少了几个，

#一般可能不符合条件，或者再结局中没有找到

dim(CHD_out) #查看筛选出多少个SNP

head(CHD_out) #查看前五行信息

##同方向纠正

Mydata <- harmonise_data(

  exposure_dat=CRP_exp_dat,

  outcome_dat=CHD_out, 

  action= 2

)

#同方向纠正

dim(Mydata)

head(Mydata)

##进行分析啦

res <- mr(Mydata)

View(Mydata)

write.csv(res, file="D:\\MDR database\\res.csv")

generate_odds_ratios(res)#算oR值

#异质性

mr_heterogeneity(Mydata, method_list=c("mr_egger_regression", "mr_ivw")) 

#异 质性检 验--IVworMRegcpleio <- mr_pleiotropy_test(Mydata)

#多效性检验-MR eggerView(pleio)

pleio <- mr_pleiotropy_test(Mydata)

pleio

write.csv(pleio, file="D:\\MDR database\\pleio.csv")

#留一法检验敏感性

single <- mr_leaveoneout(Mydata) 

mr_leaveoneout_plot(single)

write.csv(single, file="D:\\MDR database\\single.csv")

##散点图

mr_scatter_plot(res, Mydata)

###首先使用mr_singlesnp获取单个SNP的结果，之后做森林图

res_single <- mr_singlesnp(Mydata)

mr_forest_plot(res_single)#森林图

mr_funnel_plot(res_single)#漏洞图

#多变量孟德尔

id_exposure <- c("ieu-b-110", "ieu-b-109","ieu-b-111")

id_outcome <- "ebi-a-GCST005195"

exposure_dat <- mv_extract_exposures(id_exposure)

outcome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome)

mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)

res <- mv_multiple(mvdat)

res