Confronting false discoveries in single-cell differential expression

Abstract

Differential expression analysis in single-cell transcriptomics enables the dissection of cell-type-specific responses to perturbations such as disease, trauma, or experimental manipulations. While many statistical methods are available to identify differentially expressed genes, the principles that distinguish these methods and their performance remain unclear. Here, we show that the relative performance of these methods is contingent on their ability to account for variation between biological replicates. Methods that ignore this inevitable variation are biased and prone to false discoveries. Indeed, the most widely used methods can discover hundreds of differentially expressed genes in the absence of biological differences. To exemplify these principles, we exposed true and false discoveries of differentially expressed genes in the injured mouse spinal cord.

Conflict of interest statement

G.C. is a founder and shareholder of Onward Medical, a company with no direct relationships with the present work. The remaining authors declare no competing interests.

© 2021. The Author(s).

Figures

Fig. 1. A systematic benchmark of differential expression in single-cell transcriptomics.

a Schematic overview of the eighteen ground-truth datasets analyzed in this study. b Statistical methods for DE analysis employed in 500 recent scRNA-seq papers. Grey bars represent DE analysis methods included in this study. “Other” includes methods used in two or fewer studies. Inset pie chart shows the total proportion of recent scRNA-seq papers that employed DE analysis methods included in this study. Source data are provided as a Source Data file. c Area under the concordance curve (AUCC) for fourteen DE methods in the eighteen ground-truth datasets shown in a. d Mean difference in the AUCC (∆ AUCC) between the fourteen DE methods shown in c. Asterisks indicate comparisons with a two-tailed t-test p-value less than 0.05. e AUCC of GO term enrichment, as evaluated using gene set enrichment analysis, in the eighteen ground-truth datasets shown in a. f Rank and statistical significance of the GO term GO:0043330 (“response to exogenous dsRNA”) in GSEA analyses of mouse bone marrow mononuclear cells stimulated with poly-I:C, a type of synthetic dsRNA, for four h, using the output of fourteen DE methods. Source data are provided as a Source Data file.

Fig. 2. Single-cell DE methods are biased towards highly expressed genes.

a Schematic illustration of the creation of ‘pseudobulks’ from single-cell data. Top, biological replicate from which each cell was obtained. Bottom, simulated gene expression matrix. Read counts for each gene are aggregated across cells of a given type within each biological replicate. b Mean AUCCs across eighteen ground-truth datasets after dividing the transcriptome into terciles of lowly, moderately, or highly expressed genes. c Mean expression levels of the 100 top-ranked false-positive genes from each DE method. d Spearman correlation between the mean expression of 80 ERCC spike-ins expressed in at least three cells and the –log10p-value of differential expression assigned by each DE method. e Scatterplots of mean ERCC expression vs. –log10p-value for exemplary single-cell and pseudobulk DE methods. Trend lines and shaded areas show local polynomial regression and the 95% confidence interval, respectively. f Mean expression levels of the 200 top-ranked genes from each DE method in a collection of 46 scRNA-seq datasets.

Fig. 3. DE analysis of single-cell data must account for biological replicates.

a Schematic illustration of the experiment shown in b, in which the aggregation procedure was disabled and pseudobulk DE methods were applied to individual cells. b Left, AUCC of the original fourteen DE methods, plus six pseudobulk methods applied to individual cells, in the eighteen ground-truth datasets. Right, Spearman correlation between ERCC mean expression and –log10 p-value assigned by six pseudobulk DE methods, before and after disabling the aggregation procedure. c Schematic illustration of the experiment shown in d, in which the replicate associated with each cell was shuffled to produce ‘pseudo-replicates.’ d Left, AUCC of the original fourteen DE methods, plus six pseudobulk methods applied to pseudo-replicates, in the eighteen ground-truth datasets. Right, Spearman correlation between ERCC mean expression and –log10 p-value assigned by six pseudobulk DE methods, before and after shuffling replicates to produce pseudo-replicates. e Variance of gene expression in pseudobulks formed from biological replicates and pseudo-replicates in mouse bone marrow mononuclear cells stimulated with poly-I:C. Shuffling the replicate associated with each cell produced a systematic decrease in the variance of gene expression. Right, pie chart shows the proportion of genes with increased or decreased variance in pseudo-replicates, as compared to biological replicates. f Decreases in the variance of gene expression in pseudo-replicates as compared to biological replicates across 46 scRNA-seq datasets. Points show the mean variance in biological replicates; arrowheads show the mean variance in pseudo-replicates. g Left, expression of the gene c in biological replicates (points) and pseudo-replicates (arrowheads) from unstimulated cells and cells stimulated with poly-I:C, with the range of possible pseudo-replicate expression values shown as a density. Right, mean (horizontal line) and variance (shaded area) of Txnrd3 expression in biological replicates (left) and pseudo-replicates (right). P-values were calculated by edgeR-LRT.

