穿插基本概念：
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

方法都看過了较店，簡單總結(jié)：
MVMR很不方便闷堡，水平多效性、異質(zhì)性的檢驗有點沒意義，我沒看到什么方法去排除這些問題，時間有限，不想研究了丰捷。還是TSMR比較高效。所以我選擇兩部TSMR

遺憾：

1. 兩步TSMR來做中介孟德爾很方便寂汇，但是要考慮到水平多效性和異質(zhì)性的話病往，可能要選擇不同的模型來做，這樣就要求SNP數(shù)量要多骄瓣，所以我P設(shè)置在5e-6
2. 中介要求有三個數(shù)據(jù)源停巷，但現(xiàn)在除了一些疾病、性狀有單獨的數(shù)據(jù)榕栏，基本都是UKB/FINNGEN的畔勤，所以大家都搞腸道菌群、分子的QTL就很能理解了
3. 多種計算中介置信區(qū)間的方法扒磁，我只用了一種庆揪，bootstrap方法在這里其實是有局限的，畢竟沒有樣本級別的數(shù)據(jù)

代碼如下妨托，隨緣自取

#exp/med/out 是ieu_id缸榛，不過ieu很難滿足中介，大家自己讀取數(shù)據(jù)吧
#對于F值的計算兰伤，我這里是兩步的算法内颗，對于ncase比較小的性狀還是推薦(beta/se)^2
#路徑自己改哈，我是習(xí)慣每一步文件分開保存
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')