Introduction to flashier

The matrix we analyze in this example is a 1,000 x 44 matrix of z-scores quantifying support for association between genetic variants in the human genome (the rows of the matrix) and gene expression measured across 44 human tissues (the columns):

library(flashier)
data(gtex)
nrow(gtex) 
ncol(gtex)
gtex[1:2,1:2]
#> [1] 1000
#> [1] 44
#>                                       Adipose_Subcutaneous
#> ENSG00000099977.9_22_24266954_A_C_b37             8.099352
#> ENSG00000233868.1_2_222475922_A_G_b37             1.079264
#>                                       Adipose_Visceral_Omentum
#> ENSG00000099977.9_22_24266954_A_C_b37                6.4129683
#> ENSG00000233868.1_2_222475922_A_G_b37               -0.7760486

The genetic variants are linked to genes, so each row of the matrix corresponds to a gene (the names of the rows are the genes’ Ensembl ids). These z-scores were originally calculated in Urbut et al. (2019) using data made available from the Genotype Tissue Expression (GTEx) project (Lonsdale et al. 2013), so we refer to this as the “GTEx data set.”

Our motivation for factorizing the GTEx data matrix is to identify patterns (“factors”) that give insight into the genetics of the 44 tissues: in particular, which tissues tend to act in concert, and to what degree? First, to get some intuition into the structure of this data set, we visualize the tissues in a 2-d embedding generated by t-SNE (L. Van der Maaten and Hinton 2008; L. J. P. Van der Maaten 2014; Krijthe 2015):

library(Rtsne)
library(ggrepel)
library(ggplot2)
library(cowplot)
set.seed(1)
out <- Rtsne(t(gtex), dims = 2, perplexity = 10)
pdat <- data.frame(d1 = out$Y[, 1],
                   d2 = out$Y[, 2],
                   tissue = colnames(gtex))
ggplot(pdat,aes(x = d1, y = d2, label = tissue)) +
  geom_point(size = 3, color = gtex_colors) +
  geom_text_repel(size = 2, max.overlaps = Inf) +
  labs(x = "t-SNE [1]", y = "t-SNE [2]") +
  theme_classic() +
  theme(axis.ticks = element_blank(),
        axis.text = element_blank())

Visually, one of the predominant trends is the clustering of brain tissues reflecting the greater genetic similarity among brain tissues. But there may be other interesting patterns that are revealed more effectively by a matrix factorization.

The flashier package implements the empirical Bayes matrix factorization (EBMF) framework developed by W. Wang and Stephens (2021). EBMF is based on the following model of an n × p data matrix ${\bf X}$:

$$ \mathbf{X} = \mathbf{L} \mathbf{F}^\top + \mathbf{E} \\ e_{ij} \sim \mathcal{N}(0, \sigma^2) \\ \ell_{ik} \sim g_\ell^{(k)} \in \mathcal{G}_{\ell} \\ f_{jk} \sim g_f^{(k)} \in \mathcal{G}_f, $$

where ${\bf L}, {\bf F}, {\bf E}$ are, respectively, matrices of dimension n × K, p × K and n × p storing real-valued elements lik, fjk, and eij. The matrices ${\bf L}$ and ${\bf F}$ form an approximate low-rank representation of ${\bf X}$, ${\bf X} \approx {\bf L} {\bf F}^\top$, and this low-rank representation is learned from the data. K specifies the rank of the reduced representation, and is typically set to a positive number much smaller than n and p.

The flexibility of the EBMF framework comes from its ability to accommodate, in principle, any choice of prior family for 𝒢ℓ and 𝒢f. As we illustrate below, different choices of prior families can lead to very different factorizations. Some choices closely correspond to existing matrix factorization methods. For example, if 𝒢ℓ and 𝒢f are both the families of zero-centered normal priors, then the low-rank representation ${\bf L}{\bf F}^\top$ is expected to be similar to a truncated singular value decomposition (SVD) (Nakajima and Sugiyama 2011). When point-normal prior families are used instead, one obtains empirical Bayes versions of sparse SVD or sparse factor analysis (Engelhardt 2010; Yang, Ma, and Buja 2014; Witten, Tibshirani, and Hastie 2009). Thus, EBMF is a highly flexible framework for matrix factorization that includes important previous methods as special cases, but also many new combinations.

