G Durif

L Modolo

J E Mold

S Lambert-Lacroix

F Picard

Probabilistic Count Matrix Factorization for Single Cell Expression Data Analysis

Motivation:

The development of high throughput single-cell sequencing technologies now allows the investigation of the population diversity of cellular transcriptomes. The expression dynamics (gene-to-gene variability) can be quantified more accurately, thanks to the measurement of lowly-expressed genes. In addition, the cell-to-cell variability is high, with a low proportion of cells expressing the same genes at the same time/level. Those emerging patterns appear to be very challenging from the statistical point of view, especially to represent a summarized view of single-cell expression data. PCA is a most powerful tool for high dimensional data representation, by searching for latent directions catching the most variability in the data. Unfortunately, classical PCA is based on Euclidean distance and projections that poorly work in presence of over-dispersed count data with dropout events like single-cell expression data. Results: We propose a probabilistic Count Matrix Factorization (pCMF) approach for single-cell expression data analysis, that relies on a sparse Gamma-Poisson factor model. This hierarchical model is inferred using a variational EM algorithm. It is able to jointly build a low dimensional representation of cells and genes. We show how this probabilistic framework induces a geometry that is suitable for single-cell data visualization, and produces a compression of the data that is very powerful for clustering purposes. Our method is competed against other standard representation methods like t-SNE, and we illustrate its performance for the representation of single-cell expression (scRNA-seq) data.

Introduction

The combination of massive parallel sequencing with high-throughput cell biology technologies has given rise to the field of single-cell genomics, which refers to techniques that now provide genome-wide measurements of a cell's molecular profile either based on DNA [START_REF] Zong | Genome-wide detection of single-nucleotide and copy-number variations of a single human cell[END_REF], RNA [START_REF] Picelli | Smart-seq2 for sensitive full-length transcriptome profiling in single cells[END_REF], or chromatin [START_REF] Buenrostro | Single-cell chromatin accessibility reveals principles of regulatory variation[END_REF][START_REF] Rotem | Single-cell ChIP-seq reveals cell subpopulations defined by chromatin state[END_REF]. Similar to the paradigm shift of the 90s characterized by the first molecular profiles of tissues [START_REF] Golub | Molecular classification of cancer: class discovery and class prediction by gene expression monitoring[END_REF], it is now possible to characterize molecular heterogeneities at the cellular level [START_REF] Deng | Single-cell RNA-seq reveals dynamic, random monoallelic gene expression in mammalian cells[END_REF][START_REF] Saliba | Single-cell RNAseq: Advances and future challenges[END_REF]. A tissue is now viewed as a population of cells of different types, and many fields have now identified intra-tissue heterogeneities, in T cells [START_REF] Buettner | Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells[END_REF], lung cells [START_REF] Trapnell | The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells[END_REF], or intestine cells [START_REF] Grün | Single-cell messenger RNA sequencing reveals rare intestinal cell types[END_REF]. The construction of a comprehensive atlas of human cell types is now within our reach [START_REF] Wagner | Revealing the vectors of cellular identity with single-cell genomics[END_REF]. The characterization of heterogeneities in single-cell expression data thus requires an appropriate statistical model, as the transcripts abundance is quantified for each cell using read counts. Hence, standard methods based on Gaussian assumptions are likely to fail to catch the biological variability of lowly expressed genes, and Poisson or Negative Binomial distributions constitute an appropriate framework [START_REF] Chen | Detection of high variability in gene expression from single-cell RNA-seq profiling[END_REF]Riggs and Lalonde, 2017, Chap. 6). Moreover, dropouts, either technical (due to sampling difficulties) or biological (no expression or stochastic transcriptional activity), constitute another major source of variability in scRNA-seq (single-cell RNAseq) data, which has motivated the development of the so-called Zero-Inflated models [START_REF] Kharchenko | Bayesian approach to single-cell differential expression analysis[END_REF].

Principal component analysis (PCA) is one of the most widely used dimension reduction technique, as it allows the quantification and visualization of variability in massive datasets. It consists in approximating the observation matrix X [n×m] (n cells, m genes), by a factorized matrix of reduced rank, denoted UV T where U [n×K] and V [m×K] represent the latent structure in the observation and variable spaces respectively. This projection onto a lower-dimensional space (of dimension K) allows one to catch gene co-expression patterns and clusters of individuals. PCA can be viewed either geometrically or through the light of a statistical model [START_REF] Landgraf | Generalized principal component analysis: Projection of saturated model parameters[END_REF]. Standard PCA is based on the 2 distance as a metric and is implicitly based on a Gaussian distribution [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF]. Model-based PCA offers the unique advantage to be adapted to the data distribution. It consists in specifying the distribution of the data X [n×m] through a statistical model, and to factorize E(X) instead of X. In this context the 2 metric is replaced by the Bregman divergence which is adapted to maximum likelihood inference [START_REF] Collins | A generalization of principal components analysis to the exponential family[END_REF]. A probabilistic version of the Gaussian PCA was proposed by [START_REF] Pierson | ZIFA: Dimensionality reduction for zero-inflated singlecell gene expression analysis[END_REF] in the context of single cell data analysis, with the modeling of zero inflation (the ZIFA method). ScRNA-seq data may be better analyzed by methods dedicated to count data such as the Non-negative Matrix Factorization (NMF) introduced in a Poisson-based framework by [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] or the Gamma-Poisson factor model [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF][START_REF] Févotte | Nonnegative matrix factorizations as probabilistic inference in composite models[END_REF][START_REF] Landgraf | Generalized principal component analysis: Projection of saturated model parameters[END_REF]. None of the currently available dimension reduction methods fully model single-cell expression data, characterized by over-dispersed zero inflated counts [START_REF] Kharchenko | Bayesian approach to single-cell differential expression analysis[END_REF][START_REF] Zappia | Splatter: Simulation of single-cell RNA sequencing data[END_REF].

