{'fraction_non_alphanumeric': 0.030753473952278138, 'fraction_numerical': 0.012262636328406134, 'mean_word_length': 4.769243053337592, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Deep Generative Decoder: MAP estimation of representations improves modeling of single-cell RNA data

Learning low-dimensional representations of single-cell transcriptomics has become instrumental to its downstream analysis. The state of the art is currently represented by neural network models such as Variational Autoencoders (VAEs) which use a variational approximation of the likelihood for inference. We here present the Deep Generative Decoder (DGD), a simple generative model that computes model parameters and representations directly via maximum a posteriori (MAP) estimation. The DGD handles complex parametrized latent distributions naturally unlike VAEs which typically use overly simple fixed Gaussian distributions. We first show its general functionality and superiority in data generation on a commonly used benchmark set, Fashion-MNIST. Secondly, we apply the model to a single-cell dataset from peripheral blood mononuclear cells. Here the DGD learns low-dimensional, meaningful and well-structured latent representations with sub-clustering beyond the provided labels. The advantages of this approach are its simplicity and its capability to provide representations of much smaller dimensionality than a comparable VAE.


Introduction
High-throughput methods in biology and medicine produce vast amounts of high-dimensional noisy data that we seek to extract knowledge from. The first step is often to obtain a low-dimensional representation of the data through clustering, PCA analysis or other similar techniques. A good example is gene expression analysis, where we need to compare counts for each of the 10-20 thousand genes between samples or cells. Here the first analysis is often to visualize the data in two dimensions with UMAP [1] or t-SNE [2]. Although these tools are indispensable, they do not model the the uncertainty of the count data or reveal any of the underlying structure in the data. Therefore, it has become popular to make use of generative models that give low-dimensional representations mapping to a probability distribution over the data. Some of the most common are variational autoencoders (VAEs) [3], generative adversarial networks (GANs) [4], normalizing flows [5] and energy-based models (EBMs) [6]. Even though these approaches have their shortcomings [7], they have shown great success in a multitude of applications such as single-cell RNA-sequencing (scRNA-seq) [8,9], image and surface generation [10,11], and natural language processing [10,12].
The VAE is the most common type of generative model applied to biological data. In general, VAEs consist of an encoder that maps data to a lower-dimensional probabilistic representation in latent space and a decoder that maps back to the data. The model is trained using maximum likelihood estimation (MLE), which seeks to find a set of parameters that maximize the likelihood of the data. However, the likelihood calculation is intractable due to an integral over the latent space. The solution in the VAE is to approximate the likelihood and maximize a variational lower bound (the "ELBO"). This approximation is not always accurate [13] and favours larger latent dimensionalities than necessary [14].
The decoder alone specifies the generative model and the encoder of the VAE is just a convenient tool for finding representations, but it adds a whole new set of modelling choices to be made, which influence the inference gap and expressiveness [13]. There are simpler approaches without an encoder that maximise the likelihood over representations of the training data and the decoder parameters at the same time [15][16][17]. The representations are treated like the other model parameters and since gradients can be derived exactly, the model can be trained using standard gradient descent. The standard VAE and many other models use a fixed Gaussian distribution over latent space. This may result in under-expressive representations, and it has been shown that learning the parameters of a Gaussian mixture leads to improved generalization and clustering [18,19]. However, this makes the variational inference even more complex, as it has to be combined with score matching [20], or other approaches. Another alternative approach worth mentioning is the VQ-VAE [21] which learns discrete representations with more powerful priors than a standard Gaussian, but still requires an encoder.
In this work, we present a probabilistic formulation of a generative decoder using maximum a posteriori (MAP) estimation of all parameters. MAP estimation is the Bayesian analog of MLE and seeks to find the set of parameters that maximize the posterior, which is the probability of the representations and model parameters given the observed data. There is no intractable step in this estimation, so both model parameters and representations can be estimated directly. The model we are introducing consists of a decoder (or generator), such as a feed-forward neural network, and a distribution over representations with learnable parameters, which are optional. Estimating the representations of the training samples as well as the parameters is straight-forward and can be done by gradient descent as in a non-probabilistic formulation [22]. One of the main advantages of this approach over variational inference is the ease and flexibility with which distributions over representations can be learned. Our emphasis here is on a model with a parameterized distribution over representations. This is motivated by the intuition that representations are meaningful and are likely to group samples with different properties. We think of this approach as a marriage of generative models and manifold learning, such as UMAP. The beauty of the approach lies in its simplicity and utility. The decoder parameters, latent representations and latent distribution are estimated in the same way. The simplicity makes it easy to extend the model to more complex distributions and losses. It also makes it much simpler to understand for non-specialists. This is especially important, as we see great potential for this approach in fields like biology and medicine.
In the following sections we present the theory of the model and experiments on two very different types of data. We first demonstrate the general functionality of the proposed model on the Fashion-MNIST [23] benchmark dataset. This is a typical choice in deep learning and an important step, because it enables us to investigate and understand the model's generative capabilities and potential caveats. We further explore the complexity of the distribution over representations and compare latent space and generative capabilities to models using variational inference. We also showcase the extension of the approach to supervized learning in which the mixtures of the latent distribution represent specific data classes. After establishing the model's functionality, we show an application to single-cell gene expression count modelling. We chose a data set of peripheral blood mononuclear cells (PBMC) [24], for which we compare our approach in terms of data reconstruction and clustering performance to scVAE [25], which is a variational autoencoder also using an adaptive mixture of Gaussians rather than a single fixed Gaussian to model the latent representation. We show that our model learns a representation and Gaussian mixture model which cluster well according to cell type with additional, previously unobserved sub-clustering. This is achieved at a much lower dimensionality than scVAE. This is an important advantage, as lower-dimensional spaces are much more interpretable and provide meaningful spacial relationships. Fig. 1 Schematic of the DGD. The deep generative decoder consists of a decoder of any desired architecture and parameters θ mapping the latent representation Z to the data space X. The representation is modelled by a parametrized probabilistic distribution, here represented as a Gaussian mixture model (GMM) with K components and priors φ.

The model
Consider a model with observed variable x and continuous latent variable z, which we also refer to as a representation. Usually, the x-space (or sample space) is of a higher dimension than the z-space, so the model gives a low-dimensional representation of data. A neural network, the decoder, with parameters θ takes z as input and outputs f θ (z), which are the parameters for the distribution of x. These could be the mean and standard deviations of independent normal distributions. Thus, the neural network defines a conditional distribution in sample space, P (x | z, θ) = P (x | f θ (z)). We assume a distribution over representation space, P (z | φ), with parameters φ. The parameters may be fixed as in a standard VAE with a fixed Gaussian distribution over representations. It can also have adjustable parameters, such as a mixture of Gaussians with trainable means, covariances and mixture coefficients, but essentially any (differentiable) parametrized distribution can be used. When introducing priors over parameters, the joint probability of everything becomes Again, P (x | z, θ) represents the decoder neural network, P (z | φ) the distribution over representations, P (θ) the prior over the decoder parameters (neural network weights), and P (φ) the prior over the parameters in the representation space distribution. We use maximum a posteriori (MAP) estimation to find optimal parameters and representations. For a data set X = {x 1 , x 2 , . . . , x N } with N samples we thus want to find representations Z = {z 1 , z 2 , . . . , z N } and parameters φ, θ that maximize P (Z, φ, θ | X): argmax Z,φ,θ P (Z, φ, θ | X) = argmax Z,φ,θ P (X, Z, φ, θ)/P (X) = argmax Z,φ,θ P (X, Z, φ, θ).
(2) If we assume independence among training samples, the log of this joint probability can be obtained from (2), with i indexing the samples log P (X, Z, φ, θ) = i log P (x i | z i , θ) + log P (z i | φ) + log P (φ) + log p(θ).
(3) This is the quantity we want to maximise with respect to both the representations (z) for all samples and the model parameters (φ, θ). The training process can be imagined like this: A representation z for each training sample is initialized. This can be a random sample from a chosen distribution or a zero-valued vector. Representations are then passed through the decoder to give P (x | z, θ) and losses are computed from this and P (z | φ). All parameters θ, φ, and z are then updated via back-propagation. The implementation is efficient and straightforward and compatible with standard modules and loss functions in deep learning frameworks, such as Pytorch. Once the model is estimated, prediction on new data points is done by first finding an optimal representation by maximizing P (x | z, θ)P (z | φ) as above while keeping all other model parameters fixed, i.e., do gradient descent in z alone.
In this work we use a parametrized Gaussian mixture for the representations (P (z | φ)) with diagonal covariance matrix and a mollified Uniform ("softball prior") as prior over component means. We use a Dirichlet prior as a prior on the mixture weights and a Gaussian prior on the log-variances of the components. We use a weight decay in the training of the decoder, which corresponds to a Gaussian prior on the weights.

Demonstration on image data
It is customary to test generative models on image benchmarks, so we begin our experimental analysis by applying our model to the Fashion-MNIST data. It allows us to understand and evaluate model behaviour, latent space and sampling quality much better than more complex, biological data. The decoder used for the Fashion-MNIST models is based on the DCGAN [26] generator architecture with modifications using common methods in image generation and is described in detail in the Methods section and figure 7. The best results were achieved with a latent dimension of 20 and a GMM with 20 components. The remaining parameters are described in the Methods section. The model was trained for 500 epochs.
In figure 2b) we see that the 20 Gaussian components distribute well over the complex latent space. The resulting distribution over latent space (the parametrized GMM) models the latent representation much better than for GMMs with less components which fail to pick up on more subtle structures (Supplementary figure 5). This supports the hypothesis stated in many works that overly simple distributions such as standard Gaussians may not be sufficient for describing the probabilistic representations of many data. Another advantage of the GMM is the increase in control over sampling. Figure 2a) shows generated samples from individual components of the GMM. Having more components leads to each of them covering a much smaller area of the latent space, and thus offering less varied but more specific samples. In this case of a 20-dimensional latent space with 20 mixture components, each component mean (figure 2b) top) can be attributed to a (sub-)class of Fashion-MNIST. A relatively high number of components can thus enable more deliberate, more qualitative sampling and a better model of the latent space. While the latent representations in figure 2b) show clear separation for super-classes such as shoes (Ankle boot, Sneaker, Sandal) and some distinct classes like of trousers and bags, other classes are much more mixed. We evaluate the clustering with the adjusted Rand index (ARI) [27]. This metric achieves values between 0 and 1 with 1 representing a perfect clustering and is corrected for chance. The ARI of this clustering based on predicted labels derived from the GMM's component probabilities is only 0.295. As an extension of the DGD, it is also possible to perform supervized training, where each GMM component is assigned a data class a priori. In the optimization, each component is only committed to covering the latent space of its assigned classes' representations. When trained with these assigned 10 components, we arrive at the distinct latent clusters seen in figure 2c) with an ARI of 0.995. While the quality of random samples from this model is reduced slightly (see Supplementary figure 6) in comparison to those from an equivalent unsupervized model, we achieve perfect control over the sample identity.
We next investigated how the MAP approach compares to Variational Inference (VI) and amortization. We found that once the encoder is removed (using the Variational Auto-Decoder (VAD) [28]), reconstruction performance becomes much better than that of the standard VAE and similar to that of the DGD, as seen in 2d) and Supplementary table 3. Supplementary table 3 also shows that latent space normality is improved by removing the encoder. The DGD, however, achieved higher-quality generated images compared to the VAD as shown in figure 2e), and even works best for low-dimensional representations.

