MASH v FLASH workflows for GTEx

Introduction

While I previously compared MASH and FLASH fits on strong tests only and on a random subset of tests, here I propose several workflows that are analogous to the one suggested for the GTEx data in this MASH vignette. For the code used in this analysis, see below.

Fitting Methods

For MASH, I follow the workflow in the vignette linked above, except that I assume that the null tests are uncorrelated (that is, I set \(V = I\)). This is almost certainly not the case, but some more work needs to be done before we can handle the case \(V \ne I\) in FLASH.

The workflows for FLASH proceed along similar lines to the workflow for MASH:

1. I obtain “data-driven” loadings (analogous to MASH’s data-driven covariance matrices) by fitting a FLASH object to the “strong” tests using either the “OHF” method or the “Zero” method described in my simulation study.

2. I fix the loadings obtained in step 1 at their expectation \(EL\) and add them as fixed data-driven loadings to a FLASH object. Further, I add 45 fixed “canonical” loadings (a vector of all ones and a one-hot vector for each of 44 conditions), which can be viewed as analogous to MASH’s canonical covariance matrices. (See the discussion in my vignette.) Since the “Zero” method generates many more data-driven loadings than the “OHF” method, I experiment with only keeping a subset of the data-driven loadings from the “Zero” method.

3. With the loadings fixed, I backfit a FLASH object to the random subset of tests to obtain priors \(g_f\) on the factors. I use the random subset rather than the “strong” subset because I want the priors to hold generally, and not just for the strong tests.

4. Finally, using the same fixed loadings as in step 2 and fixing the priors \(g_f\) at the values obtained in step 3, I backfit a FLASH object to the strong tests to get posterior means and variances.

I fit five FLASH objects, using five variations of the above workflow:

methods <- data.frame(fit = as.character(1:5),
                      data.driven = c("OHF", "Zero (top 5)", 
                                      "Zero (top 10)", "Zero (top 20)",
                                      "Zero (full)"),
                      include.canonical = c(rep("Yes", 4), "No"))
knitr::kable(methods)

Comments on MASH fit

library(mashr)

Loading required package: ashr

fpath <- "./output/MASHvFLASHgtex2/"
m_final <- readRDS(paste0(fpath, "m.rds"))

It took 2.4 minutes to run Extreme Deconvolution to find data-driven covariance matrices. The MASH fit on the random subset of tests (to determine mixture weights) required 2.7 minutes, and the final MASH fit on the strong tests required only 16 seconds. The total fitting time was 5.3 minutes.

The estimated mixture weights were as follows. Note in particular that there are large weights on the data-driven matrices “ED_tPCA” and “ED_PCA_2,” as well as on the canonical “equal_effects” matrix. There are moderate weights on the null matrix, some unique effects (including Testis, Thyroid, and Transformed fibroblasts), and two of the canonical “simple_het” matrices (where effect sizes are assumed to be of equal variance and equally correlated, with correlation coefficients of, respectively, 0.25 and 0.5).

barplot(get_estimated_pi(m_final), las = 2, cex.names = 0.4)

Correlation plots for the data-driven matrices are as follows. The first (“ED_tPCA”) describes effects that are shared across a handful of tissues (which, notably, does not include brain tissues).

library(corrplot)

corrplot 0.84 loaded

corrplot(m_final$fitted_g$Ulist[["ED_tPCA"]], tl.cex=0.5)

The second (“ED_PCA_2”) describes a pattern of sharing among brain tissues that is, interestingly, strongly anti-correlated with whole blood.

corrplot(m_final$fitted_g$Ulist[["ED_PCA_2"]], tl.cex=0.5)

Comments on FLASH fits

devtools::load_all("/Users/willwerscheid/GitHub/flashr/")

Loading flashr

