Mr.ash2 parameter estimation

Results for 10 simulations runs, using different initiate and update strategy

show_param <- function(tags, tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  parf <- lapply(f, '[[', "par")
  param <- do.call(rbind, lapply(parf, function(x) t(read.table(x))[2:1,]))
  truth <- param[1:(nrow(param)/2)*2-1,]
  est <- param[1:(nrow(param)/2)*2,]
  outdt <- matrix(0, ncol = 2*ncol(param), nrow = nrow(param)/2)
  outdt[,c(1,3,5,7)] <- truth
  outdt[,c(2,4,6,8)] <- est
  outdt <-cbind(1:nrow(outdt),outdt)
  colnames(outdt) <- c("Simulation#", paste0(rep(c("Truth","Est."),4)))
  knitr::kable(outdt) %>%
  kable_styling("striped") %>%
  add_header_above(c(" " = 1, "Gene.pi1" = 2, "Gene.PVE" = 2, "SNP.pi1" = 2, "SNP.PVE" =2))
}

show_param(tags, tag2s[1])

show_param(tags, tag2s[2])

show_param(tags, tag2s[3])

show_param(tags, tag2s[4])

Regional mr.ash2s PIP overview

Take simulation 1 (NULL; expr-snp; expr-snp) as examples. We use region size 500kb and PIP cut off at 0.5 for SUSIE.

chrom = 18
f <- get_files(tag= "20200721-1-2" , tag2 = tag2s[1])
allchr <- read.table(f[["rpip"]], header = T)
a <- allchr[allchr["chrom"]==chrom,]
print(paste("plot for chr", chrom))

[1] "plot for chr 18"

par(mar=c(5, 4, 4, 6) + 0.1)
with(a, plot(p0, rPIP, col ='salmon', xlab = "position", ylab= "Sum of PIP", type = 'h', lwd = 2))
par(new = T)
with(a, plot(p0, nCausal, pch =19, col = "darkgreen",axes = FALSE, bty = "n", xlab = "", ylab = ""))
axis(side = 4)
mtext(side = 4, line = 3, 'No. causal signals')
legend("topleft",
       legend=c("Mr.ASH PIP", "# Causal"),
       lty=c(1,0), pch=c(NA, 19), col=c("salmon", "darkgreen"))
grid()

PIP calibration

We run 50 simulations and combine results.

#' s is pip or fdr.
cp_plot <- function(s, ifcausal, mode = c("PIP", "FDR"), main = mode[1]){
  # ifcausal:0,1
  a_bin <- cut(s, breaks= seq(0,1, by=0.1))
  if (mode == "PIP") {
     expected = c(by(s, a_bin, FUN = mean))
     observed = c(by(ifcausal, a_bin, FUN = mean))
  } else if (mode == "FDR"){
     expected = c(by(s, a_bin, FUN = mean))
     observed = 1 - c(by(ifcausal, a_bin, FUN = mean))
  }
  plot(expected, observed, xlim= c(0,1), ylim=c(0,1), pch =19, main = main)
  lines(x = c(0,1), y = c(0,1), col ="red")
}