Application to single-cell data
In order to assess the model performance on scRNA data, we chose to apply it to the peripheral blood mononuclear cells (PBMC) data from [24]. The single  Table 1 Comparison of performance metrics for scDGD, scVAE and scVI. Model names are indicated on the left, with c describing the number of Gaussians in the GMM (1 indicates a standard Gaussian) and d the number of latent dimensions. The asterisk next to scVAE indicates that its values are those reported in [25]. Performance metrics computed are indicated as the remaining columns. The NLL refers to the negative log-likelihood for scDGD and the NLL lower bound for scVAE and scVI, indicated by LB. RMSE and MAE are missing for scVAE. The best values for each metric are highlighted in bold.
cell expression counts are modeled by a Negative Binomial distribution. We achieved best clustering results with a latent dimensionality of 20 and 3 fully connected hidden layers. The model, which we will refer to as scDGD, was trained for 516 epochs after which it achieved its overall best performance in terms of reconstruction and clustering. The number of GMM components was set to 9, which represents the number of cell types in the data. Details about our implementation can be found in the Methods section.
The resulting latent space is depicted as a t-SNE projection in figure 3. We found that 2D PCA is insufficient for showing fine-grained clusters in this space, and UMAP separated larger clusters so much that fine-grained subclustering is again not visible (see Supplementary figure 8). From a visual perspective, the latent space and GMM components cluster well with the cell types. For clustering purposes, a cell is assigned to the component that gives the largest probability to its representation. The adjusted Rand index (ARI) of this clustering on the train set is 0.585. This is not directly comparable to scVAE, which reports a higher value of 0.656 when training on all data. Because negative log-likelihood (NLL) of the scDGD and the NLL lower bound (or "ELBO") of VAE-based models are not equivalent, we chose to compare the reconstruction performance as the root mean squared error (RMSE) and mean absolute error (MAE) of the held-out test set predictions. We computed these metrics for scDGD and scVI [8], which performed better in terms of reconstruction but worse in terms of latent space clustering according to [25]. Training and evaluating scVI on our data split reveals higher RMSE and MAE on the test set and a much lower ARI of 0.470 compared to scDGD (see Table  1).
We also trained a model with 18 Gaussian components, since we have learned in our explorations on image data that increasing the complexity of the distribution over representation can be beneficial to modeling the substructures in the representation if the problem is complex enough. Increasing the number of components resulted in equal reconstruction performance as seen in Table 1 and highly improved clustering performance with an ARI of 0.649. This makes the clustering performance of our approach equivalent to scVAE as shown in Table 1 and may imply better reconstruction performance as we are outperforming scVI which was superior to scVAE in this regard [25]. This is achieved with a model using only 78% of the parameters (including representations of the full dataset). As we can see in figure 3c), several cell types have one or more components specifically assigned to them. Among these are CD4+ naïve and memory T cells, NK, CD34+ and B cells. The cytotoxic T cell subtypes cannot easily be differentiated by component assignment, but are distinct from all other cell types in the data. The same goes for all CD4+ T cells except memory cells. Since t-SNE distorts the latent space, we also show a graph of the component means in Figure 3d) based on the Euclidean distances between them (Supplementary Figure 9). To go one step further in the analysis of our model, we next aimed to get an understanding of the structuring of the learned latent space. In both models with 9 and 18 components, we see that the CD34-positive cells form subclusters that are modelled by different Gaussian components. We therefore investigated the representations for CD34-positive cells learned by the 18-component scDGD in terms of some more fine-grained differentiation markers. We were able to determine a subcluster best modeled by component 17 (see 4a)) which is primarily occupied by cells expressing markers for the lymphoid lineage of hematopoetic stem cells (HSCs). This cluster further shows a distinction between bone marrow common lymphoid progenitors (CLPs) and cord blood CLPs, indicated by markers CD10 and CD7, respectively. Another subcluster, represented by components 8 and 16, is made up of several cells expressing CD18, which is a marker for myeloid cells regulating neutrophil production and a general marker for granulocytes. There is also a strong presence of CD14-positive cells in these clusters, suggesting an aggregation of monocytes and neutrophils. The biggest subcluster, modeled by component 14, suggests the presence of different subpopulations of the myeloid lineage of HSCs. We find gradual levels of CD41 to CD62L, which are markers for megakaryoblasts and myeloblasts, respectively. These are also correlated with component 7. We additionally found a concentrated site showing expression of different dendritic cell markers (CD1c, CD141, CD303), best represented by components 8 and 16.