Since fitting the EBMF model can be reduced to solving a sequence of EBNM problems (W. Wang and Stephens 2021), flashier leverages ebnm for the core model fitting routines. In brief, fitting an EBMF model involves solving an EBNM problem separately for each column of ${\bf L}$ and for each column of ${\bf F}$; the former results in a fitted prior ĝℓ(k) and posterior estimates of entries ℓik, and the latter results in a fitted prior ĝf(k) and posteror estimates of entries fjk. These EBNM models are repeatedly re-fitted until the priors and posterior estimates converge to a fixed point. Therefore, fitting an EBMF model can involve solving many hundreds or even many thousands of EBNM problems. Leveraging the efficient EBNM solvers available in ebnm makes it possible for flashier to tackle large-scale matrix factorization problems.

The flashier function flash() provides the main interface for fitting EBMF models. The following call to flash() fits a matrix factorization to the GTEx data with normal priors on ${\bf L}$ and ${\bf F}$:

flash_n <- flash(gtex, ebnm_fn = ebnm_normal, backfit = TRUE)
#> Adding factor 1 to flash object...
#> Adding factor 2 to flash object...
#> Adding factor 3 to flash object...
#> Adding factor 4 to flash object...
#> Adding factor 5 to flash object...
#> Adding factor 6 to flash object...
#> Adding factor 7 to flash object...
#> Adding factor 8 to flash object...
#> Adding factor 9 to flash object...
#> Adding factor 10 to flash object...
#> Adding factor 11 to flash object...
#> Adding factor 12 to flash object...
#> Adding factor 13 to flash object...
#> Adding factor 14 to flash object...
#> Adding factor 15 to flash object...
#> Factor doesn't significantly increase objective and won't be added.
#> Wrapping up...
#> Done.
#> Backfitting 14 factors (tolerance: 6.56e-04)...
#>   Difference between iterations is within 1.0e+02...
#>   Difference between iterations is within 1.0e+01...
#>   Difference between iterations is within 1.0e+00...
#>   Difference between iterations is within 1.0e-01...
#>   Difference between iterations is within 1.0e-02...
#> Wrapping up...
#> Done.
#> Nullchecking 14 factors...
#> Done.

The “ebnm_fn” argument specifies the prior family for ${\bf L}$ and ${\bf F}$, and accepts any of the prior family functions implemented in the ebnm package (e.g., ebnm_point_normal or ebnm_unimodal; see Table X for a near-complete list of prior families). Specifying the number of factors K is not needed because flashier automatically tries to determine an appropriate rank by checking whether the addition of new factors substantially improves the fit. (If one prefers a smaller number of factors than what is determined automatically, the maximum number of factors can be adjusted with the “greedy_Kmax” argument.) Here we also turned on backfitting by setting backfit = TRUE. This is generally recommended because it improves the quality of the fit (W. Wang and Stephens 2021), with the caveat that it may make the model fitting too slow for very large data sets. Even with backfitting, the full computation here is very fast; when running this example on a current MacBook Pro, for example, the model fitting is completed in less than one second.

The plot() method for flash objects can be used to visualize the estimated matrix ${\bf L}$. Adding color to distinguish among tissues yields an effective visualization that aids interpretation of factors:

plot(flash_n, pm_which = "factors", pm_colors = gtex_colors, plot_type = "bar")

The first factor sets the “baseline” z-score for each tissue. Other factors appear to capture expected trends: for example, factor 2 picks up a z-score pattern that is more dominant in brain tissues. More tissue-specific patterns are picked up by factor 3, which captures whole blood effects, and factor 6, which captures effects more specific to testis.