Our method is based on a probabilistic count matrix factorization (pCMF). We propose a dimension reduction method that is dedicated to over-dispersed counts with dropouts, in high dimension. In particular, gene expression can be normalized but does not require to be transformed (log, Anscombe) in our framework. Our factor model takes advantage of the Poisson Gamma representation to model counts from scRNA-seq data [START_REF] Zappia | Splatter: Simulation of single-cell RNA sequencing data[END_REF]. In particular, we use Gamma priors on the distribution of principal components. We model dropouts with a Zero-Inflated Poisson distribution [START_REF] Simchowitz | Zero-Inflated Poisson Factorization for Recommendation Systems[END_REF], and we introduce sparsity in the model thanks to a spike-and-slab approach (Malsiner-Walli and [START_REF] Lun | SingleCellExperiment: S4 Classes for Single Cell Data. R package version 1[END_REF] that is based on a two component sparsity-inducing prior on loadings [START_REF] Titsias | Spike and slab variational inference for multi-task and multiple kernel learning[END_REF]. We propose a heuristic to initialize the sparsity layer based on the variance of the recorded variables, acting as an integrated gene filtering step, which is an important issue in scRNA-seq data analysis [START_REF] Soneson | Bias, robustness and scalability in single-cell differential expression analysis[END_REF]. The model is inferred using a variational EM algorithm that scales favorably to data dimension compared with Markov Chain Monte-Carlo (MCMC) methods [START_REF] Hoffman | Stochastic Variational Inference[END_REF][START_REF] Blei | Variational inference: A review for statisticians[END_REF]. Then we propose a new criterion to assess the quality of fit of the model to the data, as a percentage of explained deviance, following a strategy that is standard in the Generalized Linear Models framework. Moreover, we show that our criterion corresponds to the percentage of explained variance in the PCA case, which makes it suitable to compare geometric and probabilistic methods.

We show the performance of pCMF on simulated and experimental datasets, in terms of visualization and quality of fit. Moreover, we show the benefits of using pCMF as a preliminary dimension reduction step before clustering or before the popular t-SNE approach (van der [START_REF] Van Der Maaten | Visualizing Data using t-SNE[END_REF][START_REF] Amir | viSNE enables visualization of high dimensional single-cell data and reveals phenotypic heterogeneity of leukemia[END_REF]. Experimental published data are used to show the capacity of pCMF to provide a better representation of the heterogeneities within scRNA-Seq datasets, which appears to be extremely helpful to characterize cell types. Finally, our approach also provides a lower space representation for genes (and not only for cells), contrary to t-SNE. pCMF is implemented in the form of a R package available at https://github.com/gdurif/pCMF.

Count Matrix Factorization for zero-inflated overdispersed data

The Poisson factor model. Our data consist of a matrix of counts (potentially normalized but not transformed), denoted by X ∈ N n×m , that we want to decompose onto a subspace of dimension K (being fixed). In a first step we suppose that the data follow a multivariate Poisson distribution of intensity Λ. Following the standard Poisson Non-Negative Matrix Factorization (Poisson NMF, [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], we approximate this intensity such that

X ∼ P Λ , Λ UV T , (1)
with factor U ∈ R +,n×K the coordinates of the n observations (cells) in the subspace of dimension K, and loadings V ∈ R +,m×K the contributions of the m variables (genes).

Modeling over-dispersion. We account for over-dispersion by using the the Negative-Binomial distribution [START_REF] Anders | Differential expression analysis for sequence count data[END_REF], through a hierarchical Gamma-Poisson representation (GaP) [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF]. In our factor model U and V are modeled as independent random latent variables with Gamma distributions such that

U ik ∼ Γ(α k,1 , α k,2 ) for any (i, k) ∈ [1 : n] × [1 : K] , V jk ∼ Γ(β k,1 , β k,2 ) for any (j, k) ∈ [1 : m] × [1 : K] . (2)
In practice, some third-party latent variables are introduced for the derivation of our inference algorithm [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF][START_REF] Zhou | Beta-negative binomial process and Poisson factor analysis[END_REF]. We consider latent vari-

ables Z = [Z ijk ] ∈ R n×m×K , defined such that X ij = k Z ijk .
These new indicator variables quantify the contribution of factor k to the data. Here Z ijk are assumed to be conditionally independent and to follow a conditional Poisson distribution, i.e. Z ijk | U ik , V jk ∼ P(U ik V jk ). Thus, the conditional distribution of X ij remains P( k U ik V jk ) thanks to the additive property of the Poisson distribution.

Dropout modeling with a zero-inflated (ZI) model. To model zero-inflation, i.e. random null observations called dropout events, we introduce a dropout indicator variable D ij ∈ {0, 1} for i = 1, . . . , n and j = 1, . . . , p (c.f. [START_REF] Simchowitz | Zero-Inflated Poisson Factorization for Recommendation Systems[END_REF]. In this context, each D ij = 0 if gene j has been subject to a dropout event in cell i, with D ij ∼ B(π D j ). We consider gene-specific dropout rates, π D j , following recommendations of the literature [START_REF] Pierson | ZIFA: Dimensionality reduction for zero-inflated singlecell gene expression analysis[END_REF]. Thus, to include zero-inflation in the probabilistic factor model, we consider that

X ij | U i , V j , D ∼ (1 -D ij ) × δ 0 + D ij × P k U ik V jk ,
where δ 0 is the Dirac mass at 0, i.e. δ 0 (X ij ) = 1 if X ij = 0 and 0 otherwise. The dropout indicators D ij are assumed to be independent from the factors U ik and V jk . Then, by integrating D ij out, the probability of observing a zero in the data illustrates the two potential sources of zeros and becomes

P X ij = 0 | U i , V j ; π = (1 -π D j ) + π D j exp -k U ik V jk .
Thus, inference will be based on the factors U ik and V jk and on probabilities π D j .

Probabilistic variable selection. Finally we suppose that our model is parsimonious. We consider that among all recorded variables, only a proportion carries the signal and the others are noise. To do so, we modify the prior of the loadings variables V jk , to consider a sparse model with a two-group sparsity-inducing prior [START_REF] Engelhardt | Bayesian Structured Sparsity from Gaussian Fields[END_REF]. The model is then enriched by the introduction of a new indicator variable S jk ∼ B(π S j ), that equals 1 if gene j contributes to loading V jk , and zero otherwise. π S j stands for the prior probability for gene j to contribute to any loading. To define the sparse GaP factor model, we modify the distribution of the loadings latent factor V jk , such that

V jk |S jk ∼ (1 -S jk ) × δ 0 + S jk × Γ(β k,1 , β k,2 ) .
This spike-and-slab formulation [START_REF] Mitchell | Bayesian variable selection in linear regression[END_REF] ensures that V jk is either null (gene j does not contribute to factor k), or drawn from the Gamma distribution (when gene j contributes to the factor). The contribution of gene j to the component k is accounted for in the conditional Poisson distribution of X ij , with

X ij |U i , V j , D, S j ∼ (1 -D ij )(1 -S jk ) × δ 0 + P D ij k U ik [S jk V jk ] ,
where

V jk = S jk V jk such that V jk ∼ Γ(β k,1 , β k,2 ).
Underlying geometry. Knowing U and V, to quantify the approximation of matrix X by UV T , we consider the Bregman divergence, that can be viewed as a generalization of the Euclidean metric to the exponential family (see [START_REF] Collins | A generalization of principal components analysis to the exponential family[END_REF][START_REF] Banerjee | Clustering with Bregman Divergences[END_REF][START_REF] Chen | Metrics defined by Bregman Divergences[END_REF]. In the Poisson model, the Bregman divergence between X and UV T is defined as [START_REF] Févotte | Nonnegative matrix factorizations as probabilistic inference in composite models[END_REF]:

D(X | UV T ) = n i=1 m j=1 x ij log x ij k U ik V jk -x ij + k U ik V jk .
Hence the geometry is induced by an appropriate probabilistic model dedicated to count data. Potential identifiability issues are addressed in Supp. Mat. (Section S.2).

In the following, we will refer to pCMF for the method implementing the model with dropout but the without sparsity layer, and to sparse pCMF (or spCMF) for the model with dropout and sparsity layers.

2.1 Quality of the reconstruction.

The Bregman divergence between the data matrix X and the reconstructed matrix U V T in our GaP factor model is related to the deviance of the Poisson model defined such as [START_REF] Landgraf | Generalized principal component analysis: Projection of saturated model parameters[END_REF] Dev(X,

U V T ) = -2 × log p(X | Λ = U V T ) -log p(X | Λ = X) ,
where log p(X | Λ) is the Poisson log-likelihood based on the matrix notation (1). We have Dev(X, U V T ) ∝ D(X | U V T ), thus the deviance can be used to quantify the quality of the model.

Regarding PCA, the percentage of explained variance is a natural and unequivocal quantification of the quality of the representation. We introduce a criterion that we call percentage of explained deviance that is a generalization of the percentage of explained variance to our GaP factor model. However, since our models are not nested for increasing K, the definition of this criterion appears non trivial. To assess the quality of our model, we propose to define the percentage of explained deviance as:

%dev = log p(X | Λ = U V T ) -log p(X | Λ = 1 n X) log p(X | Λ = X) -log p(X | Λ = 1 n X) (3)
where U V T is the predicted reconstructed matrix in our model, 1 n is a column vector filled with 1 and X is a row vector of size m storing the column-wise average of X.

We use two baselines: (i) the log-likelihood of the saturated model, i.e. log p(X | Λ = X) (as in the deviance), which corresponds to the richest model and (ii) the loglikelihood of the model where each Poisson intensities λ ij is estimated by the average of the observations in the column j, i.e. log p(X | Λ = 1 n X), which is the most simple model that we could use. This formulation ensures that the ratio %dev lies in [0; 1]. An interesting point is that if we assume a Gaussian distribution on the data, the percentage of explained deviance is exactly the percentage of explained variance from PCA (c.f. Section S.1), which makes our criterion suitable for to compare different factor models.

Choosing the dimension of the latent space

As noticed in the previous section, the GaP factor model with an increasing number K of factors are not nested (the model associated to the NMF presents the same properties). Consequently, testing different values of K requires to fit different models (contrary to PCA for instance). We choose the number of factors by fitting a model with a large K and verifying how the matrix U 1:k ( V 1:k ) T reconstructs X depending on k = 1, . . . , K with a rule of thumb based on the "elbow" shape of the fitting criterion . This approach is for instance widely used in PCA by checking the proportion of variability explained by each components, see Friguet (2010, p.96) for a review of the different criteria to choose K in this context. Here we use the deviance, or equivalently the Bregman divergence k → D X | U 1:k ( V 1:k ) T to find the latent dimension from where adding new factors does not improve D X | U 1:k ( V 1:k ) T . This determination is however not always unambiguous and may sometimes lead to some over-fitting, i.e. when considering too much factors. In addition, when focusing on data visualization, we generally set K = 2.

Model inference using a variational EM algorithm

Our goal is to infer the posterior distributions over the factors U and V depending on the data X. To avoid using the heavy machinery of MCMC [START_REF] Nathoo | Comparing variational Bayes with Markov chain Monte Carlo for Bayesian computation in neuroimaging[END_REF] to infer the intractable posterior of the latent variables in our model, we use the framework of variational inference [START_REF] Hoffman | Stochastic Variational Inference[END_REF]. In particular, we extend the version of the variational EM algorithm [START_REF] Beal | The variational Bayesian EM algorithm for incomplete data: With application to scoring graphical model structures[END_REF] proposed by [START_REF] Dikmen | Maximum marginal likelihood estimation for nonnegative dictionary learning in the Gamma-Poisson model[END_REF] in the context of the standard Gamma-Poisson factor model to our sparse and zero-inflated GaP model. Figure S.1 in Supp. Mat. gives an overview of the variational framework.

Definition of variational distributions

In the variational framework, the posterior p(Z, U, V , S, D | X) is approximated by the variational distribution q(Z, U, V , S, D) regarding the Kullback-Leibler divergence [START_REF] Hoffman | Stochastic Variational Inference[END_REF], that quantifies the divergence between two probability distributions. Since the posterior is not explicit, the inference of q is based on the optimization of the Evidence Lower Bound (ELBO), denoted by J(q) and defined as:

J(q) = E q [log p(X, Z, U, V , S, D)] -E q [log q(Z, U, V , S, D)] , (4)
that is a lower bound on the marginal log-likelihood log p(X). In addition, maximizing the ELBO J(q) is equivalent to minimizing the KL divergence between q and the posterior distribution of the model [START_REF] Hoffman | Stochastic Variational Inference[END_REF]. To derive the optimization, q is assumed to lie in the mean-field variational family, i.e. (i) to be factorisable with independence between latent variables and between observations and (ii) to follow the conjugacy in the exponential family, i.e. to be in the same exponential family as the full conditional distribution on each latent variables in the model. Thanks to the first assumption, in our model, the variational distribution q is defined as follows:

q(Z, U, V , S, D) = n i=1 m j=1 q (Z ijk ) k | (r ijk ) k × n i=1 K k=1 q(U ik | a ik ) × m j=1 K k=1 q(V jk | b jk ) × m j=1 K k=1 q(S jk | p S jk ) × n i=1 m j=1 q(D ij | p D ij ) (5)
where (r ijk ) k , a ik , b jk , p S jk and p D ij are the parameters of the variational distribution regarding (Z ijk ) k , U ik , V jk , S jk , D ij respectively. Then we need to precise the full conditional distributions of the model before defining the variational distributions more precisely.

Approximate posteriors

To approximate the (intractable) posterior distributions, variational distributions are assumed to lie in the same exponential family as the corresponding full conditionals and to be independent such that:

Z ij q ∼ M (r ijk ) k U ik q ∼ Γ(a ik,1 , a ik,2 ) V jk q ∼ Γ(b jk,1 , b jk,2 ) S jk q ∼ B(p S jk ) D ij q ∼ B(p D ij ),
where q ∼ denotes the variational distribution. The strength of our approach is the resulting explicit approximate distribution on the loadings that induces sparsity:

V jk |S jk q ∼ (1 -S jk ) × δ 0 + S jk × Γ(b jk,1 , b jk,2 ),
In the following, the derivation of variational parameters involves the moments and log-moments of the latent variables regarding the variational distribution. Since the distributions q is fully determined, these moments can be directly computed. For the sake of simplicity, we will use notation

U ik = E q [U ik ] and log U ik = E q [log U ik ]
(collected in the matrices U and log U respectively), with similar notations for other hidden variables of the model (V jk , D ij , S jk , Z ijk ).

Derivation of variational parameters

In order to find a stationary point of the ELBO, J(q) is differentiated regarding each variational parameter separately. The formulation of the ELBO regarding each parameter separately is based on the corresponding full conditional, e.g. p(U ik | -), p(V jk | -) and p (Z ijk ) k | -. The partial formulation are therefore respectively:

J(q) a ik = E q [log p(U ik | -)] -E q [log q(U ik ; a ik )] + cst J(q) b jk = E q [log p(V jk | -)] -E q [log q(V jk ; b jk )] + cst J(q) (r ijk ) k = E q log p (Z ijk ) k | - -E q log q (Z ijk ) k ; (r ijk ) k + cst
Similar formulations can be derived regarding parameters p D ij and p S jk . Therefore, the ELBO is explicit regarding each variational parameter and the gradient of the ELBO J(q) depending on the variational parameters a ik , b jk , r ijk , p D ij and p S jk respectively can be derived to find the coordinate of the stationary point (corresponding to a local optimum). In our factor model all full conditionals are tractable (c.f. Section S.4.1 in Supp. Mat.). In practice, thanks to the formulation in the exponential family, the optimum value for each variational parameter corresponds to the expectation regarding q of the corresponding parameter of the full conditional distribution (see [START_REF] Hoffman | Stochastic Variational Inference[END_REF]. Thus the coordinates of the ELBO's gradient optimal point are explicit. We mentioned that distributions with a mass at 0 (zero-inflated Poisson or sparse Gamma) lie in the exponential family [START_REF] Eggers | On Statistical Methods for Zero-Inflated Models[END_REF] and the general formulation from [START_REF] Hoffman | Stochastic Variational Inference[END_REF] remains valid. Detailed formulations of update rules regarding all variational parameters are given in Supp. Mat. (Section S.4.2).

Variational EM algorithm

We use the variational-EM algorithm [START_REF] Beal | The variational Bayesian EM algorithm for incomplete data: With application to scoring graphical model structures[END_REF] to jointly approximate the posterior distributions and to estimate the hyper-parameters Ω = (α, β, π S , π D ). In this framework, the variational inference is used within a variational E-step, in which the standard expectation of the joint likelihood regarding the posterior

E[p(X, U, V , S, D ; Ω)|X] is approximated by E q [p(X, U, V , S, D ; Ω)].
Then the variational M-step consists in maximizing E q [p(X, U, V , S, D ; Ω)] w.r.t. the hyper-parameters Ω. In the variational-EM algorithm, we have explicit formulations of the stationary points regarding variational parameters (E-step) and prior hyper-parameters (M-step) in the model, thus we use a coordinate descent iterative algorithm (see [START_REF] Wright | Coordinate descent algorithms[END_REF], for a review) to infer the variational distribution. Detailed formulations of update rules regarding all prior hyper-parameters are given in Supp. Mat. (Section S.4.3).

Initialization of the algorithm

To initialize variational and hyper-parameters of the model, we sample U and V from Gamma distributions such that X ij k U ik V jk on average. The Gamma variational and hyper parameters are initialized from these values following the update rules detailed in Supp. Mat. Section S.4.2. Dropout probabilities p D ij and π D j are initialized by 1/n i 1 {X ij >0} , i.e. the proportion of non-zero in the data for the corresponding gene. To initialize the sparsity probabilities p S jk and π S j , we use a heuristic based on the variance of the corresponding gene j. Assuming that genes with low variability will have less impact on the structure embedded in the data, we propose a starting value such that

P (0) j = 1 -exp (-s j / m j ) , (6)
where m j is the mean of the non-null observations for gene j and s j its standard deviation (including null values). The scaling is better when removing the null values to compute the mean. This quantity adapts to the empirical variance of the observations, and will be close to 0 for genes with low variance, and close to 1 for genes with high variability.

Empirical study of pCMF

All codes are available on a public repository for reproducibility2 . We compare our method with standard approaches for unsupervised dimension reduction: the Poisson-NMF [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], applied to raw counts (model-based matrix factorization approach based on the Poisson distribution); the PCA [START_REF] Pearson | LIII. On lines and planes of closest fit to systems of points in space[END_REF] and the sparse PCA [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF], based on an 1 penalty in the optimization problem defining the PCA to induce sparsity in the loadings V, both applied to log counts. We use sparse methods (sparse PCA, sparse pCMF) with a re-estimation step on the selected genes. We will refer to them as (s)PCA and (s)pCMF in the results respectively, to differentiate them from sparse PCA and sparse pCMF (without re-estimation), PCA and pCMF (without the sparsity layer). In addition, we use the Zero-Inflated Factor Analysis (ZIFA) by [START_REF] Pierson | ZIFA: Dimensionality reduction for zero-inflated singlecell gene expression analysis[END_REF], a dimension reduction approach that is specifically designed to handle dropout events in single-cell expression data (based on a zero-inflated Gaussian factor model applied to log-transformed counts). We present quantitative clustering results and qualitative visualization results on simulated and experimental scRNA-seq data. We also compare our method with t-SNE that is commonly used for data visualization (van der [START_REF] Van Der Maaten | Visualizing Data using t-SNE[END_REF]. It requires to choose a "perplexity" hyper-parameter that cannot be automatically calibrated, thus being less appropriate for a quantitative analysis. In the following, we always choose the perplexity values that gives the better clustering results.

Simulated data analysis

To generate synthetic data we follow the hierarchical Gamma-Poisson framework as adopted by others [START_REF] Zappia | Splatter: Simulation of single-cell RNA sequencing data[END_REF]. Details are provided in the Supp. Mat. (Section S.5). We generate synthetic multivariate over-dispersed counts, with n = 100 individuals and m = 800 recorded variables. We artificially create clusters of individuals (with different level of dispersion) and groups of dependent variables. We we set different levels of zero-inflation in the data (i.e. low or high probabilities of dropout events, corresponding to random null values in the data), and some part of the m variables are generated as random noise that do not induce any latent structure. Thus, we can test the performance of our method in different realistic data configurations, the range of our simulation parameters being comparable to other published simulated data [START_REF] Zappia | Splatter: Simulation of single-cell RNA sequencing data[END_REF].

We train the different methods with K = 2 to visualize the reconstructed matrices U and V (c.f. Section 2). To assess the ability of each method to retrieve the structure of cells or genes, we run a κ-mean clustering algorithm on U and V respectively (with κ = 2) and we measure the adjusted Rand Index [START_REF] Rand | Objective Criteria for the Evaluation of Clustering Methods[END_REF] the accordance between the predicted clusters and original groups of cells or genes. Regarding our approach pCMF, we use log U and log V for data visualization and clustering because the log is the canonical link function for Gamma models. In addition, we also compute the percentage of explained deviance associated to the model to assess the quality of the reconstruction. Regarding the PCA (sparse or not) and ZIFA, we use the standard explained variance criterion (c.f. Section 2.1).

Clustering in the observation space

Effect of zero-inflation and cell representation. We first assess the robustness of the different methods to the level of dropout or zero-inflation (ZI) in the data. We generate data with 3 groups of observations (c.f. Supp. Mat. Section S.5) with a wide range of dropout probabilities. Figure 1a shows that (s)pCMF adapts to the level of dropout in the data and recovers the original clusters (high adjusted Rand Index) even with dropouts. Despite comparable performance with low dropout, Poisson-NMF and ZIFA are very sensitive to the addition of zeros. In addition, methods based on transformed counts like (s)PCA and t-SNE perform poorly, as they do not account for the specificity of the data (discrete, over-dispersed, (O'Hara and Kotze, 2010)).

Effect of noisy genes and gene representation. To quantify the impact of noisy genes on the retrieval of the clusters, we consider data generated with different proportions of noisy genes (genes that do not induce any structure in the data). We generate data with three groups of genes: two groups inducing some latent structure and one group of noisy genes (c.f. Supp. Mat. Section S.5). Figure 2 shows that (s)pCMF identifies correctly the clusters of genes, including the set of noisy genes, contrary to other approaches. This point shows that our approach correctly identifies the genes that support the lower-dimensional representation.

In addition to the clustering results, we compared the selection accuracy of the only two methods (sPCA, spCMF) that perform variable selection (Supp. Mat. Figure S.2). A selected gene is a gene that contributes to any latent dimension (any V jk = 0). Sparse pCMF performs better than sparse PCA for various latent dimension K even for high levels of noisy genes. Sparse PCA shows better selection accuracy when the proportion of noisy genes is low. This point suggests that sparse pCMF would be less sensitive to gene pre-filtering when analysing scRNA-seq data, which corresponds to a removal of noisy genes and is generally crucial [START_REF] Soneson | Bias, robustness and scalability in single-cell differential expression analysis[END_REF]. Details about data generation and additional data configuration regarding Figures 1 and 2 can be found in Supp. Mat. (Section S.7, especially Figures S.3 and S.4).

Data visualization

Data visualization is central in single-cell transcriptomics for the representation of high dimensional data in a lower dimensional space, in order to identify groups of cells, or to illustrate the cells diversity (e.g. [START_REF] Llorens-Bobadilla | Single-Cell Transcriptomics Reveals a Population of Dormant Neural Stem Cells that Become Activated upon Brain Injury[END_REF][START_REF] Segerstolpe | Single-Cell Transcriptome Profiling of Human Pancreatic Islets in Health and Type 2 Diabetes[END_REF]. In the matrix factorization framework, we represent observation (cell) coordinates and variable (gene) contributions: resp. ( u i1 , u i2 ) i=1,...,n and ( v i1 , v i2 ) i=1,...,n (or their log transform for pCMF) when the dimension is K = 2 (see Section 2). We consider the same simulated data as previously (n = 100, m = 800, with three groups of cells, two groups of relevant genes and the set of noisy genes). Our visual results are consistent with the previous clustering results both regarding cell and gene visualization. In this challenging context (high zero-inflation and numerous noisy variables), by using pCMF, we are able to graphically identify the groups of individuals (cells) in the simulated zero-inflated count data (c.f. Figure 3). On the contrary, the 2-D visualization is not successful with PCA, ZIFA, Poisson-NMF and t-SNE, illustrating the interest of our data-specific approach. This point supports our claim that using data-specific model improves the quality of the reconstruction in the latent space.

In addition, linear projection methods (NMF, PCA, pCMF, ZIFA) can be used to visualize the contribution of genes to the principal axes (c.f. Figure 4). Thanks to sparsity constraints, the contribution of noisy genes are mostly set to 0 for sparse pCMF. Surprisingly, this selection is not efficient in the case of sparse PCA, indicating a lack of calibration of the sparsity constraint when data are counts. In comparison, Poisson NMF and ZIFA (not sparse) do not identity the cluster of noisy genes as clearly as sparse pCMF.

To quantify the model quality of the different methods on simulated data, we used the deviance associated to each method (c.f. Section 2.1). Figures 1b and2b shows that the dimension reduction proposed by pCMF has excellent fit to the data regardless the level of dropout or the proportion of noisy genes, as compared with other methods.

Additional comparisons of computational time show that PCA is fastest method (but with low performance), whereas (sparse) pCMF is faster than ZIFA and sparse PCA, with increased performance (cf Supp. Mat. Section S.7).

Analysis of single-cell data

We now illustrate the performance of pCMF on different recent and large single-cell RNA-seq datasets that are publicly availbale: the [START_REF] Baron | A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Interand Intra-cell Population Structure[END_REF] dataset, the goldstandard and silverstandard datasets used in [START_REF] Freytag | Comparison of clustering tools in R for medium-sized 10x Genomics single-cell RNA-sequencing data[END_REF] (we used the silverstandard dataset 5 which was the largest). We also consider an older and smaller dataset from [START_REF] Llorens-Bobadilla | Single-Cell Transcriptomics Reveals a Population of Dormant Neural Stem Cells that Become Activated upon Brain Injury[END_REF] which is interesting because it discribes a continuum of activation in Neural Stem Cells (NSC). All datasets are available here3 with the codes. More details about their origins are given in Supp. Mat. Section S.7.5. We consider large datasets with ∼ 1000 or ∼ 10000 cells to test the ability of our approach to face the expected increase of data volume in the next couple of years.

We present some quantitative results about clustering and data reconstruction (c.f. Table 1) and the corresponding qualitative results about cell visualization (c.f. Figures 5 andS.7 to S.9 in Supp. Mat.) and gene visualization (c.f. Figures S.10 to S.13 in Supp. Mat.). For each dataset (except the one from Llorens-Bobadilla et al., 2015, where we used their pre-filtering), we use the same pipeline, we filter out genes expressed in less than 5% of the cells. In a second step, we remove genes for which the variance heuristic defined in Equation ( 6) is low. In practice we removed genes for which P (0) j ≤ 0.2. Our idea was to remove uninformative genes, since pre-filtering is crucial [START_REF] Soneson | Bias, robustness and scalability in single-cell differential expression analysis[END_REF] in such data, but also to reduce the number of genes to reduce the computation cost, in particular for ZIFA. Then we compare (s)pCMF, PCA, ZIFA and t-SNE. We discarded (s)PCA because the sparse PCA is computationally expensive (c.f. Supp. Mat. Section S.7.3) due to the required cross-validation.

A general empirical property is that clustering accuracy dicreases for all methods when the number of groups of cells increases. However, our approach (s)pCMF produces a better (or as good) view of the cells regarding clustering purpose in every examples, since the adjusted Rand Index is higher (c.f. Table 1). We observe the same trend regarding the quality of the reconstruction since the explained deviance is generally also higher for (s)pCMF. Data visualization is not always clear (c.f. Figures S.7 and S.8), especially when the number of groups is large as in [START_REF] Baron | A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Interand Intra-cell Population Structure[END_REF] or [START_REF] Freytag | Comparison of clustering tools in R for medium-sized 10x Genomics single-cell RNA-sequencing data[END_REF] silverstandard, however it is possible to clearly identify large clusters of cells in the (e.g. beta cells in [START_REF] Baron | A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Interand Intra-cell Population Structure[END_REF] or CD14+ Monocyte in Freytag et al.

(2018) silverstandard) with our method (and some of the others). On the goldstandard dataset from [START_REF] Freytag | Comparison of clustering tools in R for medium-sized 10x Genomics single-cell RNA-sequencing data[END_REF], the difference regarding cells representation between the different approaches is more visual (c.f. Figure 5), where our approach is the only one that is able to distinctly identify the three populations of cells. On the Llorens-Bobadilla et al. (2015) dataset, our approach clearly highlights this continuum of activation presented in their paper, which can also be seen with ZIFA, but is not as much clear with PCA and t-SNE.

Regarding gene visualization (c.f. Figures S.10 to S.13 in Supp. Mat.), we compare the representation of sparse pCMF to PCA, ZIFA (and not t-SNE since it does not jointly learn U and V). The interest of sparsity for gene selection is to highlight more precisely the genes that contribute to the latent representation. For each dataset, it is possible to detect which genes are important for each latent dimension: some are null on both (e.g. uninformative genes), some contribute to a single dimension, some contribute to both. This pattern is not as clear with methods that do not implement any sparsity layers.

These different points show the interest of our approach to analyze recent single-cell RNA-seq datasets, even large ones. Empirical properties studied on simulations are confirmed on experimental data: providing a dimension reduction method adapted to single cell data, where the sparsity constraints is powerful to represent complex single cell data. In addition, our heuristic to assume gene importance based on their variance appears to be efficient (i) to perform a rough pre-filtering to reduce the dimension and (ii) to discriminate between noisy genes and informative ones directly in the sparse pCMF algorithm. In addition, it appears that our method is fast compared to ZIFA for instance, since it takes less than two minutes on the different examples to run (s)pCMF (sparse pCMF + re-estimation), on a 16-core machine, whereas ZIFA can take between 5 and 25 minutes on the same architecture. (adjusted Rand Index). Each scRNA-seq dataset is characterized by the number of cells, the number of genes used in the analysis (we specify the original number before the pre-filtering step) and by the number of original groups. The adjusted Rand Index compares clusters found by a κ-means algorithm (applied to U with κ = nb group) and the original groups of cells.

Conclusion

In this work, we provide a new framework for dimension reduction in unsupervised context. In particular, we introduce a model-based matrix factorization method specifically designed to analyze single-cell RNA-seq data. Matrix factorization allows to jointly construct a lower dimensional representation of cells and genes. Our probabilistic Count Matrix Factorization (pCMF) approach accounts for the specificity of these data, being zero-inflated and over-dispersed counts. In other words, we propose a generalized PCA procedure that is suitable for data visualization and clustering. The interest of our zero-inflated sparse Gamma-Poisson factor model is to replace the variance-based formulation of PCA, associated to the Euclidean geometry and the Gaussian distribution, with a metric (based on Bregman divergence) that is adapted to scRNA-seq data characteristics.

Analyzing single-cell expression profiles is a huge challenge to understand the cell diversity in a tissue/an organism and more precisely for characterizing the associated gene activity. We show on simulations and experimental data that our pCMF approach is able to catch the underlying structure in zero-inflated over-dispersed count data. In particular, we show that our method can be used for data visualization in a lower dimensional space or for preliminary dimension reduction before a clustering step. In both cases, pCMF performs as well or out-performs state-of-the-art approaches, especially the PCA (being the gold standard) or more specific methods such as the NMF (count based) or ZIFA (zero-inflation specific). In particular, pCMF data representation is less sensitive to the choice of the latent dimension K regarding clustering results, which supports the interest of our approach for data exploration. It appears (through the explained deviance criterion that we introduced) that the reconstruction learned by pCMF better represents the variability in the data (compared to PCA or ZIFA). In addition, pCMF can select genes that explain the latent structure in the data, thanks to a sparsity layer which does not require any parameter tuning.

An interesting direction to improve pCMF would be to integrate covariables or con-founding factors in the Gamma-Poisson model, for instance to account for technical effect in the data or for data normalization. A similar framework based on zero-inflated Negative Binomial distribution was proposed by [START_REF] Risso | ZINB-WaVE: A general and flexible method for signal extraction from single-cell RNA-seq data[END_REF], and it could be extended to our framework of matrix factorization

Supplementary materials

S.1 Generalization of explained variance

In the Gaussian framework, we assume that X ij ∼ N (µ ij , 1) since data are preliminary centered and scaled in PCA. Under the assumptions of independence between observations, the log-likelihood is then in matrix notation:

log p(X | M) = n i=1 m j=1 log p(x ij | µ ij ) = n i=1 m j=1 (x ij -µ ij ) 2 = X -M 2 F where M = [µ ij ]
is the matrix of Gaussian expectation and • 2 F is the squared Frobenius norm. In the generalized PCA framework [START_REF] Collins | A generalization of principal components analysis to the exponential family[END_REF], we are looking for U ∈ R n×K and V ∈ R m×K such that M = UV T . Thanks to [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] theorem, best U and V minimizing the Frobenius norm between X and UV T are given by Singular Value Decomposition (SVD) of X, and optimal U exactly corresponds to the principal components from the PCA, which highlights the link between PCA, SVD and Gaussian framework.

In this Gaussian framework, the explained deviance defined in Equation (3) can be rewritten as

%dev = log p(X | M = U V T ) -log p(X | M = 1 n X) log p(X | M = X) -log p(X | M = 1 n X) ,
since the saturated model corresponds to M = X in this case. It follows that

%dev = X -U V T 2 F -X -1 n X 2 F X -X 2 F -X -1 n X 2 F .
In addition, we have that X = 0 thanks to the pre-centering, and the formulation becomes:

%dev =1 - X -U V T 2 F X 2 F , =1 - rk(X) k=K+1 σ 2 k rk(X) =1 σ 2 = K k=1 σ 2 k rk(X) =1 σ 2
, where rk(X) is the rank of X and σ 1 > . . . > σ rk(X) the singular values of X (given by the SVD). The criterion corresponds exactly to the percentage of explained variance computed in PCA. Thus our percentage of explained deviance can be viewed as a generalization of this criterion to other distributions in the exponential family.

S.2 Identifiability issues

S.2.1 Factor order

Gamma-Poisson factor model suffers from an identifiability issue regarding the order of factors. Unlike PCA, the components of model-based factor models are not orthogonal and can not be ordered naturally since the associated likelihood is identifiable up to a permutation of factors. Thus we propose an ordering defined by the cumulative Bregman divergence:

k → D X | U 1:k ( V 1:k ) T .
In addition, we mention that the different GaP factor models are not nested when the dimension K increases (as in the NMF), thus the factor estimates should be all computed for every choice of dimension K, contrary to PCA.

sub

S.3 Scaling effect in GaP factor model

As stated in [START_REF] Dikmen | Maximum marginal likelihood estimation for nonnegative dictionary learning in the Gamma-Poisson model[END_REF], GaP factor models suffer from identifiability issues, due to the scaling of the Gamma prior parameters α and β. Indeed, considering

α * k,2 = η k α k,2 and β * k,2 = (η k ) -1 β k,2
for fixed values η k , and using the scaling property of the Gamma distribution: if

U ik ∼ Gamma(α k,1 , α k,2 ) then η k U ik ∼ Γ(α k,1 , η -1 k α k,2
). We show (c.f. below) that the joint log-likelihood regarding UH -1 and VH with H = diag(η k ) k=1:K verifies:

log p(X, UH -1 , VH | α 1 , Hα 2 , β 1 , H -1 β 2 ) = log p(X, U, V | α 1 , α 2 , β 1 , β 2 ) + (n -p) k log(η k ) (S.1)
When n = p, there is an identifiability issue regarding the scaling of the parameters α k,2 and β k,2 , because different values lead to the same joint log-likelihood. In such case, a solution will be to fix the scale parameters α k,2 and β k,2 to avoid the scaling effect. When n = p, the only problem is a potential solution with infinite norm with α k,2 → 0 and β k,2 → ∞ or vice-versa (c.f. [START_REF] Dikmen | Maximum marginal likelihood estimation for nonnegative dictionary learning in the Gamma-Poisson model[END_REF]. When considering zeroinflation or sparsity in the model, Equation (S.1) holds regarding the parameters of the Gamma prior distributions and we have to consider the same precaution. However, in practice we did not encounter such sequence of diverging parameters.

Proof of Equation (S.1). We set, α * k,2 = η k α k,2 and β * k,2 = (η k ) -1 β k,2 for fixed values η k > 0. We use the scaling property of the Gamma distribution: if

U ik ∼ Gamma(α k,1 , α k,2 ) then η k U ik ∼ Γ(α k,1 , (η k ) -1 α k,2
). The joint log-likelihood regarding UH -1 and VH with H = diag(η k ) k=1:K is then:

log p(X, UH -1 , VH | α 1 , Hα 2 , β 1 , H -1 β 2 ) = i,j,k log p x ij | {(η k ) -1 u ik , η k v jk } k=1:K + i,k log p η -1 k u ik ; α k,1 , η k α k,2 + j,k log p η k v jk ; β k,1 , (η k ) -1 β k,2 = log p(X | U, V) + log p(U ; α 1 , α 2 ) + n i=1 K k=1 log(η k ) + log p(V ; β 1 , β 2 ) + p j=1 K k=1 -log(η k ) = log p(X, U, V | α 1 , α 2 , β 1 , β 2 ) + (n -p) k log(η k )

S.4 Variational inference algorithm

Figure S.1 describes the variational framework (for the GaP factor model) that we extended to develop our approach.

S.4.1 Full conditional distributions

In our factor model all full conditionals are tractable. Thanks to the Gamma-Poisson conjugacy, the full conditionals of U ik and V jk are Gamma distributions. The proof is based on the Bayes rule and the distribution of the latent variables Z, that are actually necessary to derive p(U ik | -) and p(V jk | -). The full conditional of the vector Z ij is also explicit, being a Multinomial distribution [START_REF] Zhou | Beta-negative binomial process and Poisson factor analysis[END_REF] when D ij = 0 and deterministic null when

D ij = 0, i.e. (Z ijk ) k | -∼ D ij M X ij , (ρ ijk ) k .
Here

The model ( * )

X ij = k Z ijk Z ijk | U ik , V jk ∼ P(U ik V jk ) U ik ∼ Γ(α k,1 , α k,2 ) V jk ∼ Γ(β k,1 , β k,2 ) -→ Intractable posterior -→ Variational framework   Optimization of J(q) ←- Approximate the posterior by the distrib. q Variational distribution U ik q ∼ Γ(a ik,1 , a ik,2 ) V jk q ∼ Γ(b jk,1 , b jk,2 ) (Z ijk ) k q ∼ M X ij , (r ijk ) k Complete conditional U ik | -∼ Γ η ik (-) V jk | -∼ Γ η jk (-) (Z ijk ) k | -∼ M X ij , (ρ ijk ) k
Inference of q ( * ) with conditional independence between the Z ijk 's and independence between the U ik 's and V jk 's the Multinomial probabilities (ρ ijk ) k depend on (S jk , U ik , V jk ) k , and quantify the prior contribution of factor k to the observations X ij , i.e.

ρ ijk = S jk U ik V jk S j U i V j .
This point justifies why the variational distribution is based on the vector Z ij instead of taking each Z ijk separately. Note that if the S jk are null for all k or if D ij = 0 (i.e. X ij = 0), the vector (Z ijk ) k is deterministic and takes null values. We summarize the full conditionals in the sparse ZI-GaP factor model regarding U ik , V jk and (Z ijk ) k , that are defined such as:

U ik | -∼ Γ(α k,1 + j D ij S jk Z ijk , α k,2 + j D ij S jk V jk ) , V jk | -∼ Γ(β k,1 + i D ij S jk Z ijk , β k,2 + i D ij S jk U ik ) , (Z ijk ) k | -∼ D ij M X ij , (ρ ijk ) k , (S.2)
Zero Inflation. Regarding the zero-inflation indicators, D ij is a binary variable, its distribution is either deterministic or Bernoulli. When the entry X ij is non null, D ij is certainly equal to one. When X ij = 0, the full conditional is explicit and the Bernoulli probability only depends on the prior over D ij and the probability that X ij is null. It can be formulated as follows:

p(D ij = 1 | -) ∝ π D j e -k S jk U ik V jk .
Sparsity and variable selection. The sparsity indicator S jk is also a binary variable and its full conditional is also an explicit Bernoulli distribution. It depends on the prior over S jk and the probability that gene j contributes to the components k, quantified by the joint distribution on (Z ijk ) i , thus:

p(S jk = 1 | -) ∝ π s j × i exp(-S jk U ik V jk ) (S jk U ik V jk ) Z ijk .

S.4.2 Derivation of variational parameters

Variational parameters of factors. We derive the stationary point formulation for the variational parameters regarding U ik and V jk , being explicitly (directly derived from the partial derivatives of J(q)):

a ik = α k,1 + j D ij S jk Z ijk , α k,2 + j D ij S jk V jk T b jk = β k,1 + S jk i D ij Z ijk , β k,2 + S jk i D ij U ik T ,
which generalizes formulations from standard GaP factor model [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF]. As for variable Z ijk , its posterior distribution depends on parameter r ijk with the relation log(r ijk ) = E q [log(ρ ijk )]. Hence, the variational distribution on (Z ijk ) k naturally depends on the selection indicator S jk (since our model focuses on loadings selection).

In particular, the variational parameter r ijk depends on S jk , through a specific term E q [log(S jk V jk )] that is computed using the variational distribution of S jk (a Bernoulli distribution of parameter p S jk ). To proceed, we introduce S jk , the discretized predictor of S jk such that S jk = 1 {p S jk >τ } , where τ is a threshold specified by the user (for instance 0.5). Then, the formulation of the optimal variational parameter r ijk is approximated by:

r ijk = S jk exp log U ik + log V jk S j exp log U i + log V j .
Variational dropout proportion. Regarding the zero-inflated probabilities p D ij , when X ij = 0, the posterior is explicit since D ij = 1 with probability one. Hence, only the case X ij = 0 requires a variational inference. As stated previously, the full conditional is explicit and it is possible to derive and optimize the ELBO (based on the natural parametrization of the Bernoulli distribution in the exponential family). Eventually, p D ij is computed as:

logit(p D ij ) = logit(π D j ) -k S jk U ik V jk ,
where the Bernoulli prior probability π D j is corrected by E q [log P(X ij = 0)] to account for the probability of X ij being a true zero.

Variational Selection probability. Concerning the sparse indicator S jk , the natural parametrization of the Bernoulli distribution is based on the logit of the Bernoulli probability. Hence we can write an explicit formulation of the ELBO regarding p S jk based on the full conditional on S jk . Following this formulation, the stationary point p S jk verifies:

logit(p S jk ) = logit(π S j ) -i D ij U ik V jk + D ij Z ijk log U ik + log V jk .
This corresponds to a correction of the Bernoulli prior probability π S j , depending on the quantification of the contribution of gene j to component k in all individuals, i.e. E q [ i log p(Z ijk )].

S.4.3 Derivation of prior parameters

The hyper-parameters of priror distribution regarding U ik , V jk , D ij and S jk are updated within the M-step such that (respectively):

α k,1 =ψ -1 log α k,2 + 1 n i log U ij , α k,2 = α k,1 i U ij /n , β k,1 =ψ -1   log β k,2 + 1 p j log V ij   , β k,2 = β k,1 j V ij /p , π D j = 1 n i p D ij , π S j = 1 K k p S jk ,
where ψ is the digamma function, i.e. the derivative of the log-Gamma function. Its inverse is computed thanks to the method proposed in Minka (2000, appendix C). Recalling that, for a variable U ∼ Γ(α 1 , α 2 ), E[U ] = α 1 /α 2 and E[log U ] = ψ(α 1 )log α 2 , the update rule for the Gamma prior parameters on U ik corresponds to averaging the moments and log-moments of the variational distribution on U ik over i (similarly for V jk over j). Regarding the Bernoulli prior parameters π D j , the update rule is also an average of the corresponding variational parameter over i (similarly for π S j over k).

S.4.4 Algorithm

Our pCMF algorithm is summarized in Algorithm S.1. In the initialization step, each variational Gamma shape parameter a ik,1 and b jk,2 are sampled from a Gamma distribution [START_REF] Zhou | Beta-negative binomial process and Poisson factor analysis[END_REF]. Each variational Gamma rate parameter a ik,2 and b jk,2 are set to 1 (to avoid scaling effect between shape and rate parameters). Each variational dropout probability p D ij is initialized with the corresponding indicator δ 0 (X ij ). Each variational sparsity probability p S jk is initialized with the corresponding with the threshold value τ ∈ (0, 1). In addition, all prior hyper-parameters are initialized following update rules based on variational parameters (defined in Section 3.4 in the paper).

The convergence is assessed by computing the normalized gap between two successive parameter values across iterations. When the updates does not modify the values of the parameters, we can consider that we reach a fixed point and thus the optimum.

In addition, to overcome potential issue related to local optimum, the algorithm is run several times with different random initialization and the best seed (regarding the ELBO criterion) is kept.

To overcome this issue, we use the following generative process. The variables j = 1, . . . , p are first partitioned into two groups V 0 and V ∅ of respective sizes m 0 and m -m 0 (with m 0 ≤ m). The m 0 variables in V 0 will represent the pertinent variables for the lower dimensional representation, whereas variables in V ∅ will be considered irrelevant or noise. The matrix V will be a concatenation of two matrices V 0 and V ∅ :

V m×K =   V 0 V ∅  
The ratio m 0 /m sets the expected degree of sparsity in the model. In practice, we generate m 0 /m from a Beta distribution, so that in average m 0 /m ∈ {0.2, 0.4, 0.6, 0.8} corresponding to different proportions of noisy genes (between 20 and 80% of noisy genes).

To simulate dependency between recorded variables, we generate groups of variables in the set V 0 of pertinent variables. We use a similar strategy as the one used to simulate U. V 0 is partitioned into M balanced groups, denoted by V 1 , . . . , V M . We generate the corresponding matrix V 0 with blocks on the diagonal. Each block, denoted by B V,g contains m 0 /M rows and K/M columns: Each entry V jk in each block B V,g (g = 1, . . . , M ) is drawn from a Gamma distribution Γ(1, 1/β) with a rate parameter depending on β > 0. All entries V jk that are not in the blocks on diagonal are drawn from a Gamma distribution Γ(1, 1/((1 -θ v )β), where θ v ∈ (0, 1) quantifies how much the groups of individuals are distinct. Hence, each groups of individuals V g corresponds to a block B V,g . Again, this generation pattern requires that K > M . In practice, we fix β = 80, we use θ v = 0.8 and M = 2 groups of variables.

In addition, all V jk in V ∅ (noisy genes) are drawn from a Gamma distribution Γ(1, 1/((1θ v )β), so that E[V jk ] will not be structured according to groups.

Generation of X. The data are simulated according to their conditional Poisson distribution in the model i.e. P( k u ik v jk ). In practice, we want to consider zeroinflation in the model, thus we consider the Dirac-Poisson mixture and simulate X ij according to the following conditional distribution:

X ij | (U ik , V jk ) k , D ij ∼ (1 -D ij ) × δ 0 + D ij × P( k U ik V jk ) ,
where the dropout indicator D ij is drawn from a Bernoulli distribution B(π D j ), the proportion of dropout events is set by the probability π D j . To generate data without dropout events, we just have to set D ij = 1 for any couple (i, j), i.e. π D j = 1 for any j.

In practice, we fix K = 40, n = 100 and m = 800 to simulate our data. We generate different level of zero-inflation: π D j is drawn from a beta distribution so that in average it lies in {0.3, 0.5, 0.7, 0.9}.

S.6 Softwares

The Poisson-NMF is from the NMF R-package [START_REF] Gaujoux | A flexible R package for nonnegative matrix factorization[END_REF], ZIFA from the ZIFA Python-package [START_REF] Pierson | ZIFA: Dimensionality reduction for zero-inflated singlecell gene expression analysis[END_REF], the sparse PCA from the PMA Rpackage [START_REF] Witten | A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis[END_REF] and t-SNE from the Rtsne R-package [START_REF] Krijthe | Rtsne: T-Distributed Stochastic Neighbor Embedding using Barnes-Hut Implementation[END_REF]. Computation of adjusted Rand Index was done thanks to the mclust R-package [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. penalty parameter). t-SNE is slightly faster but requires numerous run with different values for the perplexity parameter (here the timing corresponds to a run for a single perplexity value). The PCA is the gold standard regarding running time thanks to the efficiency of its algorithm based on the Singular Value Decomposition (SVD) algorithm.

S.7 Additional results

S.7.1 Selection accuracy

Packages from where the different methods can be found are detailed in Section S.6.

Eventually, we mention that our method is available in an R-package, however our algorithms are implemented in interfaced C++ for computational efficiency.

S.7.4 Standard GaP versus our ZI sparse GaP factor model

Figure S.6 illustrates the interest of our zero-inflated sparse Gamma-Poisson factor model compared to the standard Gamma-Poisson factor model, especially in presence of dropout events and noisy genes. Our method pCMF based on our ZI sparse GaP factor model performs as well as the pCMF based on the standard GaP factor model when there is no dropout events in the data, independently from the proportion of noisy genes. In addition, when the level of zero-inflation is higher, we can see that the ZI-specific model outperforms the standard ones, highlighting the interest of our approach.

S.7.5 Additional scRNA-seq data analyses

The dataset from [START_REF] Baron | A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Interand Intra-cell Population Structure[END_REF] is available here4 . The goldstandard and silverstandard datasets used in [START_REF] Freytag | Comparison of clustering tools in R for medium-sized 10x Genomics single-cell RNA-sequencing data[END_REF] can be found here5 (we used the silverstandard dataset 5 which was the largest). The 3 previous datasets are stored based on the SingleCellExperiment R package [START_REF] Lun | SingleCellExperiment: S4 Classes for Single Cell Data. R package version 1[END_REF]. The dataset from Llorens-Bobadilla et al. We illustrate the performance of pCMF on a publicly available scRNA-seq datasets of neuronal stem cells [START_REF] Llorens-Bobadilla | Single-Cell Transcriptomics Reveals a Population of Dormant Neural Stem Cells that Become Activated upon Brain Injury[END_REF]. Neural stem cells (NSC) constitute an essential pool of adult cells for brain maintenance and repair. [START_REF] Llorens-Bobadilla | Single-Cell Transcriptomics Reveals a Population of Dormant Neural Stem Cells that Become Activated upon Brain Injury[END_REF] proposed a study to unravel the molecular heterogeneities of NSC populations based on scRNA-seq, and particularly focused on quiescent cells (qNSC). In their experiment, qNSC were transplanted in vivo in order to study their neurogenic activity. Following transplantation, 92 qNSC produced neuroblasts and olfactory neurons, whose transcriptome was compared with 21 astrocytes (CTX) and 27 transient amplifying progenitor cells (TAP). The authors used a PCA approach to reveal a continuum of "activation state", from astrocytes (low activation) to amplifying progenitor cells (TAP).

As stated before, we confront pCMF with other state-of-the-art approaches. The first visual result (c.f. Figure S.9) is that pCMF provides a slightly better representation of the continuum of activation described by Llorens-Bobadilla et al. than PCA and t-SNE, which probably reflects a better modeling of the biological variations that exist between activation states. In practice, t-SNE was not able to highlight the different clusters of cells. The results from ZIFA are consistent with pCMF representation, which is a confirmation that the signal of this continuous activation state is strong in these data. These qualitative results are confirmed by clustering quantitative results (c.f. Table 1 in the manuscript). The adjusted Rand Index computed after pCMF and ZIFA are similar and better than PCA.

S.7.5.2 Representation of cells

See Figures S.7 to S.9.

Figure 1 :

 1 Figure1: Adjusted Rand Index (1a) for the clustering on U versus the true groups of cells; and explained deviance (1b) depending on the probability used to generate dropout events. Average values and deviation are estimated across 50 repetitions.

Figure 2 :

 2 Figure2: Adjusted Rand Index (2a) for the clustering on V versus the true groups of genes; and explained deviance (2b) depending on the proportion of noisy genes. Average values and deviation are estimated across 50 repetitions.

Figure 3 :Figure 4 :

 34 Figure3: Representation of cells in a subspace of dimension K = 2. Here we have 60% of noisy variables, and dropout probabilities around 0.9.

Figure 5 :

 5 Figure 5: Analysis of the goldstandard scRNA-seq data from Freytag et al. (2018), 925 cells, 8580 genes. Visualization of the cells in a latent space of dimension 2.

Figure S. 1 :

 1 Figure S.1: Variational inference to approximate the posterior of the model, based on the optimization of the ELBO that required to derive the full conditional. The notation q ∼ refers to the variational distribution.

  Figure S.2: Selection accuracy depending on the dimension K (S.2a) with a proportion of noisy genes set to 60% and the proportion of noisy genes (S.2b) with K set to 2. Average values and deviation are estimated across 50 repetitions.

  Figure S.3: Adjusted Rand Index (S.3a) for the clustering on U versus the true groups of cells; and explained deviance (S.3b) depending on the probability used to generate dropout events, for different levels of separability between cell groups. Average values and deviation are estimated across 50 repetitions.

  Figure S.6: Adjusted Rand Index comparing clusters found by a κ-means algorithm (applied to U with κ = 2) and the original groups of individuals, depending on the number of individual groups in the data, for different levels of zero-inflation and different proportion of noisy variables in the data. The number of components is set to K = 10. Data are generated with n = 100, m = 1000. Average values and deviation are estimated across 100 repetitions.

  (2015) was available as supplementary data of their paper. They kindly shared with us the information about cell tags. S.7.5.1[START_REF] Llorens-Bobadilla | Single-Cell Transcriptomics Reveals a Population of Dormant Neural Stem Cells that Become Activated upon Brain Injury[END_REF]

  Figure S.7: Analysis of the scRNA-seq dataset from Baron et al. (2016), 1186 cells, 6080 genes. Visualization of the cells in a latent space of dimension 2.

  Figure S.8: Analysis of the silverstandard 5 scRNA-seq dataset from Freytag et al. (2018), 8352 cells, 4547 genes. Visualization of the cells in a latent space of dimension 2.

  quantifying

	adj. RI	0.00 0.25 0.50 0.75 1.00		explained deviance	0.00 0.25 0.50 0.75 1.00		method NMF (s)pCMF
		0.4	0.6	0.8	0.4	0.6	0.8
		drop-out probability	drop-out probability
		(a)

Table 1 :

 1 Performance of the different methods regarding quality of reconstruction (percentage of explained deviance) and clustering

		nb	nb genes	prop. 0	nb		(s)pCMF PCA	ZIFA t-SNE
		cells (before pre-filter.)		group
	Baron et al. (2016)	1886	6080	80.9%	13	adj. RI	21.2%	14.3% 15.4% 14.2%
			(14878)			%dev	73.2%	41.6% 53.5%	/
	Freytag et al. (2018)	925	8580	39.5%	3	adj. RI	81.3%	60.1% 56.8% 60.5%
	goldstandard		(58302)			%dev	55.7%	65.6% 48.6%	/
	Freytag et al. (2018)	8352	4547	86.3%	11	adj. RI	24.2%	16.2% 19.8% 24.8%
	silverstandard 5		(33694)			%dev	70.0%	55.1%	/	/
	Llorens-Bobadilla et al. (2015) 141	13826	64.8%	6	adj. RI	40.1%	25.3% 38.3% 29.8%
			(43309)			%dev	64.4%	34.8% 42.6%	/

https://github.com/gdurif/pCMF_experiments

https://github.com/gdurif/pCMF_experiments

https://hemberg-lab.github.io/scRNA.seq.datasets/mouse/pancreas/

https://github.com/bahlolab/cluster_benchmark_data

Funding

This work was supported by the french National Resarch Agency (ANR) as part of the "ABS4NGS" project [ANR-11-BINF-0001-06] and as part of the "MACARON" project [ANR-14-CE23-0003], and by the European Research Council as part of the ERC grant Solaris. It was performed using the computing facilities of the computing center LBBE/PRABI.

Our work is implemented in the pCMF R-package 1 .

1 https://github.com/gdurif/pCMF

Algorithm S.1: Variational EM algorithm

Data: count matrix X Result: factors U and V Initialization: random initialization of variational parameters, prior hyper-parameters are updated accordingly while No convergence do Update variational parameters (see Section 3.3 in the paper); Update prior pararmeters (see Section 3.4 in the paper); end

S.5 Data generation

We set the hyper-parameters (α k,1 , α k,2 ) k and (β k,1 , β k,2 ) k of the Gamma prior distributions on U ik and V jk to generate structure in the data, i.e. groups of individuals and groups of variables.

Generation of U.

In practice, individuals i = 1, . . . , n are partitioned into N balanced groups, denoted by U 1 , . . . , U N . To do so, we generate a matrix U with blocks on the diagonal. Each block, denoted by B U,g contains n/N rows and K/N columns. Each entry U ik in each block B U,g (g = 1, . . . , N ) is drawn from a Gamma distribution Γ(1, 1/α g ) with a rate parameter depending on α g > 0 (different for each group). All entries U ik that are not in the diagonal blocks of U are drawn from a Gamma distribution Γ(1, 1/((1 -θ u )ᾱ) where ᾱ is the average of the α g 's across g, and θ u ∈ (0, 1) quantifies how much the groups of individuals are distinct. Hence, each groups of individuals U g corresponds to a block B U,g . Thus, this generation pattern requires that K > N . In practice, we fix α g ∈ {100, 250}, we use θ = 0.5 or 0.8 (for low or high separation respectively) and N = 3 groups of individuals.

Generation of V. The question of simulating data based on a sparse representation V of the variables in our context of matrix factorization is not straightforward. Indeed, if we impose that some variables j do not contribute to any component k, i.e. that V jk is null for any k, then k U ik V jk is always null for i = 1, . . . , n. Thus, the recorded data entry X ij will be deterministic and null for any observation i (i.e. the j th column in X will be null). There is no interest to generate full columns of null values in the matrix X, since it is unnecessary to use a statistical analysis to determine that a column of zeros will not be informative. This question is not an issue about the formulation of the model, but rather concerns the generation of non informative columns in X that will correspond to null rows in the matrix V.