Analytic Pearson residuals for normalization of single-cell RNA-seq UMI data

Abstract

Standard preprocessing of single-cell RNA-seq UMI data includes normalization by sequencing depth to remove this technical variability, and nonlinear transformation to stabilize the variance across genes with different expression levels. Instead, two recent papers propose to use statistical count models for these tasks: Hafemeister and Satija (2019) recommend using Pearson residuals from negative binomial regression, while Townes et al. (2019) recommend fitting a generalized PCA model. Here, we investigate the connection between these approaches theoretically and empirically, and compare their effects on downstream processing. We show that the model of Hafemeister and Satija (2019) produces noisy parameter estimates because it is overspecified (which is why the original paper employs post-hoc regularization). When specified more parsimoniously, it has a simple analytic solution equivalent to the rank-one Poisson GLM-PCA of Townes et al. (2019). Further, our analysis indicates that per-gene overdispersion estimates in Hafemeister and Satija (2019) are biased, and that the data analyzed in that paper are in fact consistent with constant overdispersion parameter across genes. We then use negative control data without biological variability to estimate the technical overdispersion of UMI counts, and find that across several different experimental protocols, the data suggest very moderate overdispersion. Finally, we argue that analytic Pearson residuals (or, equivalently, rank-one GLM-PCA or negative binomial regression after regularization) strongly outperform standard preprocessing for identifying biologically variable genes, and capture more biologically meaningful variation when used for dimensionality reduction, compared to other methods.

Comparison of Transformations for Single-Cell RNA-Seq Data

“…Lause et al (2021) argued that neither the estimation of β is nor the estimation of one overdispersion per gene are necessary. Instead, Lause et al (2021) suggested treating the log-size factors as offsets (i.e., fixing β is = 1) and fixing the overdispersion to α = 0.01, because that is roughly the overdispersion they observed in experiments where an RNA solution is homogeneously encapsulated in droplets. Hafemeister and Satija (2020) responded that estimating a gene-wise coefficient for the size factor "allows sctransform to adapt to artifacts and biases" and that fixing the overdispersion to a small value over-emphasizes the variation of highly abundant housekeeping genes.…”

Comparison and evaluation of statistical error models for scRNA-seq

“…We next focused on the application of negative binomial error models, and considered different strategies for parameterizing the level of overdispersion associated with each gene. Recent work (22) suggested that a negative binomial model with a fixed parameterization (for example, inverse overdispersion parameter θ = 100) could be applied to all scRNAseq datasets to achieve effective variance stabilization. To explore whether a single value of θ could be applied to diverse scRNA-seq datasets, we first independently fit θ estimates for each gene in each dataset using a GLM with negative binomial errors (NB GLM), using library size as an offset to account for variation in cellular sequencing depth.…”

“…We also considered the findings from (22), which proposed that θ values should not vary as a function of gene abundance, and suggested that the relationship between these two variables was driven entirely by biases in the parameter estimation procedure, especially when analyzing lowly expressed genes. We first confirmed that lowly expressed genes, particularly those with average abundance < 0.1 UMI/cell, posed difficulties for parameter estimation.…”

“…In Figure 2C-D, as well as associated supplementary figures (S3 -S10) we model each of the scRNA-seq datasets using a NB GLM with a fixed value of θ for each gene in each dataset. To test this, we utilize the 'offset' model as described by Lause et al in (22). We repeated the analysis with three different values for the fixed overdispersion parameter, θ = ∞, θ = 100, and θ = 10.…”

“…Even for methods that assume a NB distribution, different groups propose different methods to parameterize their model. For example, a recent study (22) argued that fixing the NB inverse overdispersion parameter θ to a single value is an appropriate estimate of technical overdispersion for all genes in all scRNAseq datasets, while others (23) propose learning unique parameter values for each gene in each dataset. This lack of consensus is further exemplified by the scvi-tools (11,24) suite, which supports nine different methods for parameterizing error models.…”