setwd("D:\\MDR database")

library(TwoSampleMR)

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

#參數(shù)默認(rèn)為P=5*10-8，r方=0.001，kb=10000

#查看暴露SNP基本情況

dim(CRP_exp_dat)? #查看篩選出多少個(gè)SNP

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

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

##結(jié)局指標(biāo)：在暴露因素的SNP里面找結(jié)局指標(biāo)的（從56--52個(gè)）

CHD_out <- extract_outcome_data(

??

? ? snps=CRP_exp_dat$SNP, #根據(jù)暴露SNP編號(hào)在結(jié)局中查找對用SNP

? ? outcomes="ieu-a-7", #查找數(shù)據(jù)庫ID號(hào)為ebi-a-GCST006906的SNP

? ? proxies = TRUE,#在結(jié)局中找不到的SNP是否用1000組的代替，默認(rèn)TRUE，R方0.8

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

? ? access_token = NULL#是否通過谷歌訪問，國內(nèi)選NULL

? )

#查看暴露SNP基本情況，發(fā)現(xiàn)少了幾個(gè)，

#一般可能不符合條件，或者再結(jié)局中沒有找到

dim(CHD_out) #查看篩選出多少個(gè)SNP

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

##同方向糾正

Mydata <- harmonise_data(

? exposure_dat=CRP_exp_dat,

? outcome_dat=CHD_out,?

? action= 2

)

#同方向糾正

dim(Mydata)

head(Mydata)

##進(jìn)行分析啦

res <- mr(Mydata)

View(Mydata)

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

generate_odds_ratios(res)#算oR值

#異質(zhì)性

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

#異 質(zhì)性檢 驗(yàn)--IVworMRegcpleio <- mr_pleiotropy_test(Mydata)

#多效性檢驗(yàn)-MR eggerView(pleio)

pleio <- mr_pleiotropy_test(Mydata)

pleio

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

#留一法檢驗(yàn)敏感性

single <- mr_leaveoneout(Mydata)?

mr_leaveoneout_plot(single)

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

##散點(diǎn)圖

mr_scatter_plot(res, Mydata)

###首先使用mr_singlesnp獲取單個(gè)SNP的結(jié)果，之后做森林圖

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