glmGamPoi: fitting Gamma-Poisson generalized linear models on single cell count data

1The statistical distribution of sequencing counts from single-cell RNA-Seq can be modeled by, where y are the observed counts for a particular gene across a set of sufficiently similar cells (acting as replicates) and μ represents the underlying, true expression level of the gene (the expectation value). The parameter determines the dispersion of the distribution; the tightest case is θ = 0, in which case the distribution coincides with the Poisson distribution. Larger values of θ correspond to wider distributions.

2Biological interest is added by extending the model beyond (conceptual) replicates and letting μ vary across the cells. This can be done in different ways: via a generalized linear model,, as in the differential expression methods edgeR and DESeq; via a factor analysis model; or via a matrix decomposition analogous to PCA. The model fits then provide biological insight about ‘significant’ variations in gene expression across cells, above and beyond the sampling noise.

3A popular alternative approach is to transform the counts using the shifted logarithm, with some choice of c > 0, and then proceed with analysis methods that are based on the least squares error, such as used for normal distributed data. However, this approach is fundamentally inferior as it overemphasizes the influence of small count fluctuations and deals poorly with variable sequencing depth across cells.

4With the Gamma-Poisson generalized linear model, parameter estimation proceeds by minimizing the deviance, a generalization of the sum of squares of residuals used in the least squares method. There are already a number of implementations to this end, including the R packages MASS, edgeR and DESeq2. These all follow a similar approach: for each gene, the parameter vector β is estimated using an iterative reweighted least squares algorithm, and the dispersion θ is found by likelihood maximization. After years of development, and with tens of thousands of users, edgeR and DESeq2 provide robust implementations and are a de facto standard for the analysis of bulk RNA-seq data Application of these implementations to single-cell RNA-seq data, however, suffers from several issues. First, their runtime becomes excessive as the number of cells gets large. Second, their functionality—fitting a Gamma-Poisson generalized linear model for a fixed, known design matrix X—is only part of what users need: with single-cell RNA-seq data, important research questions include identification of latent factors, dimension reduction, clustering and classification. These limitations hamper the development and uptake of statistical models based on the Gamma-Poisson distribution and appear to be driving analysts toward the transformation approach.

5The R package glmGamPoi provides inference of Gamma-Poisson generalized linear models (details of the algorithm in Supplementary Appendix S1) with the following improvements over edgeR and DESeq2:

6Substantially higher speed of the overdispersion estimation, by using an efficient data representation that makes uses of the fact that most entries in the count matrix are from a small set of integer numbers.

7Better estimates (i.e. larger likelihood) of the overdispersion on datasets with many small counts.

8No size limitations for the datasets. glmGamPoi supports fitting the model without loading all data into RAM simultaneously (i.e. working ‘on-disk’), by using the HDF5Array and beachmat packages.

9Small number of R package dependencies to facilitate use as a building-block for higher-level methods, such as factor analysis, dimension reduction or clustering/classification.

10Like edgeR, glmGamPoi also provides a quasi-likelihood ratio test with empirical Bayesian shrinkage to identify differentially expressed genes. In addition, it provides the option to form a pseudobulk sample, which found to be an effective way to identify differential expression between samples for replicated single cell experiments.

11To demonstrate how glmGamPoi can be integrated into other tools, we forked the DESeq2 package and integrated glmGamPoi as an alternative inference engine (github.com/mikelove/DESeq2/pull/24).

12We compared the runtime of glmGamPoi to other methods on the four single cell datasets summarized in Supplementary Table S1. The timing results are shown in Figure 1 and Supplementary Figure S1. The speedup by glmGamPoi compared to edgeR and DESeq2 was 6× to 13×. When the data were accessed directly from disk, the calculations took about twice as long. Omitting the estimation of θ (by setting it to zero) reduced the runtime to about a half. The forked version of DESeq2 that calls glmGamPoi was about as fast as calling glmGamPoi directly, indicating that inference carried out by glmGamPoi uses the largest part of the compute resources, while the additional steps carried out by DESeq2 make relatively small demands. Although all methods theoretically scale linearly with the number of genes and cells, we find empirically that glmGamPoi scales sub-linearly with the number of cells, which explains the observed performance benefit (Supplementary Fig. S2).

13On the PBMC68k dataset, the calculations of DESeq2 and edgeR aborted because they ran out of memory (250 GB of RAM available). In contrast, glmGamPoi completed after ca. 45 min (Supplementary Fig. S1).

14Supplementary Figures S4–S6 show that glmGamPoi’s gain in performance does not come at a cost of accuracy. On the contrary, Supplement Figure S3 shows that glmGamPoi provides better estimates (in the sense of larger likelihood) than DESeq2 for 72% of the genes and 10% of the genes in comparison with edgeR. Those differences with edgeR, seem to be of minor importance for assessing differential expression: the bottom rows of Supplementary Figures S5 and S6 show that the P-values from glmGamPoi and edgeR are very similar, consistent with the fact that they use the same statistical test. In Supplementary Figure S7, we provide a more detailed comparison for which genes the P-values of glmGamPoi and DESeq2 are similar and for which genes they are different.

15This work was supported by the EMBL International PhD Programme. In addition, this work has received funding from the European Research Council Synergy Grant DECODE under grant agreement No. 810296.

16Conflict of Interest: none declared.