Mash

Mash

mash.data = mash_set_data(Bhat = data$beta, pval = data$p)

Top genes:

# find strong genes
m.1by1 = mash_1by1(mash.data, alpha=0)
strong = get_significant_results(m.1by1)
Z = mash.data$Bhat/mash.data$Shat
Z.strong = Z[strong,]
# center
Z.center = apply(Z.strong, 2, function(x) x - mean(x))

flash.data = flash_set_data(Z.center)
fmodel = flash(flash.data, greedy = TRUE, backfit = TRUE)

fitting factor/loading 1

fitting factor/loading 2

fitting factor/loading 3

fitting factor/loading 4

factors = flash_get_ldf(fmodel)$f
row.names(factors) = colnames(data$beta)
pve.order = order(flash_get_pve(fmodel), decreasing = TRUE)

par(mfrow=c(1,3))
for(i in pve.order){
  barplot(factors[,i], main=paste0('Factor ',i, ' pve= ', round(flash_get_pve(fmodel)[i],3)), las=2, cex.names = 0.7)
}

par(mfrow=c(1,1))

loading = fmodel$EL[,1:3]
colnames(loading) = paste0('Factor',seq(1,3))
flash.loading = flash_set_data(loading)
flmodel = flash(flash.loading, greedy = TRUE, backfit = TRUE)

fitting factor/loading 1

mod = Mclust(loading)
summary(mod$BIC)

Best BIC values:
             VVE,6        VVE,7        EVE,6
BIC      -34778.07 -34815.86606 -34816.67977
BIC diff      0.00    -37.79266    -38.60637

U_list = alply(mod$parameters$variance$sigma,3)
mu_list = alply(mod$parameters$mean,2)
ll = list()
for (i in 1:length(U_list)){
  ll[[i]] = U_list[[i]] + mu_list[[i]] %*% t(mu_list[[i]])
}

Factors = fmodel$EF[,1:3]
U.loading = lapply(ll, function(U){Factors %*% (U %*% t(Factors))})
names(U.loading) = paste0('Load', "_", (1:length(U.loading)))

# rank 1
Flash_res = flash_get_lf(fmodel)
U.Flash = c(mashr::cov_from_factors(t(as.matrix(factors)), "Flash"), 
            list("tFlash" = t(Flash_res) %*% Flash_res / nrow(Z.center)))

U.pca = cov_pca(mash_set_data(Z.center), 3)

U.c = cov_canonical(mash_set_data(Z.center))

U.dd = c(U.pca, U.loading, U.Flash, list('XX' = t(Z.center) %*% Z.center / nrow(Z.center)))

U.ed = cov_ed(mash.data, U.dd, strong)

mash.model = mash(mash.data, c(U.c, U.ed))

 - Computing 11200 x 651 likelihood matrix.
 - Likelihood calculations took 7.98 seconds.
 - Fitting model with 651 mixture components.

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol

Warning in mixIP(matrix_lik = structure(c(0.000549917574123474,
0.107339052191479, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.

 - Model fitting took 19.68 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.69 seconds.

Result

The log-likelihood of fit is

get_loglik(mash.model)

[1] 74461.66

Here is a plot of weights learned:

options(repr.plot.width=12, repr.plot.height=4)
barplot(get_estimated_pi(mash.model), las = 2, cex.names = 0.7)

Check XX covariance matrix:

x           <- mash.model$fitted_g$Ulist[["ED_XX"]]
colnames(x) <- colnames(data$beta)
rownames(x) <- colnames(x)
corrplot(x, method='color', cl.lim=c(-0.2,1), type='upper', addCoef.col = "black", tl.col="black", tl.srt=45, col=colorRampPalette(c("blue","white","red"))(200))

layout(matrix(c(1,2,3,4), 2, 2, byrow=TRUE))
svd.out = svd(mash.model$fitted_g$Ulist[["ED_XX"]])
v = svd.out$v
colnames(v) = colnames(data$beta)
rownames(v) = colnames(v)
options(repr.plot.width=10, repr.plot.height=5)
for (j in 1:4)
  barplot(v[,j]/v[,j][which.max(abs(v[,j]))], cex.names = 0.7,
          las = 2, main = paste0("EigenVector ", j, " for XX"))

There are 4504 diferentially expressed genes.

Check pairwise sharing by sign:

x = get_pairwise_sharing(mash.model, factor=0)
x[x > 1]    <- 1
x[x < -1]   <- -1
colnames(x) <- colnames(data$beta)
rownames(x) <- colnames(x)
corrplot.mixed(x, tl.pos="d",upper='color', cl.lim=c(0,1), upper.col=colorRampPalette(rev(c("#D73027","#FC8D59","#FEE090","#FFFFBF",
                               "#E0F3F8","#91BFDB","#4575B4")))(40),
               tl.cex=1.2)

Check pairwise sharing by magnitude and sign:

x = get_pairwise_sharing(mash.model)
x[x > 1]    <- 1
x[x < -1]   <- -1
colnames(x) <- colnames(data$beta)
rownames(x) <- colnames(x)
corrplot.mixed(x, tl.pos="d",upper='color', cl.lim=c(0,1), upper.col=colorRampPalette(rev(c("#D73027","#FC8D59","#FEE090","#FFFFBF",
                               "#E0F3F8","#91BFDB","#4575B4")))(40),
               tl.cex=1.2)

• Effects compare with SCZ

CompareSCZ = diag(5)
CompareSCZ[,2] = -1
CompareSCZ = CompareSCZ[-2,]
row.names(CompareSCZ) = colnames(data$beta)[-2]

mash.model.SCZ = mash.model
mash.model.SCZ$result = mash_compute_posterior_matrices(mash.model, mash.data, A=CompareSCZ, algorithm.version = 'R')

Check pairwise sharing by sign:

x = get_pairwise_sharing(mash.model.SCZ, factor=0)
x[x > 1]    <- 1
x[x < -1]   <- -1
colnames(x) <- row.names(CompareSCZ)
rownames(x) <- colnames(x)
corrplot.mixed(x, tl.pos="d",upper='color', cl.lim=c(0,1), upper.col=colorRampPalette(rev(c("#D73027","#FC8D59","#FEE090","#FFFFBF",
                               "#E0F3F8","#91BFDB","#4575B4")))(40),
               tl.cex=1.2)

Check pairwise sharing by magnitude and sign:

x = get_pairwise_sharing(mash.model.SCZ)
x[x > 1]    <- 1
x[x < -1]   <- -1
colnames(x) <- row.names(CompareSCZ)
rownames(x) <- colnames(x)
corrplot.mixed(x, tl.pos="d",upper='color', cl.lim=c(0,1), upper.col=colorRampPalette(rev(c("#D73027","#FC8D59","#FEE090","#FFFFBF",
                               "#E0F3F8","#91BFDB","#4575B4")))(40),
               tl.cex=1.2)

Session information

sessionInfo()

R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

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] ggplot2_2.2.1  corrplot_0.84  colorRamps_2.3 flashr_0.5-6  
[5] plyr_1.8.4     mclust_5.4     mashr_0.2-6    ashr_2.2-7    
[9] limma_3.34.9  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.16             pillar_1.2.2            
 [3] compiler_3.4.4           git2r_0.21.0            
 [5] workflowr_1.0.1          R.methodsS3_1.7.1       
 [7] R.utils_2.6.0            iterators_1.0.9         
 [9] tools_3.4.4              digest_0.6.15           
[11] tibble_1.4.2             gtable_0.2.0            
[13] evaluate_0.10.1          lattice_0.20-35         
[15] rlang_0.2.0              Matrix_1.2-14           
[17] foreach_1.4.4            yaml_2.1.19             
[19] parallel_3.4.4           mvtnorm_1.0-7           
[21] ebnm_0.1-11              stringr_1.3.0           
[23] knitr_1.20               REBayes_1.3             
[25] rprojroot_1.3-2          grid_3.4.4              
[27] rmarkdown_1.9            rmeta_3.0               
[29] magrittr_1.5             whisker_0.3-2           
[31] scales_0.5.0             backports_1.1.2         
[33] codetools_0.2-15         htmltools_0.3.6         
[35] MASS_7.3-50              assertthat_0.2.0        
[37] softImpute_1.4           colorspace_1.3-2        
[39] stringi_1.2.2            Rmosek_8.0.69           
[41] lazyeval_0.2.1           munsell_0.4.3           
[43] doParallel_1.0.11        pscl_1.5.2              
[45] truncnorm_1.0-8          SQUAREM_2017.10-1       
[47] ExtremeDeconvolution_1.3 R.oo_1.22.0