SVD and SVD-like matrix factorizations may provide compact approximations of a data matrix, but they are generally known to have poor interpretability. “Sparse” factorizations can produce factors that are much more individually interpretable. In the EBNM framework, a sparse matrix factorization can be achieved by choosing a flexible prior family that better encourages sparsity in ${\bf L}$ and/or ${\bf F}$. The flashier interface makes it straightforward to experiment with different prior families; for example, we can use point-normal priors by simply modifying the “ebnm_fn” argument:

flash_pn <- flash(gtex, ebnm_fn = ebnm_point_normal, backfit = TRUE)
plot(flash_pn, pm_which = "factors", pm_colors = gtex_colors, plot_type = "bar")

#> Adding factor 1 to flash object...
#> Adding factor 2 to flash object...
#> Adding factor 3 to flash object...
#> Adding factor 4 to flash object...
#> Adding factor 5 to flash object...
#> Adding factor 6 to flash object...
#> Adding factor 7 to flash object...
#> Adding factor 8 to flash object...
#> Adding factor 9 to flash object...
#> Adding factor 10 to flash object...
#> Adding factor 11 to flash object...
#> Adding factor 12 to flash object...
#> Adding factor 13 to flash object...
#> Adding factor 14 to flash object...
#> Adding factor 15 to flash object...
#> Adding factor 16 to flash object...
#> Adding factor 17 to flash object...
#> Adding factor 18 to flash object...
#> Adding factor 19 to flash object...
#> Adding factor 20 to flash object...
#> Adding factor 21 to flash object...
#> Adding factor 22 to flash object...
#> Factor doesn't significantly increase objective and won't be added.
#> Wrapping up...
#> Done.
#> Backfitting 21 factors (tolerance: 6.56e-04)...
#>   Difference between iterations is within 1.0e+02...
#>   Difference between iterations is within 1.0e+01...
#>   Difference between iterations is within 1.0e+00...
#>   Difference between iterations is within 1.0e-01...
#>   Difference between iterations is within 1.0e-02...
#>   Difference between iterations is within 1.0e-03...
#> Wrapping up...
#> Done.
#> Nullchecking 21 factors...
#> Done.

The point-normal priors require estimation of slightly more parameters than the normal priors, but the call to flash() is still fast; it took less than 5 seconds to run on our MacBook Pro.

As before, the first factor acts as a “baseline” and the second factor captures z-scores largely specific to brain tissues. However, there are also some major differences: some factors are much more clearly tissue-specific (e.g., factors 3 and 4 for, respectively, whole blood and testis); other patterns include cerebellar-specific patterns (factor 9) and heart-specific patterns (factor 12). These same patterns are also captured in the matrix factorization with normal priors, but they often appear in combination with other patterns in unintuitive ways.

Another matrix factorization approach that has been recently proposed is semi-nonnegative matrix factorization (Ding, Li, and Jordan 2010; M. Wang, Fischer, and Song 2019; He et al. 2020), in which the elements of either ${\bf L}$ or ${\bf F}$ (but not both) are constrained to be nonnegative. Supposing that ${\bf L}$ is constrained to be nonnegative, the elements ℓik ≥ 0 can be viewed as “weights” or “memberships,” and each row of ${\bf X}$ is then interpretable as a weighted combination of factors: $x_{i \cdot} \approx \sum_{k=1}^K \ell_{ik} f_{\cdot k}$.

To obtain a semi-nonnegative matrix factorization via flashier, one need only choose appropriate prior families. Here we use point-normal priors for ${\bf L}$ and point-exponential priors for ${\bf F}$; thus we combine the benefits of both semi-nonnegative and sparse matrix factorization. In flashier, we can specify different priors for ${\bf L}$ and ${\bf F}$ by setting the “ebnm_fn” argument to be a list in which the first element is the prior family for ${\bf L}$ and the second the prior family for ${\bf F}$:

flash_snn <- flash(gtex, 
                   ebnm_fn = c(ebnm_point_normal, ebnm_point_exponential), 
                   backfit = TRUE)
plot(flash_snn, pm_which = "factors", pm_colors = gtex_colors, plot_type = "bar")