fl_final <- list()
fl_final[[1]] <- readRDS(paste0(fpath, "OHF.rds"))
fl_final[[2]] <- readRDS(paste0(fpath, "top5.rds"))
fl_final[[3]] <- readRDS(paste0(fpath, "top10.rds"))
fl_final[[4]] <- readRDS(paste0(fpath, "top20.rds"))
fl_final[[5]] <- readRDS(paste0(fpath, "zero.rds"))

missing.tissues <- c(7, 8, 19, 20, 24, 25, 31, 34, 37)
gtex.colors <- read.table("https://github.com/stephenslab/gtexresults/blob/master/data/GTExColors.txt?raw=TRUE", sep = '\t', comment.char = '')[-missing.tissues, 2]

par(mar=c(1,1,1,1))
par(mfrow=c(2,1))
for (i in 46:47) {
  barplot(fl_final[[1]]$EL[, i], main=paste('Loading', i), las=2,
          cex.names = 0.4, col=as.character(gtex.colors), names="")
}

t <- readRDS(paste0(fpath, "t.rds"))
t_random <- data.frame(fit = as.character(1:5),
                       fixed.loadings = c("Canonical + OHF", "Canonical + Top 5 Zero", "Canonical + Top 10 Zero", "Canonical + Top 20 Zero", "Zero only"), 
                       number.of.factors = sapply(fl_final, function(x) {flash_get_nfactors(x)}),
                       minutes.to.fit = sapply(t$random, function(x) {as.numeric(x, units="mins")})
)
knitr::kable(t_random, digits=1)

par(mfrow = c(1, 2))
plot_gf.w <- function(fl, names, colors, main) {
  w <- 1 - sapply(fl$gf, function(g) {g$pi0})
  barplot(w, ylim=c(0, 1), ylab="w", xlab="factor/loading",
          cex.names=0.4, names=names, col=colors, main=main)
}
canonical.names <- c("EE", as.character(1:44))
c.col <- "blue"
d.col <- "red"
plot_gf.w(fl_final[[1]], 
          c(canonical.names, paste0("D", 1:2)),
          c(rep(c.col, 45), rep(d.col, 2)),
          main="Canonical + OHF")
plot_gf.w(fl_final[[2]], 
          c(canonical.names, paste0("D", 1:5)),
          c(rep(c.col, 45), rep(d.col, 5)),
          main="Canonical + Top 5")

plot_gf.w(fl_final[[3]], 
          c(canonical.names, paste0("D", 1:10)),
          c(rep(c.col, 45), rep(d.col, 10)),
          main="Canonical + Top 10")
plot_gf.w(fl_final[[4]], 
          c(canonical.names, paste0("D", 1:20)),
          c(rep(c.col, 45), rep(d.col, 20)),
          main="Canonical + Top 20")

plot_gf.w(fl_final[[5]], 
          paste0("D", 1:37),
          rep(d.col, 37),
          main="Zero")

par(mfrow = c(1, 2))
plot_gf.s <- function(fl, main) {
  s <- sqrt(1 / sapply(fl$gf, function(g) {g$a}))
  barplot(s[1:45], ylim=c(0, 6), ylab="sigma", xlab="factor/loading",
          cex.names=0.4, names=canonical.names, col="blue", main=main)
}
plot_gf.s(fl_final[[1]], 
          main="Canonical + OHF")
plot_gf.s(fl_final[[2]],
          main="Canonical + Top 5")

plot_gf.s(fl_final[[3]],
          main="Canonical + Top 10")
plot_gf.s(fl_final[[4]],
          main="Canonical + Top 20")

t_final <- data.frame(fit = as.character(1:5),
                      fixed.loadings = c("Canonical + OHF", "Canonical + Top 5 Zero", "Canonical + Top 10 Zero", "Canonical + Top 20 Zero", "Zero only"), 
                       number.of.factors = sapply(fl_final, function(x) {flash_get_nfactors(x)}),
                       minutes.to.fit = sapply(t$final, function(x) {as.numeric(x, units="mins")})
)
knitr::kable(t_final, digits=1)