Discussion
In this work, we have presented the DGD, a Bayesian formulation of a deep generative neural network using MAP estimation for model parameters and representations. This is a fully tractable and simple approach, for which no encoder is needed. This gives the DGD the advantage of having fewer parameters and fewer modelling choices, such as encoder architecture. We have shown that the DGD achieves good reconstruction, clustering and generative performance on much smaller latent spaces than some of the other generative approaches such as VAEs. In this sense it is comparable to flow-based models [29], but with less architectural restrictions. Additionally, we see clear advantages of this approach due to the simplicity of incorporating more complex distributions for modelling the latent representation. Analogous to our results on Fashion-MNIST, others have reported much better clustering and generative performances with a mixture of Gaussians compared to classic VAEs with standard Gaussians as priors over latent space [25].
In our application to single-cell expression counts, we showed that we can achieve reconstruction and clustering performances equivalent to those of comparable models, but with much smaller and better structured latent spaces, fewer parameters and more detailed distributions. Our choice for comparison is the scVAE, a well-performing VAE with a GMM modelling the representation of single-cell expression data [25]. We compared scVI and the reported losses of scVAE to our models with 9 and 18 components. We find that our model with 18 Gaussian components performs as well as scVAE in clustering and potentially better in reconstruction. The scDGD may, however, have learned a superior representation with much lower dimensionality in addition to the number of parameters being nearly half of those of scVAE. When compared to scVI, our model even performs better on both reconstruction and clustering. We find that the structure of the learned latent space shows meaningful relationships and sub-clustering by Gaussian components, which can be linked to different subpopulations in some cell types. This type of analysis was not performed for scVAE and not possible with scVI.
A caveat of this model is the inference time of new data points. The lack of an encoder means that the representations of new data points need to be found by back-propagation (with fixed model parameters). This process is more time-consuming and computationally expensive than passing new data points through an encoder as in the VAE-based model. However, in many applications such as single-cell genomics, computation time is not a serious issue. If it is an issue, one can train an encoder afterwards as described in [30], which still results in a more expressive representation than the traditional AE setup.
Altogether, we find that the DGD is a simple and efficient approach to learning low-dimensional representations. Unlike some other generative models, it is fully tractable and bears no architectural restrictions. Although the model can be supervised as we showed on Fashion-MNIST, our scRNA model is completely unsupervized and can thus lend itself to data without existing labels, while achieving a meaningful clustering of the representation on its own. We also observed that a more complex distribution relating to more GMM components can be benefitial for discovering more substructures in the data than provided by class labels. This makes the DGD an ideal candidate for modelling a multitude of biological and medical data. For these types of data, interpretability of the latent space is crucial. Several downstream applications such as perturbations and pseudo-time require a low-dimensional and continuous representation. We plan to expand the model to a variety of biological tasks and to include semi-supervized learning and other features that are of interest to generative models for biological and medical data. We further aim to apply this approach to more extensive single-cell expression datasets and multi-omics data that allow a more thorough investigation of the latent space and can provide meaningful tools for data integration and downstream analysis.

