穿插基本概念：
https://zhuanlan.zhihu.com/p/383529186
https://zhuanlan.zhihu.com/p/643496483
中介
https://www.zhihu.com/question/520320330
https://zhuanlan.zhihu.com/p/657940139
https://zhuanlan.zhihu.com/p/348349222
https://zhuanlan.zhihu.com/p/604939175
https://zhuanlan.zhihu.com/p/671120939
https://zhuanlan.zhihu.com/p/660363579
方法都看过了，简单总结：
MVMR很不方便，水平多效性、异质性的检验有点没意义，我没看到什么方法去排除这些问题，时间有限，不想研究了。还是TSMR比较高效。所以我选择两部TSMR
遗憾：
- 两步TSMR来做中介孟德尔很方便，但是要考虑到水平多效性和异质性的话，可能要选择不同的模型来做，这样就要求SNP数量要多，所以我P设置在5e-6
- 中介要求有三个数据源，但现在除了一些疾病、性状有单独的数据，基本都是UKB/FINNGEN的，所以大家都搞肠道菌群、分子的QTL就很能理解了
- 多种计算中介置信区间的方法，我只用了一种，bootstrap方法在这里其实是有局限的，毕竟没有样本级别的数据
代码如下，随缘自取
#exp/med/out 是ieu_id，不过ieu很难满足中介，大家自己读取数据吧
#对于F值的计算，我这里是两步的算法，对于ncase比较小的性状还是推荐(beta/se)^2
#路径自己改哈，我是习惯每一步文件分开保存
for(i in 1:length(exp)){
exp_dt<-exp[i]
exp_exp<-extract_instruments(outcomes = exp_dt, p1=5e-6,access_token=NULL)
exp_exp<-transform(exp_exp,R2=2*((beta.exposure)^2)*eaf.exposure*(1-eaf.exposure))
exp_exp<-transform(exp_exp,F_index=(samplesize.exposure[1]-2)*R2/(1-R2))
exp_exp<-filter(exp_exp,F_index>10)
Sys.sleep(3)
for(j in 1:length(med)){
med_dt<-med[j]
A_dat<-out_harmon(exp_exp,med_dt)
saveRDS(A_dat,paste0('D:/MR2/dats_med/A_dats/',exp_dt,'_',med_dt))
Sys.sleep(5)
}
}
for(i in 1:length(med)){
med_dt<-med[i]
med_exp<-extract_instruments(outcomes = med_dt, p1=5e-6,access_token=NULL)
med_exp<-transform(med_exp,R2=2*((beta.exposure)^2)*eaf.exposure*(1-eaf.exposure))
med_exp<-transform(med_exp,F_index=(samplesize.exposure[1]-2)*R2/(1-R2))
med_exp<-filter(med_exp,F_index>10)
Sys.sleep(3)
for(j in 14:length(out)){
out_dt<-out[j]
B_dat<-out_harmon(med_exp,out_dt)
saveRDS(B_dat,paste0('D:/MR2/dats_med/B_dats/',med_dt,'_',out_dt))
Sys.sleep(5)
}
}
for(i in 1:length(exp)){
exp_dt<-exp[i]
exp_exp<-extract_instruments(outcomes = exp_dt, p1=5e-6,access_token=NULL)
exp_exp<-transform(exp_exp,R2=2*((beta.exposure)^2)*eaf.exposure*(1-eaf.exposure))
exp_exp<-transform(exp_exp,F_index=(samplesize.exposure[1]-2)*R2/(1-R2))
exp_exp<-filter(exp_exp,F_index>10)
Sys.sleep(3)
for(j in 1:length(out)){
out_dt<-out[j]
C_dat<-out_harmon(exp_exp,out_dt)
saveRDS(C_dat,paste0('D:/MR2/dats_med/C_dats/',exp_dt,'_',out_dt))
Sys.sleep(5)
}
}
for(abc in c('A','B','C')){
if(abc=='A'){
dats<-list.files('D:/MR2/dats_med/A_dats/')
save_paths<-paste0('D:/MR2/res_med/A_res/',dats,'.res')
dats<-paste0('D:/MR2/dats_med/A_dats/',dats)
}else if(abc=='B'){
dats<-list.files('D:/MR2/dats_med/B_dats/')
save_paths<-paste0('D:/MR2/res_med/B_res/',dats,'.res')
dats<-paste0('D:/MR2/dats_med/B_dats/',dats)
}else if(abc=='C'){
dats<-list.files('D:/MR2/dats_med/C_dats/')
save_paths<-paste0('D:/MR2/res_med/C_res/',dats,'.res')
dats<-paste0('D:/MR2/dats_med/C_dats/',dats)
}
for (i in 1:length(dats)){
print(i)
print(i)
dat<-read_rds(dats[i])
res<-try(mr_p_re<-mr_presso(BetaOutcome ="beta.outcome", BetaExposure = "beta.exposure",
SdOutcome ="se.outcome",SdExposure = "se.exposure",
OUTLIERtest = TRUE,DISTORTIONtest = TRUE, data = dat,
NbDistribution = 1500,SignifThreshold = 0.05),silent = T)
if(!('try-error' %in% class(res))) {
if(mr_p_re$`MR-PRESSO results`$`Global Test`$Pvalue<0.05){
if(!is.character(mr_p_re$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`)){
dat<-dat[-mr_p_re$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`,]
}
}
}
res<-try(res_h<-mr_heterogeneity(dat, method_list = c("mr_ivw")),silent = T)
if(!('try-error' %in% class(res)) & dim(res_h)[1]>0) {
res_m<-mr_pleiotropy_test(dat)
if(res_m$pval>0.05){
if(res_h$Q_pval<0.05){
res<-mr(dat,method_list = c('mr_ivw_mre'))
}else{
res<-mr(dat,method_list = c('mr_ivw_fe'))
}
}else{
res<-mr(dat,method_list = c('mr_egger_regression'))
}
}else{
res<-mr(dat,method_list = c('mr_ivw_fe'))
}
saveRDS(res,save_paths[i])
}
}
A_res<-list()
A_res_list<-paste0('D:/MR2/res_med/A_res/',list.files('D:/MR2/res_med/A_res/'))
for(i in A_res_list){
A_res<-append(A_res,list(read_rds(i)))
}
re1<-do.call(rbind.data.frame,A_res)
A_or<-generate_odds_ratios(re1)
B_res<-list()
B_res_list<-paste0('D:/MR2/res_med/B_res/',list.files('D:/MR2/res_med/B_res/'))
for(i in B_res_list){
B_res<-append(B_res,list(read_rds(i)))
}
re1<-do.call(rbind.data.frame,B_res)
B_or<-generate_odds_ratios(re1)
C_res<-list()
C_res_list<-paste0('D:/MR2/res_med/C_res/',list.files('D:/MR2/res_med/C_res/'))
for(i in C_res_list){
C_res<-append(C_res,list(read_rds(i)))
}
re1<-do.call(rbind.data.frame,C_res)
C_or<-generate_odds_ratios(re1)
res_med_all<-data.frame(matrix(ncol=8))
colnames(res_med_all)<-c('Exposures','Mediations','Outcomes',
'Total effect','Direct effect A','Direct effect B',
'Mediation effect','P')
for(i in 1:dim(A_or)[1]){
exp_temp<-A_or$id.exposure[i]
for(j in 1:19){
med_temp<-A_or$id.outcome[i]
out_temp<-out[j]
index_A<-A_or$id.exposure==exp_temp & A_or$id.outcome==med_temp
index_B<-B_or$id.exposure==med_temp & B_or$id.outcome==out_temp
index_C<-C_or$id.exposure==exp_temp & C_or$id.outcome==out_temp
if(A_or$pval[index_A]<0.05 &
B_or$pval[index_B]<0.05 &
C_or$pval[index_C]<0.05){
med_res <- medci(mu.x=A_or$b[index_A],
mu.y=B_or$b[index_B],
se.x=A_or$se[index_A],
se.y=B_or$se[index_B],
rho=0, alpha=.05,type='asymp')
Z <- (med_res$Estimate - 0) / med_res$SE
p_value <- 2 * pnorm(-abs(Z))
temp<-data.frame('Exposures'=exp_temp,'Mediations'=med_temp,'Outcomes'=out_temp,
'Total effect'=paste0(round(C_or$b[index_C],3),'(',round(C_or$lo_ci[index_C],3),' to ',round(C_or$up_ci[index_C],3),')'),
'Direct effect A'=paste0(round(A_or$b[index_A],3),'(',round(A_or$lo_ci[index_A],3),' to ',round(A_or$up_ci[index_A],3),')'),
'Direct effect B'=paste0(round(B_or$b[index_B],3),'(',round(B_or$lo_ci[index_B],3),' to ',round(B_or$up_ci[index_B],3),')'),
'Mediation effect'=paste0(round(med_res$Estimate,3),'(',round(med_res$`95% CI`[1],3),' to ',round(med_res$`95% CI`[2],3),')'),
'P'=p_value)
colnames(temp)<-str_replace_all(colnames(temp),'\\.',' ')
res_med_all<-rbind(res_med_all,temp)
}
}
}
res_med_all<-res_med_all[-1,]
res_med_all$Exposures<-work_list$...1[match(res_med_all$Exposures,work_list$...2)]
res_med_all$Mediations<-work_list$...1[match(res_med_all$Mediations,work_list$...2)]
res_med_all$Outcomes<-work_list$...1[match(res_med_all$Outcomes,work_list$...2)]
tab<-as_flextable(xtable(res_med_all))
doc<-read_docx()
doc = body_add_flextable(doc,tab)
print(doc,"D:/tab.docx")
write.csv(C_or[C_or$pval<0.05,],'C_or_p.csv')