par(mfrow = c(1, 2))
plot_pve <- function(fl, names, colors, main) {
  pve <- flash_get_pve(fl)
  barplot(pve[2:length(pve)], ylim=c(0, .05), ylab="PVE", 
          xlab="factor/loading", cex.names=0.4,
          names=names, col=colors, main=main)
}
canonical.names <- as.character(1:44)
plot_pve(fl_final[[1]], 
         c(canonical.names, paste0("D", 1:2)),
         c(rep(c.col, 44), rep(d.col, 2)),
         main="Canonical + OHF")
plot_pve(fl_final[[2]], 
         c(canonical.names, paste0("D", 1:5)),
         c(rep(c.col, 44), rep(d.col, 5)),
         main="Canonical + Top 5")

plot_pve(fl_final[[3]], 
         c(canonical.names, paste0("D", 1:10)),
         c(rep(c.col, 44), rep(d.col, 10)),
         main="Canonical + Top 10")
plot_pve(fl_final[[4]], 
         c(canonical.names, paste0("D", 1:20)),
         c(rep(c.col, 44), rep(d.col, 20)),
         main="Canonical + Top 20")

plot_pve(fl_final[[5]], 
         paste0("D", 2:37),
         rep(d.col, 36),
         main="Zero")

Total fitting time

The total time needed for each workflow is as follows.

fl_t <- (t_random$minutes.to.fit + t_final$minutes.to.fit
         + c(t$strong$OHF, t$strong$zero * c(5, 10, 20, 37) / 37))
total_time <- data.frame(method = c("OHF", "Top 5", "Top 10", 
                                    "Top 20", "Zero", "MASH"),
                         total.time.to.fit = c(fl_t, 5.3))
knitr::kable(total_time, digits=1)

MASH v FLASH initial observations

There is substantial disagreement among fits. The correlation matrix for posterior means is as follows.

all_fitted <- sapply(fl_final, flash_get_fitted_values)
all_fitted <- cbind(all_fitted, as.vector(t(get_pm(m_final))))
colnames(all_fitted) <- c("OHF", "Top 5", "Top 10", 
                          "Top 20", "Zero", "MASH")
round(cor(all_fitted), digits=3)

         OHF Top 5 Top 10 Top 20  Zero  MASH
OHF    1.000 0.984  0.983  0.983 0.971 0.966
Top 5  0.984 1.000  0.995  0.994 0.966 0.961
Top 10 0.983 0.995  1.000  0.999 0.968 0.963
Top 20 0.983 0.994  0.999  1.000 0.969 0.963
Zero   0.971 0.966  0.968  0.969 1.000 0.979
MASH   0.966 0.961  0.963  0.963 0.979 1.000

Entry \((i, j)\) in the following matrix gives, of the values that method \(j\) finds to be significant at 5%, the proportion of values that method \(i\) also finds to be significant at 5%. For example, of the values that MASH finds to be significant, OHF also finds 90% to be significant. But of the values that OHF finds to be significant, MASH only agrees in 81% of cases.

fl_lfsr <- readRDS(paste0(fpath, "fllfsr.rds"))

all_signif <- sapply(fl_lfsr, function(x) {as.numeric(x <= 0.05)})
all_signif <- cbind(all_signif, as.vector(t(get_lfsr(m_final) <= 0.05)))
m <- ncol(all_signif)
agree_mat <- matrix(0, nrow=m, ncol=m)
for (i in 1:m) {
  for (j in 1:m) {
    agree_mat[i, j] <- (sum(all_signif[, i] * all_signif[, j]) /
                          sum(all_signif[, j]))
  }
}
rownames(agree_mat) <- colnames(agree_mat) <- colnames(all_fitted)
round(agree_mat, digits=2)

        OHF Top 5 Top 10 Top 20 Zero MASH