Datasets
In this work we made use of two publicly available datasets. The first is the natural image Fashion-MNIST [23] set. We used the data set's implemented train-test split. The second is a single-cell gene expression count data set of peripheral blood mononuclear cells (PBMC) presented by [24]. The data is provided by 10x Genomics under "Single Cell 3' Paper: Zheng et al. 2017 (v1 Chemistry)" and consists of data from the following 9 cell types: CD4+/CD45RA+/CD25-naïve T cells, CD4+ helper T cells, CD4+/CD25+ regulatory T cells, CD4+/CD45RO+ memory T cells, CD8+/CD45RA+ naïve cytotoxic T cells, CD8+ cytotoxic T cells, CD56+ natural killer cells, CD34+ cells, and CD19+ B cells. As [25] we used the filtered gene-cell matrices. This data is extremely sparse with 98% of the data being zero [25]. Of the 32738 genes covered by this assay, only 21812 are expressed in the whole dataset at least once. The 92043 samples are randomly split into train, validation and test set using percentages 81-9-10%. The validation set is used to finding hyperparameters, and the test set is used for final evaluation.

The DGD
The DGD consists of three components. Representation Z is learned directly as a set of trainable parameters. For each of the N samples of X, there exists an m-dimensional vector z that receives gradients. The values are initialized with zero for training but can be initialized with any other tensor of shape (N, m). The decoder is characterized by the input dimension m and output dimension n, which represents the number of data features. The decoder can be any type of Neural Network that fulfils this constraint. The choice of reconstruction loss is equally flexible but should ideally represent P (x | z). The distribution over latent space, which makes this a generative model, is here presented by a GMM. The mixture model consists of K mixture components each with an m-dimensional mean vector µ µ µ, and diagonal covariance Σ and a mixture coefficient c for each component. The mixture coefficients are transformed into component weights w through softmax activation. Since the diagonal covariance matrix can only take positive values, the according parameter is learned as the negative log-covariance.
The priors over the GMM parameters φ are as follows. The prior on the mixture weights is an m-dimensional Dirichlet distribution with uniform parameters, α. The prior on mixture means is a mollified uniform which we call the "softball" distribution with hyperparameters scale as the radius of the m-ball and sharpness determining the slope of the boundary. Details on initialization and formulation of the log-probability can be found in Supplementary Section 4.5. For the negative log-variance, we use a Gaussian prior N (−2 log(σ), 1) with the same mean and standard deviation for all dimensions. Mixture coefficients and negative log-variances are initialized with default values 1 (so that components weights are uniformly K −1 ) and −2 log(σ) with σ = 0.2 × scale × K −1 , respectively. The default values for α, scale, and sharpness are all 1. Component means are initialized by sampling from the softball prior. The distributional component of the loss is given as with P (z i | µ µ µ k , c k , −log Σ k ))) as the density of the kth multivariate Gaussian component for representation z i . The training procedure is described in Algorithm 1. Since representations are updated once per epoch and decoder and GMM parameters every batch pass, we use different instances of the optimizer for the different parameter sets. We have also observed that the GMM requires much larger steps than neural network parameters, so we choose to have an optimizer per parameter set with individually selected learning rates.  little longer, usually between 10-20 epochs. Learning of new representations is of course done with frozen decoder and GMM parameters.