Fig. 4. False discoveries in single-cell DE.

a Schematic illustration of simulation experiments. Single-cell RNA-seq datasets were simulated with varying degrees of heterogeneity between replicates. Replicates were then randomly assigned to either a ‘treatment’ or ‘control’ group, and DE analysis was performed between groups. b Number of DE genes detected in stimulation experiments with varying degrees of heterogeneity between replicates by a representative single-cell DE method, a representative pseudobulk method, and the same pseudobulk method applied to pseudo-replicates. Points and error bars show the mean and standard deviation across ten independent simulations. c Number of DE genes detected by the tests shown in b for genes divided into deciles by the magnitude of the change in variance between biological replicates and pseudo-replicates (∆-variance). d Volcano plots showing DE between T cells from random groups of unstimulated controls drawn from Kang et al. using a representative pseudobulk method, edgeR-LRT, applied to biological replicates or pseudo-replicates. Discarding information about biological replicates leads to the appearance of false discoveries. e Number of DE genes detected in comparisons of random groups of unstimulated controls from 14 scRNA-seq studies with at least six control samples. f Number of DE genes detected within spinal cord regions from control mice profiled by spatial transcriptomics using a representative pseudobulk method, edgeR-LRT (points), or a representative single-cell method, the Wilcoxon rank-sum test (arrowheads). g Mean change in variance between biological replicates and pseudo-replicates for 18 human and 20 mouse scRNA-seq datasets.

Fig. 5. True and false discoveries in the injured mouse spinal cord.

a Synapses and projections of Vglut2ON neurons in the mouse lumbar spinal cord before and after SCI. Experiments were repeated in n = 3 mice. b Top, chronophotography of mice before and after SCI. Bottom left, principal component analysis of gait parameters for each condition (small circles). Large circles show the average per group. Bottom right, bar plot showing the average scores on principal component 1 (PC1), which quantify the locomotor performance of injured and uninjured mice. c snRNA-seq experimental design. Inset demonstrates persistent amplitude at increasing frequency of the H-reflex in response to plantar stimulation, reflecting hyperexcitability. d Uniform manifold approximation and projection (UMAP) visualization of 19,237 nuclei, revealing the major cell types of the mouse lumbar spinal cord. e UMAP visualization colored by the AUC of cell type prioritizations, as calculated by Augur. Endothelial cells are highlighted as the cell type with the highest AUC. f Pseudobulk expression levels in lumbar spinal cord endothelial cells by snRNA-seq (in counts per 10,000), left, and mean counts per Pecam1ON cell by RNA-scope, center, for Igfbp6, a gene identified as DE only by edgeR-LRT. Horizontal line and shaded area show the mean and standard deviation, respectively. Right, colocalization of Igfbp6 with the endothelial marker gene Pecam1 in RNAscope in situ hybridizations from injured and uninjured mouse spinal cords. g As in f, but for Prex2, a gene identified as DE only by the Wilcoxon rank-sum test. h Summary of in vivo screen results. Filled bars represent genes validated by RNAscope. i, Immunohistochemistry reveals Pecam1ON cells with atrophic features after SCI. Experiments were repeated in n = 3 mice. Source data are provided as a Source Data file.