OHF    1.00  0.93   0.92   0.92 0.90 0.90
Top 5  0.93  1.00   0.96   0.96 0.92 0.92
Top 10 0.92  0.96   1.00   0.98 0.93 0.92
Top 20 0.92  0.95   0.97   1.00 0.93 0.91
Zero   0.89  0.91   0.92   0.92 1.00 0.91
MASH   0.81  0.83   0.83   0.83 0.83 1.00

I will offer some explanations for the differences among fits in a subsequent analysis.

FLASH loadings

The FLASH loadings (as fitted on the strong tests) are as follows.

par(mar=c(1,1,1,1))
par(mfrow=c(3,2))
for(i in 1:flash_get_nfactors(fl_final[[5]])){
  barplot(fl_final[[5]]$EL[, i], main=paste('Loading', i), las=2,
          cex.names = 0.4, col=as.character(gtex.colors), names="")
}

Code

Click “Code” to view the code used to obtain the above results.

# Make sure to use branch "trackObj" when loading flashr.

# devtools::install_github("stephenslab/flashr", ref="trackObj")
devtools::load_all("/Users/willwerscheid/GitHub/flashr/")
# devtools::install_github("stephenslab/ebnm")
devtools::load_all("/Users/willwerscheid/GitHub/ebnm/")

library(mashr)

source("./code/fits.R")
source("./code/utils.R")

fpath <- "./output/MASHvFLASHgtex2/"

gtex <- readRDS(gzcon(url("https://github.com/stephenslab/gtexresults/blob/master/data/MatrixEQTLSumStats.Portable.Z.rds?raw=TRUE")))

strong <- t(gtex$strong.z)
random <- t(gtex$random.z)


# FLASH fits ------------------------------------------------------------

strong_data <- flash_set_data(strong, S = 1)
random_data <- flash_set_data(random, S = 1)

# methods to get DD: OHF; Zero(top5, top10, top20, full)
# LLs to use: canonical + OHF; Zero; canonical + top5/top10/top15; full

# Step 1. Learn data-driven loadings using strong tests.
t0 <- Sys.time()
fl_strong.OHF <- fit_flash_OHF(strong_data, Kmax=50, backfit=FALSE)
t_strong.OHF <- Sys.time() - t0

t0 <- Sys.time()
fl_strong.zero <- fit_flash_zero(strong_data, Kmax=50, backfit=FALSE)
t_strong.zero <- Sys.time() - t0

# Step 2. Fit the model to random tests to learn priors on factors.
LL <- list()
LL[[1]] <- fl_strong.OHF$f$EL

n <- nrow(strong)
canonical <- cbind(rep(1, n), diag(rep(1, n)))
LL[[2]] <- cbind(canonical, fl_strong.zero$f$EL[, 1:5])
LL[[3]] <- cbind(canonical, fl_strong.zero$f$EL[, 1:10])
LL[[4]] <- cbind(canonical, fl_strong.zero$f$EL[, 1:20])
LL[[5]] <- fl_strong.zero$f$EL

fl_random <- list()
t_random <- list()
for (i in 1:length(LL)) {
  t0 <- Sys.time()
  fl <- flash_add_fixed_l(random_data, LL[[i]])
  fl <- flash_backfit(random_data,
                      fl,
                      var_type="zero",
                      ebnm_fn = "ebnm_pn",
                      nullcheck = FALSE,
                      warmstart = TRUE,
                      verbose = TRUE)
  t_random[[i]] <- Sys.time() - t0
  fl_random[[i]] <- fl$f
}

# Step 3. Compute posteriors on strong tests, using priors from step 2.

fl_final <- list()
t_final <- list()
for (i in 1:length(LL)) {
  t0 <- Sys.time()
  fl <- flash_add_fixed_l(strong_data, LL[[i]])
  ebnm_param_f = lapply(fl_random[[i]]$gf, function(g) {list(g=g, fixg=TRUE)})
  fl <- flash_backfit(strong_data,
                      fl,
                      var_type="zero",
                      ebnm_fn = "ebnm_pn",
                      ebnm_param = list(f = ebnm_param_f, l = list()),
                      nullcheck = FALSE,
                      verbose = TRUE)
  t_final[[i]] <- Sys.time() - t0
  fl_final[[i]] <- fl$f
}