Fashion-MNIST DGD
For the image-generating DGD, we choose a network architecture based on convolutions and tricks from other image generation models and image segmentation techniques [31][32][33]. The decoders start with two fully connected hidden layers fed with the latent representations. The first hidden layer has 100 units, the second capacity × 3 × 3 units. The capacity represents the minimum number of channels of the convolutional layers except for the grayscale output channel. The output of the fully connected hidden layers is reshaped into (batch, 3, 3, capacity) and fed into a series of NN blocks. These blocks consist of a transposed convolutional layer, swish activation and a Squeeze and Excitation layer as applied by [33]. Input-and output channels as well as kernel size, strike and padding depend on the position of the block in the scheme and can be seen in Supplementary Figure 7. There are four main blocks and two skip connection blocks. The outputs of combined blocks go through the activation function after summation. The series of these NN blocks is followed by a PixelCNN with 5 layers and mask size 5 [32] and a last simple transposed convolutional layer reducing the number of channels to 1. The output is scaled using sigmoid activation.
Latent dimensionality and convolutional capacity vary depending on the experiment. For the models with 10 and 20 Gaussian components, we use a latent dimensionality of 20 and a capacity of 64. For the models involving a standard Gaussian (including VAD and VAE), we test latent dimensionalities of 20, 50, and 100 with a capacity of 64. The supplementaries contain results of a model with latent dimensionality 2, capacity 32 and varying number of Gaussian components (1, 10, and 20). The number of parameters of the DGD models thus vary between roughly 1.5 and 3.5 million parameters.
For training, we use the binary cross entropy (BCE) loss and Adam optimization with betas 0.5 and 0.7. For decoder, representation and GMM we choose learning rates 1e − 3, 1e − 2 and 1e − 1, respectively. The GMM hyperparameters for different models can be found in Table 2.

