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

We estimate the covariance using column-centered Z scores

D.center = apply(as.matrix(data$max$z), 2, function(x) x - mean(x))
mash_data_center = mashr::set_mash_data(Bhat = as.matrix(D.center))

Generate covariance matrices

Besides the 3 types of data-driven covariance matrices created in Immune_Cov_EZ, we have 2 more types of covariance: \[ \tilde{Z}_{n\times 7} = LF' + E \] \[l_{i} \sim \sum_{m=1}^{28} \hat{\pi}_{m} N_{5}(\hat{\mu}_{m}, \hat{U}_{m})\] The estimated covariance structures for LF’ are \(F\hat{U}_{1}F', \cdots, F\hat{U}_{28}F'\).

The results are in Immune_Flash_MVN.

U.loading = readRDS('../output/Immune.flash.load.reduce.cov.rds')

We apply Extreme Deconvolution to refine these 4 types of matrices.

Additionally we include 2 other types of covariance matrices:

• canonical configurations (aka bmalite)

• simple heterogeneity models

FlashResult = readRDS('../output/Immune.flash2.center.greedy.K10.rds')
Flash_res = flash_get_lf(FlashResult$f)

U.flash = c(mashr::cov_from_factors(t(as.matrix(FlashResult$F_flash[,1:5])), "Flash"), 
            list("tFlash" = t(Flash_res) %*% Flash_res / nrow(data$max$z)))

# PCA matrices
U.pca = cov_pca(mash_data_center, P)


# Emperical data matrices
# Denoised data-driven matrices
U.dd = c(U.flash, U.pca, list("XX" = t(D.center) %*% D.center / nrow(data$max$z)), U.loading)

mash_data = mashr::set_mash_data(Bhat = as.matrix(data$max$z))

U.ed = cov_ed(mash_data, U.dd)

# Canonical
U.can = cov_canonical(mash_data) 

# save Covariance matrix
saveRDS(list(Ulist = c(U.ed, U.can), DD_raw = U.dd), 
        paste0('../output/ImmuneEZ.U.center.xtx.Load.K',K,'.P',P,'.rds'))

Session information

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.2

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.4-0            ExtremeDeconvolution_1.3
[3] mashr_0.2-4              ashr_2.1-27             

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.14      knitr_1.17        magrittr_1.5     
 [4] MASS_7.3-47       doParallel_1.0.11 pscl_1.5.2       
 [7] SQUAREM_2017.10-1 lattice_0.20-35   foreach_1.4.4    
[10] plyr_1.8.4        stringr_1.2.0     tools_3.4.3      
[13] parallel_3.4.3    grid_3.4.3        rmeta_2.16       
[16] git2r_0.20.0      htmltools_0.3.6   iterators_1.0.9  
[19] assertthat_0.2.0  yaml_2.1.16       rprojroot_1.2    
[22] digest_0.6.13     Matrix_1.2-12     codetools_0.2-15 
[25] evaluate_0.10.1   rmarkdown_1.8     stringi_1.1.6    
[28] compiler_3.4.3    backports_1.1.2   mvtnorm_1.0-6    
[31] truncnorm_1.0-7