MR_BMA分析记录

setwd('H:\\Puberty_SEM\\MR_GWAS\\demo_AMD-master')
library(TwoSampleMR)
source("./demo_AMD-master/summary_mvMR_SSS.R")
source("./demo_AMD-master/summary_mvMR_BF.R")
library(ggplot2)
library(TwoSampleMR)
library(tibble)
library(dplyr)
library(Hmisc)
exps= mv_extract_exposures_local(
  filenames_exposure = c('Alanine_aminotransferase.vcf','ApoA.vcf','HbA1c.vcf','HDL.vcf'),
  sep = "\t",
  phenotype_col = "Phenotype",
  snp_col = "SNP",
  beta_col = "BETA",
  se_col = "SE",
  eaf_col = "Freq",
  effect_allele_col = "A1",
  other_allele_col = "A2",
  pval_col = "P",
  samplesize_col = "N",
  min_pval = 1e-200,
  log_pval = FALSE,
  pval_threshold = 5e-08,
  plink_bin = 'C:/Users/DELL/Documents/R/win-library/4.1/plinkbinr/bin/plink_Windows.exe',
  bfile = 'H:\\Puberty_SEM\\MR_GWAS\\EUR\\EUR',
  clump_r2 = 0.001,
  clump_kb = 10000,
  pop = "EUR",
  harmonise_strictness = 2
)
write.table(exps,file = 'exps.txt',row.names = F,sep = '\t',quote = F)
exps = read.table('exps.txt',header = T)
outcome_id = c("Coronary_artery_disease.txt", "Heart_failure.txt", "mvAge.txt")
for (i in 1:length(outcome_id)){
    outcome_dat <- read_outcome_data(
    snps = exps$SNP,
    filename = outcome_id[i],
    sep = "\t",
    snp_col = "SNP",
    beta_col = "BETA",
    se_col = "SE",
    effect_allele_col = "A1",
    other_allele_col = "A2",
    pval_col = "P",
    samplesize_col = "N",
    phenotype_col = 'Phenotype')
    mvdat = mv_harmonise_data(exposure_dat=exps, outcome_dat, harmonise_strictness = 2)
    betaX=as.data.frame(mvdat$exposure_beta)
    colnames(betaX)=capitalize(sapply(strsplit(mvdat$expname$exposure,fixed = TRUE, split= " ||"),"[",1))
    obesity_beta = mvdat$outcome_beta
    obesity_se = mvdat$outcome_se
    rs = row.names(betaX)
    genes = row.names(betaX)
    dim(betaX)
    colnames(betaX)
    rf = colnames(betaX)
    betaX_ivw = betaX / obesity_se
    obesity_beta_ivw = obesity_beta / obesity_se
    dim(betaX_ivw)
    length(obesity_beta_ivw)
    set.seed(12345)
    obesity_wbc_input=new("mvMRInput", betaX = as.matrix(betaX_ivw), betaY = as.matrix(obesity_beta_ivw), snps=rs, exposure=rf, outcome = "obesity")
    BMA_output=summarymvMR_SSS(obesity_wbc_input,kmin=1,kmax=length(rf), prior_prob=0.1, max_iter=1000)
    filename1 = paste('amd_best_',outcome_id[i])
    best.model.out = sss.report.best.model(BMA_output, top = 4, write.out = TRUE, csv.file.name=filename1)
    best.model.out
    filename2 = paste('amd_mr_',outcome_id[i])
    mr.bma.out = sss.report.mr.bma(BMA_output, top = 4, write.out = TRUE, csv.file.name=filename2)    
    mr.bma.out  
    filename3 = paste('permutation_',outcome_id[i])
    permute_bma = create.permutations(BMA_output, nrepeat = 10000, save.matrix=TRUE, file.name = filename3)
    empirical.p = calculate.p(BMA_output, permute_bma)
    empirical.p = as.data.frame(empirical.p)
    filename4 = paste('empirical.p_',outcome_id[i])
    write.table(empirical.p,file = filename4)
}