VAD and VAE
The VAD and VAE decoder implementations tested here are architecturally identical to the DGDs with a standard Gaussian. The VAE encoder was for simplicity chosen to mirror the decoder, with normal convolutional layers instead of the transposed convolutions. The capacity for all models is 32 and latent dimensions tested are 20, 50, and 100. Instead of the negative log density of the GMM, the loss for the distribution over latent space is given as the Kullback-Leibler divergence as typical for VI [3].

Supervized learning
In the unsupervized model, the log likelihood log P (z i | φ) of an individual representation z i being drawn from the parametrized distribution P (φ) is given by the log sum of the probability densities of z i per component times the corresponding component probabilities, which are given by the softmax of the weights. In the supervized setting, we can ignore all components except the one that has been assigned to the sample's class. We thus only have to calculate the probability density of z i for the given component multiplied with the component probability. The losses for all other components will be zero and they will thus not receive gradients.

Single-cell DGD
The architecture of scDGD is much simpler than that of the Fashion-MNIST DGD. The decoder consists of the input layer with units defined by the latent dimensionality, and 3 hidden layers of each 100 units, connected to the output layer of 32728 units (representing all transcripts in the data). The number of output units is given by the number of genes with non-zero expression counts in the data. The latent dimensionality was empirically found best to be 20.
The normalized expression counts are modelled with a Negative Binomial distribution. Normalization is achieved by dividing the true expression count with the samples largest count. We can therefore use Sigmoid activation on the output layer. The reconstruction loss is given as the probability density of the rescaled expression count and a gene-specific, learned dispersion parameter. Since this parameter is positive definite, we learn its logarithmic counterpart. For training, we use Adam optimization with betas 0.5 and 0.7. For decoder, representation and GMM we choose learning rates 1e − 3, 1e − 2 and 1e − 2, respectively. Two models are trained with 9 and 18 Gaussian components, respectively.
The component means of the GMM are distributed according to the mollified Uniform, for which we set scale and sharpness to 1 each. This applies to both 9-component and 18-component model. The standard deviation of the components is initialized with a mean of 0.02 for the 9-component model and 0.01 for the 18-component model. Dirichlet alphas are set to 2 and 1 for models with 9 and 18 components, respectively.

