Motivation

Differential expression (DE) analysis in single-cell transcriptomics reveals cell-type-specific responses. Recent studies have raised concerns about current methods used to identify differentially expressed genes in this context.

• scRNA sequencing provides the absolute abundances of RNA molecules in single cells, but normalization - a pre-processing step inherited from the bulk RNA-seq era - reduces this information and returns data as relative abundances, which may mask important differences among cell types

• The majority of single-cell DE analysis methods are susceptible to generating false discoveries. This is mainly due to the lack of accounting for variations between biological replicates, commonly referred to as “donor effects”.

  • Batch effects are often estimated from leading principal components, representing a consensus from most genes.
  • Pseudo-bulk analysis ignores within-sample heterogeneity by treating donor effects as a fixed effect and assuming that each cell from the same donor is equally affected.

• Clustering and DE analysis are different problems. The current commonly workflow works well in clustering, but cannot guarantee success in downstream analysis.

• Excessive zeros are usually considered as “drop-outs”, while they are actually informative in cell-type heterogeneity. Ignoring zeros in single-cell gene expression data discards valuable information for any analysis.

We provided a generalized linear mixed model framework (LEMUR: Linear mixed model-based differential Expression analysis accounting for Multiple samples Using Random effects) to detect differentially expressed genes (DEGs) between two given cell types. The model takes donor-specific variations as random effects, and uses raw UMI counts to prevent biases in DE analysis.

LEMUR framework

By default, we excluded any genes detected in fewer than 5% cells in the compared groups from differential testing. The GLMMs were implemented with glmmPQL function of the MASS package. We calculated adjusted p-values by using Benjamini-Hochberg correction. Each model fitting was applied on one gene and the two compared groups.

Benchmarked methods

In this project, we compared a few methods performing the DE analysis results. Our comparison encompassed Poisson-glmm and Binomial-glmm from the new paradigm, as well as pseudo-bulk approaches including DESeq2 and edgeR. Additionally, we assessed the performance of single cell specific tools including MAST, Wilcox in Seurat, and linear mixed models in Muscat. More details can be found here.

New criteria on determining DEGs

We proposed new criteria that based on the convention and also the gene mean and the difference in mean. If the log2 gene mean in two groups is lower than a certain value (-2.25 in case study 1) and the log2 mean difference is smaller than a threshold (-1 in case study 1), the genes would not be considered as a DEGs. These can also be used as a filter before any DE analysis to speed up the computation. Both of these criteria are adjustable, depending on the dataset’s performance and characteristics. An examination in heatmaps and mean difference against mean plot in advanced can be helpful to determine the thresholds when analyzing a new dataset. More details can be found here.

Application