Set up the data

# Load required packages
library(mashr); library(ExtremeDeconvolution); library(flashr2)

Loading required package: ashr

# read data
data = readRDS('../data/ImmuneQTLSummary.4MASH.rds')
data$max$se = data$max$beta/data$max$z
data$null$se = data$null$beta / data$null$z

# set parameters
K = 10
P = 5
vhat = 1

Fit mash model

Using a random null set, the cell below computes the weights for input covariance matrices (Covariance matrices) in MASH mixture. The output contains matrix of log-likelihoods as well as weights learned from the hierarchical model.

if (vhat == 1) {
  V = cor(data$null$z[which(apply(abs(data$null$z),1, max) < 2),])
} else {
  V = diag(ncol(data$null$z))
}

mash_data = mashr::set_mash_data(Bhat = as.matrix(data$null$beta), 
                                 Shat = as.matrix(data$null$se), 
                                 V=as.matrix(V), alpha=1)

Ulist = readRDS('../output/ImmuneEZ.U.center.xtx.K10.P5.rds')$Ulist

saveRDS(mashr::mash(mash_data, Ulist, outputlevel = 1), 
        paste0('../output/ImmuneEZ.V',vhat,'.center.mash_model.K',K,'.P',P,'.rds'))

 - Computing 193301 x 476 likelihood matrix.
 - Likelihood calculations took 11.69 seconds.
 - Fitting model with 476 mixture components.
 - Model fitting took 325.46 seconds.

Posterior inference

Applying hyperparameters learned from the training (the null) set to the test (the top eQTL) set, the cell below computes posterior quantities.

U.m = readRDS('../output/ImmuneEZ.V1.center.mash_model.K10.P5.rds')

if (vhat == 1) {
  V = cor(data$null$z[which(apply(abs(data$null$z),1, max) < 2),])
} else {
  V = diag(ncol(data$null$z))
}

mash_data = mashr::set_mash_data(Bhat = as.matrix(data$max$beta), 
                                 Shat = as.matrix(data$max$se), 
                                 V = as.matrix(V), 
                                 alpha = 1)

saveRDS(mash_compute_posterior_matrices(U.m, mash_data), 
        paste0('../output/ImmuneEZ.V',vhat,'.center.mash_posterior.K',K,'.P',P,'.rds'))

FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE

Session information

sessionInfo()

R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.1

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] flashr2_0.2-3            ExtremeDeconvolution_1.3
[3] mashr_0.2-4              ashr_2.1-27             

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.13      compiler_3.4.2    git2r_0.19.0     
 [4] plyr_1.8.4        iterators_1.0.8   tools_3.4.2      
 [7] digest_0.6.12     evaluate_0.10.1   tibble_1.3.4     
[10] gtable_0.2.0      lattice_0.20-35   rlang_0.1.2      
[13] Matrix_1.2-11     foreach_1.4.3     yaml_2.1.14      
[16] parallel_3.4.2    mvtnorm_1.0-6     stringr_1.2.0    
[19] knitr_1.17        REBayes_0.85      rprojroot_1.2    
[22] grid_3.4.2        rmarkdown_1.7     rmeta_2.16       
[25] ggplot2_2.2.1     magrittr_1.5      backports_1.1.1  
[28] scales_0.5.0      codetools_0.2-15  htmltools_0.3.6  
[31] MASS_7.3-47       assertthat_0.2.0  colorspace_1.3-2 
[34] stringi_1.1.5     Rmosek_8.0.69     lazyeval_0.2.1   
[37] pscl_1.5.2        doParallel_1.0.11 munsell_0.4.3    
[40] truncnorm_1.0-7   SQUAREM_2017.10-1