Mr.ash2 parameter estimation

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

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-3" , 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()

Past versions of “rpip 1-1.png”

PIP calibration

We run 50 simulations and combine results.

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

a <- caliFDR_plot(tags = tags_ext, tag2 = tag2)

FDR at bonferroni corrected p = 0.05:  0.645933

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

a <- 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(tags, tag2s[1])

scatter_plot_PIP(tags, tag2s[2])

scatter_plot_PIP(tags, tag2s[3])

scatter_plot_PIP(tags, tag2s[4])

ROC curve

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

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

ROC_plot(tags, tag2s[2])

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

ROC_plot(tags, tag2s[3])

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

ROC_plot(tags, tag2s[4])

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

PIP vs p value

a <- scatter_plot_PIP_p(tags, "lassoes-es")

Pretty plots

1. PIP calibration plot; B) FDP inflation of TWAS; C) Scatter plot of individual genes: PIPs vs. -log10-p values from TWAS.

We use SNP genotype data from chr 17 to chr 22 combined. These genomic regions represents 12.5% of the genome. SNPs are downsampled to 1/10 (randomly), eQTLs used in building the expression model were added back. We used lasso to train expression model using GTEx adipose tissue v7. In this simulation, we simulate phenotype with PVE explained by gene as 0.01, by SNP 0.05. pi1 for gene 0.05, pi1 for SNP 0.0025.

For our method, we run mr.ash2s for both genes and SNPs (initiated with lasso). Then apply SUSIE for regions with sum of PIP > 0.5. This gives us PIP for individual genes.

For TWAS, we use linear regression for each gene and obtain p values.

par(mfrow=c(1,3))
cp_plot(res1$SUSIE.w0_PIP, res1$ifcausal, main = "SUSIE PIP calibration", col = "firebrick",cex.lab=1.4, cex.axis=1.4, cex.main=1.4, cex.sub=1.4)
grid()
cp_plot(res2$FDR, res2$ifcausal, mode ="FDR", main = "TWAS FDP", col= "darkgreen",cex.lab=1.4, cex.axis=1.4, cex.main=1.4, cex.sub=1.4)
grid()
plot(res3$TWAS, res3$SUSIE.w0, col = ifelse(res3$ifcausal=="Causal",  "salmon", "cyan3"), xlab = expression('TWAS -log'[10]*'(p)'), ylab = "SUSIE PIP",pch =19,cex.lab=1.4, cex.axis=1.4, cex.main=1.4, cex.sub=1.4)
grid()
legend(8.5, 0.6, 
       legend = c("Causal", "Non causal"), 
      cex=1.1, pch=19, col=c("salmon", "cyan3"))

Past versions of “pretty plot-1.png”

Session information

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.2.1    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] promises_1.0.1    scales_1.0.0      backports_1.1.2  
[31] webshot_0.5.1     jsonlite_1.6.1    mime_0.6         
[34] fs_1.3.1          digest_0.6.25     stringi_1.3.1    
[37] shiny_1.2.0       dplyr_1.0.0       grid_3.5.1       
[40] rprojroot_1.3-2   tools_3.5.1       magrittr_1.5     
[43] lazyeval_0.2.1    tibble_3.0.1      crayon_1.3.4     
[46] whisker_0.3-2     pkgconfig_2.0.2   ellipsis_0.3.1   
[49] Matrix_1.2-15     xml2_1.2.0        rmarkdown_1.10   
[52] httr_1.4.1        rstudioapi_0.11   R6_2.3.0         
[55] git2r_0.26.1      compiler_3.5.1