#> Adding factor 1 to flash object...
#> Adding factor 2 to flash object...
#> Adding factor 3 to flash object...
#> Adding factor 4 to flash object...
#> Adding factor 5 to flash object...
#> Adding factor 6 to flash object...
#> Adding factor 7 to flash object...
#> Adding factor 8 to flash object...
#> Adding factor 9 to flash object...
#> Adding factor 10 to flash object...
#> Adding factor 11 to flash object...
#> Adding factor 12 to flash object...
#> Adding factor 13 to flash object...
#> Adding factor 14 to flash object...
#> Adding factor 15 to flash object...
#> Adding factor 16 to flash object...
#> Adding factor 17 to flash object...
#> Adding factor 18 to flash object...
#> Adding factor 19 to flash object...
#> Adding factor 20 to flash object...
#> Adding factor 21 to flash object...
#> Adding factor 22 to flash object...
#> Factor doesn't significantly increase objective and won't be added.
#> Wrapping up...
#> Done.
#> Backfitting 21 factors (tolerance: 6.56e-04)...
#>   Difference between iterations is within 1.0e+02...
#>   Difference between iterations is within 1.0e+01...
#>   Difference between iterations is within 1.0e+00...
#>   Difference between iterations is within 1.0e-01...
#>   Difference between iterations is within 1.0e-02...
#>   Difference between iterations is within 1.0e-03...
#> Wrapping up...
#> Done.
#> Nullchecking 21 factors...
#> Done.

This code returned results in less than 3 seconds on our MacBook Pro.

The result is in many ways strikingly similar to the sparse factorization obtained using point-normal priors. In some cases, the same tissue-specific effects are recovered, up to a difference in signs (e.g., factor 5 for fibroblasts and factor 8 for thyroid). But the semi-nonnegative factorization identifies factors that are arguably more intuitive; for example, point-normal factor 10, which combines artery and adipose tissues, is split into two factors in the semi-nonnegative factorization (factor 9 for artery and factor 11 for adipose tissue).

The flashier package has many more options which we did not explore here, including other choices of prior family, different options for estimating the variances of the residuals ${\bf E}$, and options for fine-tuning the model fitting to improve the speed and quality of EBMF model fits for larger data sets. Many of these options are discussed in the other vignettes.

Session information

The following R version and packages were used to generate this vignette:

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
#>  [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       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggrepel_0.9.6   Rtsne_0.17      dplyr_1.1.4     cowplot_1.1.3  
#> [5] ggplot2_3.5.1   flashier_1.0.54 ebnm_1.1-34     rmarkdown_2.29 
#> 
#> loaded via a namespace (and not attached):
#>  [1] softImpute_1.4-1     gtable_0.3.6         xfun_0.49           
#>  [4] bslib_0.8.0          htmlwidgets_1.6.4    lattice_0.22-6      
#>  [7] quadprog_1.5-8       vctrs_0.6.5          tools_4.4.2         
#> [10] generics_0.1.3       parallel_4.4.2       Polychrome_1.5.1    
#> [13] tibble_3.2.1         fansi_1.0.6          pkgconfig_2.0.3     
#> [16] Matrix_1.7-1         data.table_1.16.2    SQUAREM_2021.1      
#> [19] RColorBrewer_1.1-3   RcppParallel_5.1.9   scatterplot3d_0.3-44
#> [22] lifecycle_1.0.4      truncnorm_1.0-9      farver_2.1.2        
#> [25] compiler_4.4.2       progress_1.2.3       munsell_0.5.1       
#> [28] RhpcBLASctl_0.23-42  htmltools_0.5.8.1    sys_3.4.3           
#> [31] buildtools_1.0.0     sass_0.4.9           yaml_2.3.10         
#> [34] lazyeval_0.2.2       plotly_4.10.4        crayon_1.5.3        
#> [37] tidyr_1.3.1          pillar_1.9.0         jquerylib_0.1.4     
#> [40] uwot_0.2.2           cachem_1.1.0         trust_0.1-8         
#> [43] gtools_3.9.5         tidyselect_1.2.1     digest_0.6.37       
#> [46] purrr_1.0.2          ashr_2.2-63          labeling_0.4.3      
#> [49] maketools_1.3.1      splines_4.4.2        fastmap_1.2.0       
#> [52] grid_4.4.2           colorspace_2.1-1     cli_3.6.3           
#> [55] invgamma_1.1         magrittr_2.0.3       utf8_1.2.4          
#> [58] withr_3.0.2          prettyunits_1.2.0    scales_1.3.0        
#> [61] horseshoe_0.2.0      httr_1.4.7           fastTopics_0.6-192  
#> [64] deconvolveR_1.2-1    hms_1.1.3            pbapply_1.7-2       
#> [67] evaluate_1.0.1       knitr_1.49           viridisLite_0.4.2   
#> [70] irlba_2.3.5.1        rlang_1.1.4          Rcpp_1.0.13-1       
#> [73] mixsqp_0.3-54        glue_1.8.0           jsonlite_1.8.9      
#> [76] R6_2.5.1

