Identifying Differentially Expressed Genes of Zero Inflated Single Cell RNA Sequencing Data Using Mixed Model Score Tests

Abstract

Single cell RNA sequencing (scRNA-seq) allows quantitative measurement and comparison of gene expression at the resolution of single cells. Ignoring the batch effects and zero inflation of scRNA-seq data, many proposed differentially expressed (DE) methods might generate bias. We propose a method, single cell mixed model score tests (scMMSTs), to efficiently identify DE genes of scRNA-seq data with batch effects using the generalized linear mixed model (GLMM). scMMSTs treat the batch effect as a random effect. For zero inflation, scMMSTs use a weighting strategy to calculate observational weights for counts independently under zero-inflated and zero-truncated distributions. Counts data with calculated weights were subsequently analyzed using weighted GLMMs. The theoretical null distributions of the score statistics were constructed by mixed Chi-square distributions. Intensive simulations and two real datasets were used to compare edgeR-zinbwave, DESeq2-zinbwave, and scMMSTs. Our study demonstrates that scMMSTs, as supplement to standard methods, are advantageous to define DE genes of zero-inflated scRNA-seq data with batch effects.

Keywords: differential expression analyses; generalized linear mixed model; observational weights; score test; single cell RNA sequencing; zero inflation.

FIGURE 1

False positive rate control on simulated null Usoskin datasets and Tung datasets. (A) Boxplot of PCER for 30 simulated null Usoskin datasets generated by splatter for each of 12 DE methods. scMMSTs are marked in blue. (B) Histogram of uncorrected p-values for one dataset in panel A. (C) Boxplot of PCER for 30 simulated null Tung datasets generated by splatter for each of 12 DE methods. scMMSTs are marked in blue. (D) Histogram of uncorrected p-values for one dataset in panel C. PCER, per-comparison error rate; DE, differential expression; scMMST, single cell mixed model score test.

FIGURE 2

FDP-TPR curves of DE methods on simulated Usoskin datasets and Tung datasets. (A) Line plot of the FDP-TPR curves for simulated Usoskin datasets generated by splatter for each of 12 DE methods. (B) Line plot of the FDP-TPR curves for simulated Tang datasets generated by splatter for each of 12 DE methods. Circles represent values at a 0.05 nominal FDR threshold and are filled in if the FDP (i.e., empirical FDR) is less than 0.05. DE, differential expression; TPR, true positive rate; FDP, false discovery proportion; FDR, false discovery rate.

FIGURE 3

FDP-TPR curves of DE methods on simulated datasets generated by GLMMs with μ_π = 0. (A) Line plot of the FDP-TPR curves for simulated datasets based on NB GLMMs for each of 12 DE methods with the dispersion parameter θ = 0.5. (B) Line plot of the FDP-TPR curves for simulated datasets based on negative binomial (NB) GLMMs for each of 12 DE methods with θ = 1. (C) Line plot of the FDP-TPR curves for simulated datasets based on NB GLMMs for each of 12 DE methods with θ = 2. (D) Line plot of the FDP-TPR curves for simulated datasets based on Poisson GLMMs for each of 12 DE methods with ${\mathrm{\beta }}_0={\mathrm{\sigma }}_{\mathrm{\beta }}^2=$ 0.01. Circles represent values at a 0.05 nominal FDR threshold and are filled in if the FDP (i.e., empirical FDR) is less than 0.05. DE, differential expression; GLMM, generalized linear mixed model; NB, negative binomial; TPR, true positive rate; FDP, false discovery proportion; FDR, false discovery rate.

FIGURE 4

AUCs of DE methods for simulated datasets generated by GLMMs with μ_π = 0. Adjusted p-values are used as predictors. (A) Bar plot of AUCs for simulated datasets generated by NB GLMMs for each of 12 DE methods with the dispersion parameter θ = 0.5. (B) Bar plot of AUCs for simulated datasets generated by NB GLMMs for each of 12 DE methods with θ = 1. (C) Bar plot of AUCs for simulated datasets generated by NB GLMMs for each of 12 DE methods with θ = 2. (D) Bar plot of AUCs for simulated datasets generated by Poisson GLMMs for each of 12 DE methods. AUC, area under curve; DE, differential expression; GLMM, generalized linear mixed model; NB, negative binomial.

FIGURE 5

Computational times for differential expression methods on the simulated null Usoskin and Tung datasets, which were generated by splatter. The number of cores were set to be 1 and 8 on a cluster with 24 Intel Xeon Processor (Skylake, IBRS) at 2.60 GHz (2593 MHz) and 128 GB RAM.