scVI
An scVI model was trained on our PBMC train set. The model is implemented in Python version 3.8 with the scvi-tools [34] package version 0.17.4. as described in [8] with one hidden layer of 128 hidden units in both encoder and decoder. We chose a latent dimensionality of 20 for comparability with our model. The dropout rate is set to 0.1. The model was trained for 107 epochs. For the analysis and comparison to scDGD, the latent space is clustered with Kmeans clustering and then evaluated by the ARI. The lower bound is computed on our heldout test set. Latent representations and ELBO were computed using scVI's implemented functions.

Clustering metric
We use the adjusted Rand index (ARI) [27] for evaluating and comparing clustering performance of our models. The value of the metric ranges from 0 to 1, representing random and identical clustering, respectively. This metric is relatively robust to different clustering approaches and diverging numbers of clusters, which makes it a good metric for comparing our method to others.

Reconstruction metrics
For Fashion-MNIST, we use the binary cross-entropy (BCE) loss in order to evaluate the reconstruction of image pixels.
In scDGDs, we calculate the negative log likelihood as the summed logdensity of parametrized Negative Binomials over all genes. For each gene, the log density is calculated based on the rescaled model output and a genespecific, learned, dispersion factor.

Software, statistics and visualization methods
Python 3.8 is used as the programming basis for all methods. Neural Networks are implemented in Pytorch 1.10 [35]. Dimensionality reduction is either done with our described model or common methods such as PCA, t-SNE [2] or UMAP [1]. Graph visualizations are achieved using the NetworkX package [36]. Requirements can be found in the corresponding repository here.

The softball prior
In the softball prior used on the GMM means µ µ µ, the means are initialized as µ µ µ = scale × l × u u with u ∼ N (0, 1) and l ∼ U(0, 1) m −1 .
The log-probability of this prior is given as log P (µ µ µ) = A − log 1 + e sharpness( µ µ µ scale −1) with A = log Γ(1 + 0.5m) − m(log (scale) + 0.5 log π) The normalization constant can be understood as the volume of the m-ball and represents an approximation of the true normalization constant.    Table 3 Metrics of the quantitative analysis of DGD, VAD and VAE trained and evaluated on Fashion-MNIST. The BCE loss represents the averaged Binary Cross-Entropy loss of the reconstructed images. The lower it is, the better. The Shapiro-Wilk test statistic describes the Normality of the latent representation. The test statistic lies within the range of 0 to 1, the higher the better.