直接复现一篇5分多的中介孟德尔随机化文章

今天给大家分享如何复现一篇5分多的中介孟德尔随机化文章，这篇文章发表在Frontiers in Endocrinology期刊上，影响因子：5.2，中科院2区，文章题目如下：

作者采用two step的方法来计算中介效应，具体分析思路和步骤如下：

第一步：RA to AAS TWMR(得到总效应，beta_all)

第二步：AAS to RATWMR(确定可以做中介)

第三步：RA to CRP TWMR(得到beta1)

第四步：CRP to AAS TWMR(得到beta2)

中介效应：beta12=beta1*beta2

直接效应：beta_dir=beta_all-beta12

具体分析代码如下：

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


#install.packages("ggplot2")


#library
library(TwoSampleMR)
library(ggplot2)


#set up your work directory
setwd("")




#######two step###########
#step 1 RA to AAS
expo_rt<-extract_instruments(outcome="ukb-b-11874",clump = FALSE)
expo_rt<-clump_data(expo_rt,clump_kb = 10000,
                    clump_r2 = 0.001,
                    clump_p1 = 5e-08,
                    clump_p2 = 5e-08,
)


outc_rt <- extract_outcome_data(expo_rt$SNP, outcomes = "finn-b-M13_ATLOAXSUBLUX")


#####################################


harm_rt <- harmonise_data(
  exposure_dat =  expo_rt, 
  outcome_dat = outc_rt,action=2)


mr_result<- mr(harm_rt)
View(mr_result)
result_or=generate_odds_ratios(mr_result)
write.table(result_or[,5:ncol(result_or)],"OR1.txt",row.names = F,sep = "	",quote = F)
#####################################


# step 2 AAS to RA
expo_rt2<-extract_instruments(outcome="finn-b-M13_ATLOAXSUBLUX",clump = FALSE)
expo_rt2<-clump_data(expo_rt,clump_kb = 10000,
                    clump_r2 = 0.001,
                    clump_p1 = 5e-08,
                    clump_p2 = 5e-08,
)


outc_rt2 <- extract_outcome_data(expo_rt2$SNP, outcomes = "ukb-b-11874")


harm_rt2 <- harmonise_data(
  exposure_dat =  expo_rt2, 
  outcome_dat = outc_rt2,action=2)


mr_result2<- mr(harm_rt2)
View(mr_result2)
result_or2=generate_odds_ratios(mr_result2)
write.table(result_or2[,5:ncol(result_or2)],"OR2.txt",row.names = F,sep = "	",quote = F)
#####################################


#step 3 RA to CRP
expo_rt3<-expo_rt
outc_rt3 <- extract_outcome_data(expo_rt3$SNP, outcomes = "ieu-b-4764")


harm_rt3 <- harmonise_data(
  exposure_dat =  expo_rt3, 
  outcome_dat = outc_rt3,action=2)
write.table(harm_rt3, "harmonise3.txt",row.names = F,sep = "	",quote = F)


mr_result3<- mr(harm_rt3)
View(mr_result3)
result_or3=generate_odds_ratios(mr_result3)
write.table(result_or3[,5:ncol(result_or3)],"OR3.txt",row.names = F,sep = "	",quote = F)
#####################################


#step 4 CRP to AAS
expo_rt4<-extract_instruments(outcome="ieu-b-4764",clump = FALSE)
expo_rt4<-clump_data(expo_rt4,clump_kb = 10000,
                     clump_r2 = 0.001,
                     clump_p1 = 5e-08,
                     clump_p2 = 5e-08,
)


outc_rt4 <- extract_outcome_data(expo_rt4$SNP, outcomes = "finn-b-M13_ATLOAXSUBLUX")
#####################################
####关注微信公众号生信狂人团队
###遇到代码报错等不懂的问题可以添加微信scikuangren进行答疑
###作者邮箱：sxkrteam@shengxinkuangren.com
harm_rt4 <- harmonise_data(
  exposure_dat =  expo_rt4, 
  outcome_dat = outc_rt4,action=2)
write.table(harm_rt4, "harmonise4.txt",row.names = F,sep = "	",quote = F)


mr_result4<- mr(harm_rt4)
View(mr_result4)
result_or4=generate_odds_ratios(mr_result4)
write.table(result_or4[,5:ncol(result_or4)],"OR4.txt",row.names = F,sep = "	",quote = F)
#####################################






#中介效应
beta12=beta1*beta2
#直接效应
beta_dir=beta_all-beta12

通过上面的代码就可以计算出总效应、中介效应和直接效应，当然还有一些细节问题没有办法一一放在代码中展示，需要大家自行探索，希望上面的参考文章和代码对大家有所帮助。