Bibliography

Ding, Chris H. Q., Tao Li, and Michael I. Jordan. 2010. “Convex and Semi-Nonnegative Matrix Factorizations.” IEEE Transactions on Pattern Analysis and Machine Intelligence 32 (1): 45–55.
Engelhardt, Matthew, Barbara E. AND Stephens. 2010. “Analysis of Population Structure: A Unifying Framework and Novel Methods Based on Sparse Factor Analysis.” PLOS Genetics 6 (9): 1–12.
He, Yuan, Surya B. Chhetri, Marios Arvanitis, Kaushik Srinivasan, François Aguet, Kristin G. Ardlie, Alvaro N. Barbeira, et al. 2020. “sn-spMF: Matrix Factorization Informs Tissue-Specific Genetic Regulation of Gene Expression.” Genome Biology 21: 235.
Krijthe, Jesse H. 2015. : t-Distributed Stochastic Neighbor Embedding Using a Barnes-Hut Implementation. https://github.com/jkrijthe/Rtsne.
Lonsdale, John, Jeffrey Thomas, Mike Salvatore, et al. 2013. “The Genotype-Tissue Expression (GTEx) Project.” Nature Genetics 45 (6): 580–85.
Nakajima, Shinichi, and Masashi Sugiyama. 2011. “Theoretical Analysis of Bayesian Matrix Factorization.” Journal of Machine Learning Research 12: 2583–2648.
Urbut, Sarah M., Gao Wang, Peter Carbonetto, and Matthew Stephens. 2019. “Flexible Statistical Methods for Estimating and Testing Effects in Genomic Studies with Multiple Conditions.” Nature Genetics 51: 187–95.
Van der Maaten, L. J. P. 2014. “Accelerating t-SNE Using Tree-Based Algorithms.” Journal of Machine Learning Research 15: 3221–45.
Van der Maaten, Laurens, and Geoffrey Hinton. 2008. “Visualizing Data Using t-SNE.” Journal of Machine Learning Research 9 (11): 2579–2605.
Wang, Miaoyan, Jonathan Fischer, and Yun S Song. 2019. “Three-Way Clustering of Multi-Tissue Multi-Individual Gene Expression Data Using Semi-Nonnegative Tensor Decomposition.” The Annals of Applied Statistics 13 (2): 1103–27.
Wang, Wei, and Matthew Stephens. 2021. “Empirical Bayes Matrix Factorization.” Journal of Machine Learning Research 22 (120): 1–40.
Witten, Daniela M, Robert Tibshirani, and Trevor Hastie. 2009. “A Penalized Matrix Decomposition, with Applications to Sparse Principal Components and Canonical Correlation Analysis.” Biostatistics 10 (3): 515–34.
Yang, Dan, Zongming Ma, and Andreas Buja. 2014. “A Sparse Singular Value Decomposition Method for High-Dimensional Data.” Journal of Computational and Graphical Statistics 23 (4): 923–42.