# Encoding Domain Information with Sparse Priors for Inferring Explainable
Latent Variables

Arber Qoku  German Cancer Consortium (DKTK) and German Cancer Research Center
(DKFZ)  Goethe University Frankfurt, Germany  Florian Buettner  German Cancer
Consortium (DKTK) and German Cancer Research Center (DKFZ)  Goethe University
Frankfurt, Germany  Frankfurt Cancer Insititute, Germany
{arber.qoku<EMAIL_ADDRESS>
###### Abstract

Latent variable models are powerful statistical tools that can uncover
relevant variation between patients or cells, by inferring unobserved hidden
states from observable high-dimensional data. A major shortcoming of current
methods, however, is their inability to learn sparse and interpretable hidden
states. Additionally, in settings where partial knowledge on the latent
structure of the data is readily available, a statistically sound integration
of prior information into current methods is challenging. To address these
issues, we propose spex-LVM, a factorial latent variable model with sparse
priors to encourage the inference of explainable factors driven by domain-
relevant information. spex-LVM utilizes existing knowledge of curated
biomedical pathways to automatically assign annotated attributes to latent
factors, yielding interpretable results tailored to the corresponding domain
of interest. Evaluations on simulated and real single-cell RNA-seq datasets
demonstrate that our model robustly identifies relevant structure in an
inherently explainable manner, distinguishes technical noise from sources of
biomedical variation, and provides dataset-specific adaptations of existing
pathway annotations. Implementation is available at https://github.com/MLO-
lab/spexlvm.

##  1 Introduction

Recent advances in omics profiling technologies from proteomics to single-cell
RNA-sequencing (scRNA-seq) provide an ever-finer molecular characterization of
patients. In particular, scRNA-seq technologies facilitate the investigation
of biological processes at a cellular level. Such a significant increase in
data resolution paints a much more comprehensive picture of the underlying
drivers of heterogeneity not only between patients, but also at the inter-
cellular level. This allows the discovery of new patient sub-populations as
well as novel cell sub-populations present across patients, which in turn has
the potential to accelerate progress in personalized medicine. However,
extracting biologically and clinically meaningful structures from thousands or
millions of cells and hundreds of patients presents an array of computational
challenges[Lähnemann et al., 2020]. Latent variable modeling is a general
approach suitable for finding major axes of variation between samples in terms
of unobserved numerical components. However, understanding the output of such
models usually involves an ongoing exchange between data analysts and domain
experts, which can be inflexible and time-consuming. A much more efficient
alternative is taking advantage of established domain knowledge to
automatically infer interpretable latent factors, thereby facilitating a
faster transition to conclusive insights. Biomedical gene set annotations, for
instance, serve as a reliable source of information for encoding a relevant
structure of the latent factors, providing an additional layer of
understanding to numerically encoded outputs.

Figure 1: Schematic overview of spex-LVM. A gene expression matrix is
decomposed as a sum of latent factors informed by potentially noisy gene set
annotations. Spex-LVM learns the true underlying axes of variation by pruning
redundant factors and refining noisy annotations. On the right we depict a
typical downstream analysis involving imputation of incomplete observations,
factor ranking and inspection of refined features, mapping of samples onto
latent axes for visualization and modeling of technical noise via additional
unannotated factors.

We propose spex-LVM, a flexible and efficient latent variable model capable of
utilizing prior information of well-known biomedical gene set annotations to
guide the inference of explainable sources of variation in scRNA-seq data.
More specifically, spex-LVM learns interpretable factors that capture
meaningful subpopulations of cells or patients, automatically tagged with the
biomedical pathway characterizing the corresponding axis of variation. We
account for unwanted technical noise[Hicks et al., 2018] by explicitly
modeling irrelevant sources of variation as additional factors that capture
perturbations affecting large groups of features. In addition, our approach
considers potentially noisy annotations and provides refined alternatives
supported by the data. Finally, our method provides identifiable latent
factors. Experiments show that the order of the learned factors matches the
order of the provided annotations, thereby circumventing a known drawback
inherent to the conventional factor analysis[Shapiro, 1985], and yielding
consistent and reproducible results. We should note that our method is not
limited to biomedical data analysis, and it can be easily applied to other
domains where partial information on the latent structure is available.

##  2 Method

Spex-LVM decomposes high dimensional gene expression data into distinct
factors that jointly contribute to the data generation process. One major
advantage of our model is its ability to incorporate existing knowledge of
biomedical pathways to encourage learning of interpretable factors as depicted
schematically in Figure 1. Furthermore, the model assumes incomplete prior
knowledge by modeling additional unannotated sparse and dense factors, which
serve as placeholders for encoding sparse unidentified pathways, and external
or technical perturbations that affect the expression pattern of nearly all
genes.
Let $\mathbf{Y}\in\mathbb{R}^{N\times G}$ denote the gene expression matrix
comprising $N$ samples across $G$ genes. spex-LVM decomposes $\mathbf{Y}$ into
a sum of $A$ annotated and $H$ unannotated factors, where the latter can be
further differentiated between sparse and dense unannotated factors.

$\mathbf{Y}=\sum_{a}^{A}\mathbf{p}_{a}\otimes\mathbf{v}_{a}+\sum_{h}^{H}\mathbf{q}_{h}\otimes\mathbf{u}_{h}+\mathbf{\Psi},$
(1)

where $\mathbf{p}_{a}$, $\mathbf{q}_{h}\in\mathbb{R}^{N}$ describe the state
of the annotated factor $a$ and the unannotated factor $h$ across $N$ cells,
and $\mathbf{v}_{a}$, $\mathbf{u}_{h}\in\mathbb{R}^{G}$ are the corresponding
regulatory weights or factor loadings for factors $a$ and $h$ across $G$ genes
in $\mathbf{Y}$. The residuals are denoted by
$\mathbf{\Psi}=diag(\bm{\sigma}^{2})$, a diagonal matrix storing the marginal
variances $\sigma_{g}^{2}$ of each gene $g$. We may rewrite 1 concisely in
matrix notation by defining $\mathbf{X}\in\mathbb{R}^{N\times K}$ as
$\left[\mathbf{p}_{1},\cdots,\mathbf{p}_{A},\mathbf{q}_{1},\cdots\mathbf{q}_{H}\right]$
and $\mathbf{W}\in\mathbb{R}^{G\times K}$ as
$\left[\mathbf{v}_{1},\cdots,\mathbf{v}_{A},\mathbf{u}_{1},\cdots\mathbf{u}_{H}\right]$.

$\mathbf{Y}=\mathbf{X}\cdot\mathbf{W}^{T}+\mathbf{\Psi}$ (2)

Factor analysis models[Bishop, 2006, Murphy, 2012] typically assume the
projection of a sample $x_{i,k}\sim\mathcal{N}(0,1)$ in the latent space, to a
sample $y_{i,g}\sim\mathcal{N}(\sum_{k}w_{g,k}x_{i,k},\sigma_{g}^{2})$ in the
original high-dimensional space via the corresponding factor loadings
$w_{g,k}\sim\mathcal{N}(0,1)$. We extend the conventional factor analysis by
introducing a _global-local_ shrinkage prior[Polson and Scott, 2010] on the
scales of each weight $w_{g,k}$.

$w_{g,k}|\tau,\lambda_{g,k}\sim\mathcal{N}(0,\tau^{2}\cdot\lambda_{g,k}^{2})$
(3)

The consequences of this addition are twofold. First, the global scale $\tau$
shrinks all the weights towards zero, encouraging the model to learn sparse
representations of the underlying factors. This is especially useful when
learning sparse unannotated factors. Second, the local scales act as a
regulator on each individual weight. A small local scale $\lambda_{g,k}$
confines $w_{g,k}$ to values very close to zero. On the other hand, a large
local scale allows its corresponding parameter to escape the global
regularization penalty. As an enhancement to the original horseshoe
prior[Carvalho et al., 2009], the regularized version provides means to
control the slab width of each individual scale[Piironen et al., 2017],
allowing partial knowledge to guide the model, e.g. annotated gene sets.

$\displaystyle w_{g,k}|\tau,\lambda_{g,k},c_{g,k}$
$\displaystyle\sim\mathcal{N}(0,\tau^{2}\cdot\tilde{\lambda}_{g,k}^{2}),$ (4)
$\displaystyle\tilde{\lambda}_{g,k}$
$\displaystyle=\frac{c_{g,k}^{2}\lambda_{g,k}^{2}}{c_{g,k}^{2}+\tau^{2}\lambda_{g,k}^{2}}$
(5)

The slab width effectively encodes our confidence about the given gene set
annotations. We accommodate potentially noisy annotations by providing soft
bounds for the local scales, allowing the adaptation of incorrect signals
given sufficient evidence from the data. Experiments show that a slab width of
$c_{g,k}\approx 1.0$ for genes annotated to a pathway and $0<c_{g,k}<0.1$ for
unannotated genes provide the best results. We assign a relatively large slab
width across all genes for the unannotated dense factors, and vice-versa for
the unannotated sparse factors, allowing the model to autoregulate active
signals based on the data. Alternative implementations of sparse factor
analysis models rely on the spike-and-slab prior [Buettner et al., 2017,
Mitchell and Beauchamp, 1988]; however, due to their discrete nature, these
models cannot be trained using blackbox stochastic variational inference and
therefore tend to be inefficient. We implement the alternative
parameterization of the regularized horseshoe as suggested in [Piironen et
al., 2017] and described above. The resulting algorithm is amenable to
efficient and scalable stochastic variational inference[Ranganath et al.,
2014, Hoffman et al., 2013, Wingate and Weber, 2013], and optimizes the
evidence lower bound via standard gradient methods such as Adam[Kingma and Ba,
2014].

##  3 Results and Discussion

We first evaluate several properties of spex-LVM on synthetic data. We show
that spex-LVM learns informed and relevant factors, while redundant annotated
factors that do not explain variation in a specific dataset are pruned.
Second, we illustrate that spex-LVM is robust against noisy gene set
annotations. Third, we demonstrate that spex-LVM scales sub-linearly with the
number of samples and handles millions of samples and tens of thousands of
genes.
We then show on three real-world datasets that spex-LVM explicitly models and
separates unwanted variation from identified pathway factors that are
biologically meaningful. For each experiment, we derive gene set annotations
from established public databases such as MSigDB[Liberzon et al., 2015] and
Reactome[Fabregat et al., 2018].

###  3.1 Synthetic Data

Figure 2: Learning from noisy information. Models are trained independently on
different levels of deliberate noise. Each line summarizes the progress of the
$F_{1}$ score computed from the true annotation signals before introducing
noise and the learned signals. Below, a subset of the learned weights after
specific iterations as heatmaps. Colors encode absolute gene weight (black:
inactive, white: active). Noisy factor loadings are indicated with colored
edges (red: false positives, green: false negatives).

We assess the ability of our model to discern valuable prior information from
noisy annotations. We simulate a gene expression toy dataset of $N=$ 10,000
cells and $G=$ 8,000 genes as a linear transformation of $K=100$ latent
factors generated from a standard normal distribution
$\mathcal{N}(0,\mathbf{I})$, where $\mathbf{I}$ denotes the $K$-dimensional
identity matrix. The factor loadings were generated from a zero-centered
normal distribution $\mathcal{N}(0,4)$, where loadings with an absolute value
of less than $0.5$ were set to zero, to emphasize the gap between active and
inactive signals. In addition, we randomly set $85\%$-$95\%$ of the loadings
to zero to more closely resemble the true size of gene set collections. After
extracting the true factor mask, indicating active and inactive genes for each
factor, we introduce additional noise in the form of false positive, false
negative annotations and redundant factors. Figure 2 shows a summary of the
results for several models trained independently on different noise fractions.
A noise fraction of $0.1$ means that $10\%$ of true active genes have been
switched off, i.e. false negatives, as well as $10\%$ of true inactive genes
have been incorrectly switched on, i.e. false positives. For each noise level
we train four different models with $0\%$, $10\%$, $20\%$ and $50\%$
additional random factors, irrelevant to the data generation process. The
validation on the toy data demonstrates that spex-LVM learns relevant
structure from noisy information and converges to the true latent factors. We
identify three distinct stages during model training (Fig. 2). In stage one,
the model eagerly discovers all of the factors provided by the annotation
mask, as indicated by multiple false positive factor loadings illustrated with
red colored edges in the first heatmap. The redundant factors 6 and 7 are
weakly identified, and quickly vanish during the second learning stage. We
call the third stage the denoising stage as the model iteratively flips false
positives and false negatives into their correct state. The fourth and last
heatmap in Fig. 2 shows the state of the weights after the model has fully
converged.
We conclude the synthetic evaluation by providing a summary of the running
times (Fig. 3) recorded across different dataset sizes with respect to the
number of samples, features and latent factors. While sweeping one parameter,
the default values for the other parameters are as follows: $N=$ 10,000
samples, $D=$ 8,000 features and $K=$ 100 factors. Models were trained on a
single NVIDIA Quadro RTX 5000 GPU with 16 GB memory.

Figure 3: Scalability of spex-LVM in terms of sample, feature and factor size.
The $y$-axes match on all plots to help gauge the relative impact of each data
dimension.

The available memory can fit up to 20 thousand samples by approximately 10
thousand features concurrently. The significant jump in the elapsed time for
the data with more than 20 thousand samples is due to the introduction of
batches in the training process. Consequently, the number of samples presents
the largest impact in the amount of time required for the algorithm to
converge. However, for more than 100 thousand samples, we see a sub-linear
increase in the runtime as the number of the total batch updates necessary for
convergence does not increase as rapidly with the sample size. On the other
hand, due to the parallelized computation, spex-LVM efficiently handles up to
20 thousand features and fit hundreds or even thousands of latent factors for
tens of thousands of samples in less than 30 minutes on average.

###  3.2 Mouse Embryonic Stem Cells

Figure 4: Bivariate visualization of mouse embryonic stem cells in terms of
the two main factors, the G2/M checkpoint (H) and the cell cycle (R). Each
cell is color-coded according to its corresponding stage in the cell cycle.

We evaluate our model on 182 mouse embryonic stem cells transcriptomes
measured across 6 thousand genes, where the stage of the cell cycle is
determined a priori and provided for each cell[Buettner et al., 2015]. During
training, we provide annotations from 29 hallmark gene sets and 54 reactome
gene sets with a median size of 68 genes. As expected, the model robustly
identifies biological pathways responsible for the cell cycle. Cell cycle (R),
G2/M checkpoint (H) and P53 pathway (H) comprise three of the top five
candidates in terms of variance explained. Figure 4 depicts the cell
population mapped onto the latent axes corresponding to the G2/M checkpoint
and the cell cycle. While both factors contribute in finding sub-populations
of cells during the same cell cycle phase, the G2/M checkpoint is a stronger
indicator of cells in the G2/M phase as supported by their corresponding area
under the receiver operating characteristic curve (_AUROC_) of 0.76 and 0.82,
respectively.

###  3.3 PBMCs from Lupus Patients

Next, we investigate a much larger dataset of peripheral blood mononuclear
cells from lupus patients profiled from droplet single-cell RNA-sequencing
(dscRNA-seq)[Kang et al., 2018], where cells were either stimulated with
Interferon-$\beta$ or left untreated, yielding two similarly sized groups of
cells. We adopt the preprocessing steps from [Rybakov et al., 2020], and train
our model on 13,576 samples, each comprising 979 features or genes. As prior
information we consider annotations provided by the reactome gene sets, where
at least 15 genes were present in the data, resulting in 131 annotations with
a median size of 25. In Figure 5 we summarize the top 10 most relevant factors
based by their corresponding relevance in terms of the coefficient of
determination of each factor for the complete dataset.

Figure 5: PBMCs from lupus patients results. Top 10 most relevant factors
discovered by spex-LVM (left). Two interferon signaling related factors
capture the differences between the expression profiles of the control group
and IFN-$\beta$ stimulated group (right).

We mainly observe pathways related to immune system and interferon response.
In particular, interferon signaling and interferon alpha beta signaling
capture the transcriptomic differences between the control and IFN-$\beta$
stimulated group, as anticipated. Pathways involved in the immune response
explain the main drivers of variation along the different cell types as shown
in Figure 6.

Figure 6: Relevant factors related to the immune system response. Factors of
general pathways, e.g. innate immune system, neutrophil degranulation, explain
variation within the larger clusters of monocytes (left). More specific
pathways infer axes of variation for smaller clusters such as the dendritic
cells and the B cells (right).

The two larger clusters of monocytes, CD14 and CD16, which comprise nearly
$40\%$ of the samples, exhibit significant differences from the rest of the
population in terms of the innate immune system and the neutrophil
degranulation. In addition, we observe relatively smaller, yet more specific
pathways, which explain less represented cell types. MHC class II antigen
presentation describes the process of binding and transporting antigenic
peptides to the membrane of professional antigen-presenting cells (APCs), such
as dendritic cells[Roche and Furuta, 2015]. Our experiments show that the
corresponding factor reliably separates dendritic cells from the other cell
types with an _AUROC_ score of $0.95$. Another candidate is the signaling by
the B cell receptor, consisting of only 19 gene set annotations, which
consistently infers the axis of variation along the B cells ($0.87$ _AUROC_
score).

###  3.4 Adult Human Retina

Finally, we investigate a transcriptome atlas of the human retina[Lukowski et
al., 2019]. The dataset consists of 20 thousand cells collected from the
neural retina layers of three different donors, and belong to seven
transcriptionally distinct clusters. We retain the top 8 thousand most
variable features during training. Due to the increased feature size of the
dataset, we derive 35 annotations from hallmark and 161 from reactome gene set
collections. Additionally, we allow the model to infer 2 unannotated sparse
and 2 dense factors, bringing the total number of latent factors to 200. The
results in Figure 7 reveal several insights in accordance with previous
findings.

Figure 7: Adult human retina atlas results. Top 10 most relevant factors
discovered by spex-LVM (left). Sparse factor 1 explains unwanted variation
characterized by upregulated low-quality _MALAT1_ reads (right).

First, the unannotated factors clearly explain the majority of the variance in
the data as supported by their accumulative $R^{2}$ score of approximately
$13\%$, which we attribute to potential sources of technical heterogeneity.
Sparse factor 1 represents unwanted technical variation in form of an
upregulation for low-quality _MALAT1_ and mitochondrial reads, as it consists
of a highly expressed _MALAT1_ gene followed by amplified signals of
mitochondrial-encoded genes (_MT-ATP6_ , _MT-CYB_ , _MT-ND4_ , _MT-
CO3_)[Ilicic et al., 2016].

Figure 8: Metabolism of RNA (R) overwritten significantly by the data (left).
The new dense factor captures the major axes of variation pertaining to
individual donors (right).

An additional dense factor is introduced during training, as a large portion
of the prior annotations of the reactome gene set metabolism of RNA was
refined significantly by the model. The new dense factor explains the
variation introduced by the different batches, revealing the corresponding
donor that provided the data as shown in Figure 8.

Figure 9: Reactome visual phototransduction factor. Top 10 absolute weights
color-coded emphasize added genes (left). Refined visual phototransduction
factor serves as a strong indicator of cone photoreceptors (right). Rod cells
have altered transparency ($\alpha=0.5$), due to their large number in the
dataset.

Among the annotated factors, we observe active pathways that are responsible
for critical processes occurring inside the cell such as the translation of
mRNA sequences for protein synthesis, or the citric acid TCA cycle, also known
as the Krebs cycle. Importantly, the model consistently identifies factors
characterized by specific processes such as visual phototransduction (R) -
i.e. the conversion of light into electrical signals in rod and cone
photoreceptors - as a highly relevant factor. Although rod photoreceptors
dominate the dataset with over 15 thousand $(75\%)$ samples, the original
factor is refined by cone specific markers such as _PDE6H_ and _ARR3_[Lukowski
et al., 2019], rendering the inferred factor also suitable for explaining the
variation within the significantly smaller cluster of $627$ cone cells as
illustrated in Figure 9.

##  4 Conclusion

We proposed spex-LVM, a flexible and interpretable factor analysis approach
that supports domain-specific interpretation of discovered latent factors. By
leveraging an efficient stochastic variational inference algorithm we
demonstrated the scalability of our model to high-dimensional data. Multiple
experimental scenarios on synthetic data showed that our model consistently
pinpoints the true axes of variation, regardless whether the prior information
was severely perturbed. As an additional benefit, spex-LVM learns identifiable
factors as the designated order of the annotated factors is preserved. Further
validation of our model on real datasets demonstrated its ability to uncover
known sources of heterogeneity, isolate unwanted technical noise, and refine
existing gene set annotations in a data driven manner.
In the following, we address some limitations of our work and potential future
extensions. First, while linear models promote interpretability, their limited
expressiveness cannot capture complex, non-linear dependencies. Unsupervised
approaches in deep learning such as the variational autoencoders (VAEs)[Kingma
and Welling, 2013] provide an attractive alternative to traditional linear
models, but risk being opaque in favor of more expressive power, and may fail
in settings where data samples are scarce.
Second, a basic understanding of the domain knowledge is essential to
determining the proper amount of bias that should be introduced by the prior
information during training. The balance between the prior belief and the
likelihood based on the data is a much more complex issue that is prevalent in
most statistical methods, and especially in approaches that integrate prior
knowledge which may be limited, or otherwise noisy.
As future work, the ideas presented in this work can be extended to
accommodate multiple views, e.g. GFA[Klami et al., 2014], MOFA[Argelaguet et
al., 2018], which would utilize the prior information of a subset of views to
indirectly provide explanatory labels for all the modalities involved in the
analysis, thereby facilitating the discovery of novel feature sets.

## References

  * Argelaguet et al. [2018] Ricard Argelaguet, Britta Velten, Damien Arnol, Sascha Dietrich, Thorsten Zenz, John C Marioni, Florian Buettner, Wolfgang Huber, and Oliver Stegle.  Multi-omics factor analysis—a framework for unsupervised integration of multi-omics data sets.  _Molecular systems biology_ , 14(6):e8124, 2018\.
  * Bishop [2006] Christopher M Bishop.  _Pattern recognition and machine learning_.  springer, 2006.
  * Buettner et al. [2015] Florian Buettner, Kedar N Natarajan, F Paolo Casale, Valentina Proserpio, Antonio Scialdone, Fabian J Theis, Sarah A Teichmann, John C Marioni, and Oliver Stegle.  Computational analysis of cell-to-cell heterogeneity in single-cell rna-sequencing data reveals hidden subpopulations of cells.  _Nature biotechnology_ , 33(2):155–160, 2015\.
  * Buettner et al. [2017] Florian Buettner, Naruemon Pratanwanich, Davis J McCarthy, John C Marioni, and Oliver Stegle.  f-sclvm: scalable and versatile factor analysis for single-cell rna-seq.  _Genome biology_ , 18(1):1–13, 2017.
  * Carvalho et al. [2009] Carlos M Carvalho, Nicholas G Polson, and James G Scott.  Handling sparsity via the horseshoe.  In _Artificial Intelligence and Statistics_ , pages 73–80. PMLR, 2009\.
  * Fabregat et al. [2018] Antonio Fabregat, Steven Jupe, Lisa Matthews, Konstantinos Sidiropoulos, Marc Gillespie, Phani Garapati, Robin Haw, Bijay Jassal, Florian Korninger, Bruce May, et al.  The reactome pathway knowledgebase.  _Nucleic acids research_ , 46(D1):D649–D655, 2018\.
  * Hicks et al. [2018] Stephanie C Hicks, F William Townes, Mingxiang Teng, and Rafael A Irizarry.  Missing data and technical variability in single-cell rna-sequencing experiments.  _Biostatistics_ , 19(4):562–578, 2018.
  * Hoffman et al. [2013] Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 14(5), 2013.
  * Ilicic et al. [2016] Tomislav Ilicic, Jong Kyoung Kim, Aleksandra A Kolodziejczyk, Frederik Otzen Bagger, Davis James McCarthy, John C Marioni, and Sarah A Teichmann.  Classification of low quality cells from single-cell rna-seq data.  _Genome biology_ , 17(1):1–15, 2016.
  * Kang et al. [2018] Hyun Min Kang, Meena Subramaniam, Sasha Targ, Michelle Nguyen, Lenka Maliskova, Elizabeth McCarthy, Eunice Wan, Simon Wong, Lauren Byrnes, Cristina M Lanata, et al.  Multiplexed droplet single-cell rna-sequencing using natural genetic variation.  _Nature biotechnology_ , 36(1):89–94, 2018.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Klami et al. [2014] Arto Klami, Seppo Virtanen, Eemeli Leppäaho, and Samuel Kaski.  Group factor analysis.  _IEEE transactions on neural networks and learning systems_ , 26(9):2136–2147, 2014.
  * Lähnemann et al. [2020] David Lähnemann, Johannes Köster, Ewa Szczurek, Davis J McCarthy, Stephanie C Hicks, Mark D Robinson, Catalina A Vallejos, Kieran R Campbell, Niko Beerenwinkel, Ahmed Mahfouz, et al.  Eleven grand challenges in single-cell data science.  _Genome biology_ , 21(1):1–35, 2020.
  * Liberzon et al. [2015] Arthur Liberzon, Chet Birger, Helga Thorvaldsdóttir, Mahmoud Ghandi, Jill P Mesirov, and Pablo Tamayo.  The molecular signatures database hallmark gene set collection.  _Cell systems_ , 1(6):417–425, 2015.
  * Lukowski et al. [2019] Samuel W Lukowski, Camden Y Lo, Alexei A Sharov, Quan Nguyen, Lyujie Fang, Sandy SC Hung, Ling Zhu, Ting Zhang, Ulrike Grünert, Tu Nguyen, et al.  A single-cell transcriptome atlas of the adult human retina.  _The EMBO journal_ , 38(18):e100811, 2019.
  * Mitchell and Beauchamp [1988] Toby J Mitchell and John J Beauchamp.  Bayesian variable selection in linear regression.  _Journal of the american statistical association_ , 83(404):1023–1032, 1988.
  * Murphy [2012] Kevin P Murphy.  _Machine learning: a probabilistic perspective_.  MIT press, 2012.
  * Piironen et al. [2017] Juho Piironen, Aki Vehtari, et al.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  _Electronic Journal of Statistics_ , 11(2):5018–5051, 2017.
  * Polson and Scott [2010] Nicholas G Polson and James G Scott.  Shrink globally, act locally: Sparse bayesian regularization and prediction.  _Bayesian statistics_ , 9(501-538):105, 2010.
  * Ranganath et al. [2014] Rajesh Ranganath, Sean Gerrish, and David Blei.  Black box variational inference.  In _Artificial intelligence and statistics_ , pages 814–822. PMLR, 2014.
  * Roche and Furuta [2015] Paul A Roche and Kazuyuki Furuta.  The ins and outs of mhc class ii-mediated antigen processing and presentation.  _Nature Reviews Immunology_ , 15(4):203–216, 2015\.
  * Rybakov et al. [2020] Sergei Rybakov, Mohammad Lotfollahi, Fabian J Theis, and F Alexander Wolf.  Learning interpretable latent autoencoder representations with annotations of feature sets.  _bioRxiv_ , 2020.
  * Shapiro [1985] Alexander Shapiro.  Identifiability of factor analysis: Some results and open problems.  _Linear Algebra and its Applications_ , 70:1–7, 1985.
  * Wingate and Weber [2013] David Wingate and Theophane Weber.  Automated variational inference in probabilistic programming.  _arXiv preprint arXiv:1301.1299_ , 2013.