saveRDS(fl_final[[1]], paste0(fpath, "OHF.rds"))
saveRDS(fl_final[[2]], paste0(fpath, "top5.rds"))
saveRDS(fl_final[[3]], paste0(fpath, "top10.rds"))
saveRDS(fl_final[[4]], paste0(fpath, "top20.rds"))
saveRDS(fl_final[[5]], paste0(fpath, "zero.rds"))

t <- list()
t$strong <- list()
t$strong$OHF <- t_strong.OHF
t$strong$zero <- t_strong.zero
t$random <- t_random
t$final <- t_final


# Sample from posterior to get LFSR for FLASH fits ----------------------

nsamp <- 200
fl_lfsr <- list()
for (i in 1:length(fl_final)) {
  sampler <- flash_sampler(strong_data, fl_final[[i]], fixed="loadings")
  samp <- sampler(200)
  fl_lfsr[[i]] <- flash_lfsr(samp)
}

saveRDS(fl_lfsr, paste0(fpath, "fllfsr.rds"))


# MASH fit --------------------------------------------------------------

strong_data <- mash_set_data(t(strong), Shat = 1)
random_data <- mash_set_data(t(random), Shat = 1)

t$mash <- list()

# Step 1. Learn data-driven loadings using strong tests.
t0 <- Sys.time()
U.pca <- cov_pca(strong_data, 5)
U.ed <- cov_ed(strong_data, U.pca)
t$mash$ed <- Sys.time() - t0

# 2. Fit the model to random tests to learn mixture weights
t0 <- Sys.time()
U.c <- cov_canonical(random_data)
m_random <- mash(random_data, Ulist = c(U.ed,U.c))
t$mash$random <- Sys.time() - t0

# 3. Compute posterior summaries on the strong tests
t0 <- Sys.time()
m_final <- mash(strong_data, g=get_fitted_g(m_random), fixg=TRUE)
t$mash$final <- Sys.time() - t0

saveRDS(m_final, paste0(fpath, "m.rds"))
saveRDS(t, paste0(fpath, "t.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.6

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] flashr_0.5-12 corrplot_0.84 mashr_0.2-7   ashr_2.2-10  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17        pillar_1.2.1        compiler_3.4.3     
 [4] git2r_0.21.0        highr_0.6           plyr_1.8.4         
 [7] workflowr_1.0.1     R.methodsS3_1.7.1   R.utils_2.6.0      
[10] iterators_1.0.9     tools_3.4.3         testthat_2.0.0     
[13] digest_0.6.15       tibble_1.4.2        gtable_0.2.0       
[16] evaluate_0.10.1     memoise_1.1.0       lattice_0.20-35    
[19] rlang_0.2.0         Matrix_1.2-12       foreach_1.4.4      
[22] commonmark_1.4      yaml_2.1.17         parallel_3.4.3     
[25] ebnm_0.1-12         mvtnorm_1.0-7       xml2_1.2.0         
[28] roxygen2_6.0.1.9000 withr_2.1.1.9000    stringr_1.3.0      
[31] knitr_1.20          devtools_1.13.4     rprojroot_1.3-2    
[34] grid_3.4.3          R6_2.2.2            rmarkdown_1.8      
[37] rmeta_3.0           ggplot2_2.2.1       magrittr_1.5       
[40] whisker_0.3-2       scales_0.5.0        backports_1.1.2    
[43] codetools_0.2-15    htmltools_0.3.6     MASS_7.3-48        
[46] assertthat_0.2.0    softImpute_1.4      colorspace_1.3-2   
[49] stringi_1.1.6       lazyeval_0.2.1      munsell_0.4.3      
[52] doParallel_1.0.11   pscl_1.5.2          truncnorm_1.0-8    
[55] SQUAREM_2017.10-1   R.oo_1.21.0