caliPIP_plot <- function(tags, tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags

  .tagname <- function(x, flist){
    a <- read.table(flist[[x]], header =T)
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
 
  res <- merge(mrashres, susieres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  res <- rename(res, c("PIP" = "mr.ash_PIP", "pip" = "SUSIE_PIP", "pip.null" = "SUSIE_PIP_null", "pip.w0" = "SUSIE_PIP_w0"))
  par(mfrow=c(1,4), mar=c(5, 6, 4, 1))
  cp_plot(res$mr.ash_PIP, res$ifcausal, main = "Mr.ash PIP")
  cp_plot(res$SUSIE_PIP, res$ifcausal, main = "SUSIE PIP")
  cp_plot(res$SUSIE_PIP_null, res$ifcausal, main = "SUSIE PIP null")
  cp_plot(res$SUSIE_PIP_w0, res$ifcausal, main = "SUSIE PIP weighted null")
}

caliFDR_plot <- function(tags, tag2){
  
  f <- lapply(tags, get_files, tag2 = tag2)
  gwasf <- lapply(f, '[[', "ggwas")
  names(gwasf) <- tags

  .tagname <- function(x, flist, colnames = NULL){
    a <- read.table(flist[[x]], header =T)
    if (!is.null(colnames)){
      colnames(a) <- colnames
    }
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  
  .addFDR <- function(res){
    res$FDR <- p.adjust(res$PVALUE, method = "fdr")
    res
  }
  
  
  gwasres <- do.call(rbind, lapply(lapply(tags, .tagname, flist = gwasf, 
                                   colnames =  c("chr", "p0",   "p1", "name", 
                                                 "Estimate", "Std.Error", "t-value", "PVALUE")), .addFDR))
  
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags

  .tagname <- function(x, flist){
    a <- read.table(flist[[x]], header =T)
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
 
  res <- merge(mrashres, susieres, by = "name", all = T)
  
  res <- merge(res, gwasres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  
  cp_plot(res$FDR, res$ifcausal, mode ="FDR", main = "TWAS FDR")
  cat("FDR at bonferroni corrected p = 0.05: ", 1 - mean(res[res$PVALUE < 0.05 /dim(res)[1], "ifcausal"]))
}

tag2 = "zeroes-es"
tags_ext <- Reduce(intersect, get_tags(tagglob, tagextr, tag2 = tag2)['gsusie'])
res <- caliPIP_plot(tags = tags_ext, tag2 = tag2)

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

caliFDR_plot(tags = tags_ext, tag2 = tag2)

FDR at bonferroni corrected p = 0.05:  0.645933

caliPIP_plot(tags = tags_ext, tag2 = "lassoes-es")

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

Warning in if (mode == "PIP") {: the condition has length > 1 and only the
first element will be used

caliFDR_plot(tags = tags_ext, tag2 = "lassoes-es")

FDR at bonferroni corrected p = 0.05:  0.647343

PIP scatter plot

mr.ash2s PIP vs. susie PIP.

scatter_plot_PIP<- function(tags, tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags

  .tagname <- function(x, flist){
    a <- read.table(flist[[x]], header =T)
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
 
  res <- merge(mrashres, susieres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  res <- rename(res, c("PIP" = "mr.ash_PIP", "pip" = "SUSIE_PIP", "pip.null" = "SUSIE_PIP_null", "pip.w0" = "SUSIE_PIP_w0") )
  res$ifcausal <- mapvalues(res$ifcausal, 
          from=c(0,1), 
          to=c("Non causal", "Causal"))
  
  fig1 <- plot_ly(data = res, x = ~ mr.ash_PIP, y = ~ SUSIE_PIP, color = ~ ifcausal, 
                 colors = c( "salmon", "darkgreen"))
  
  fig2 <- plot_ly(data = res, x = ~ mr.ash_PIP, y = ~ SUSIE_PIP_null, color = ~ ifcausal, 
                 colors = c( "salmon", "darkgreen"))
  
  fig3 <- plot_ly(data = res, x = ~ mr.ash_PIP, y = ~ SUSIE_PIP_w0, color = ~ ifcausal, 
                 colors = c( "salmon", "darkgreen"))
  
  fig <- subplot(fig1, fig2, fig3, titleX = TRUE, titleY = T, margin = 0.05)
  fig
}

scatter_plot_PIP(tags, tag2s[1])

Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

scatter_plot_PIP(tags, tag2s[2])

scatter_plot_PIP(tags, tag2s[3])

scatter_plot_PIP(tags, tag2s[4])

ROC curve

ROC_plot<- function(tags, tag2){
  f <- lapply(tags, get_files, tag2 = tag2)
  mrashf <- lapply(f, '[[', "gmrash")
  names(mrashf) <- tags
  
  susief <- lapply(f, '[[', "gsusie")
  names(susief) <- tags
  
  gwasf <- lapply(f, '[[', "ggwas")
  names(gwasf) <- tags

  .tagname <- function(x, flist, colnames = NULL){
    a <- read.table(flist[[x]], header =T)
    if (!is.null(colnames)){
      colnames(a) <- colnames
    }
    a[, "name"] <- paste0(x, ":", a[, "name"])
    a
  }
  mrashres <- do.call(rbind, lapply(tags, .tagname, flist = mrashf))
  susieres <- do.call(rbind, lapply(tags, .tagname, flist = susief))
  gwasres <- do.call(rbind, lapply(tags, .tagname, flist = gwasf, 
                                   colnames =  c("chr", "p0",   "p1", "name", "Estimate", "Std.Error", "t-value", "PVALUE")))

  res <- merge(mrashres, susieres, by = "name", all = T)
  res <- merge(res, gwasres, by = "name", all = T)
  
  res <- res[complete.cases(res),]
  res <- rename(res, c("PIP" = "mr.ash", "pip" = "SUSIE", "pip.null"= "SUSIE.null", "pip.w0" = "SUSIE.w0", "PVALUE" = "TWAS") )
  res[,"TWAS"] <- -log10(res[, "TWAS"])
  
  roccolors <-  c("red", "green", "orange", "pink", "blue")
  methods <- c("mr.ash", "SUSIE", "SUSIE.null", "SUSIE.w0", "TWAS")
  plot(0, xlim=c(0,1), ylim=c(0,1), col="white", xlab = "FPR", ylab = "TPR")
  for (i in 1:length(methods)){
    method <- methods[i]
    bordered <- res[order(res[,method]),] 
    actuals <- bordered$ifcausal == 1
    sens <- (sum(actuals) - cumsum(actuals))/sum(actuals)
    spec <- cumsum(!actuals)/sum(!actuals)
    lines(1 - spec, sens, type = "l", col = roccolors[i])
    abline(c(0,0),c(1,1))
    auc <- sum(spec*diff(c(0, 1 - sens)))
    cat("AUC for ", method, ": ", auc)
  }
  legend(0.6,0.3, legend= methods, col=roccolors, lty=1, cex=0.5 )
  grid()
}

tags <- paste0('20200721-1-', c(2,4:9))
ROC_plot(tags, tag2s[2])

AUC for  mr.ash :  0.6894131AUC for  SUSIE :  0.6824739AUC for  SUSIE.null :  0.6762515AUC for  SUSIE.w0 :  0.6835701AUC for  TWAS :  0.7488581

ROC_plot(tags, tag2s[2])

AUC for  mr.ash :  0.6894131AUC for  SUSIE :  0.6824739AUC for  SUSIE.null :  0.6762515AUC for  SUSIE.w0 :  0.6835701AUC for  TWAS :  0.7488581

ROC_plot(tags, tag2s[3])

AUC for  mr.ash :  0.6778905AUC for  SUSIE :  0.6857911AUC for  SUSIE.null :  0.6818999AUC for  SUSIE.w0 :  0.6799617AUC for  TWAS :  0.7489982

ROC_plot(tags, tag2s[4])

AUC for  mr.ash :  0.7589313AUC for  SUSIE :  0.7235847AUC for  SUSIE.null :  0.7043028AUC for  SUSIE.w0 :  0.7258811AUC for  TWAS :  0.768614

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] kableExtra_1.1.0    stringr_1.4.0       plyr_1.8.6         
[4] tidyr_0.8.3         plotly_4.9.2.9000   ggplot2_3.3.1      
[7] data.table_1.12.7   mr.ash.alpha_0.1-34

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.0  purrr_0.3.4       lattice_0.20-38  
 [4] colorspace_1.3-2  vctrs_0.3.1       generics_0.0.2   
 [7] htmltools_0.3.6   viridisLite_0.3.0 yaml_2.2.0       
[10] rlang_0.4.6       later_0.7.5       pillar_1.4.4     
[13] glue_1.4.1        withr_2.1.2       lifecycle_0.2.0  
[16] munsell_0.5.0     gtable_0.2.0      workflowr_1.6.2  
[19] rvest_0.3.2       htmlwidgets_1.3   evaluate_0.12    
[22] knitr_1.20        crosstalk_1.0.0   httpuv_1.4.5     
[25] highr_0.7         Rcpp_1.0.4.6      xtable_1.8-3     
[28] readr_1.3.1       promises_1.0.1    scales_1.0.0     
[31] backports_1.1.2   webshot_0.5.1     jsonlite_1.6.1   
[34] mime_0.6          fs_1.3.1          hms_0.4.2        
[37] digest_0.6.25     stringi_1.3.1     shiny_1.2.0      
[40] dplyr_1.0.0       grid_3.5.1        rprojroot_1.3-2  
[43] tools_3.5.1       magrittr_1.5      lazyeval_0.2.1   
[46] tibble_3.0.1      crayon_1.3.4      pkgconfig_2.0.2  
[49] ellipsis_0.3.1    Matrix_1.2-15     xml2_1.2.0       
[52] rmarkdown_1.10    httr_1.4.1        rstudioapi_0.11  
[55] R6_2.3.0          git2r_0.26.1      compiler_3.5.1