Normalization and gene selection for single-cell RNA-seq UMI data using sampling-adjusted sums of squares of Pearson residuals with a Poisson model

SCTransform in Seurat and scanpy.experimental.pp.recipe pearson residuals (scanpy henceforth) normalize UMI counts as Pearson residuals of negative binomial models. Residual variance scores genes for downstream analysis. Although we observed that both methods usually assign the highest scores to the same genes, for many highly ranked genes (e.g. among the top 2,000) scores may be unstable - not robust to the selection of cells used to calculate residuals. As an alternative, we consider the Poisson model, for which a natural score is the mean sum of squares of Pearson residuals. We show that these scores can be unstable if a genes nonzero UMI counts are concentrated on a small number of cells. This explains the instability for scanpy because of its similarity to the Poisson model. We define a metric for genes instability and observe that for all three methods it is negatively correlated with the number of cells on which genes counts are nonzero. To reduce the instability of scores based on the Poisson model, we score each gene using multiple random samples of approximately half of the cells. The minimum of these values defines a "sampling-adjusted" score. For data that we analyzed, these are more stable than scores from SCTransform and scanpy while generally agreeing with them on the highest ranked genes. As a second criterion to compare our proposal with SCTransform, we use differential expression analysis. For genes with high scores, the residuals Kruskal-Wallis H-statistics are generally greater for our method than for SCTransform and are more highly correlated with our methods scores.