Gene expression variability and the analysis of large-scale RNA-seq studies with the MDSeq

Abstract

Rapidly decreasing cost of next-generation sequencing has led to the recent availability of large-scale RNA-seq data, that empowers the analysis of gene expression variability, in addition to gene expression means. In this paper, we present the MDSeq, based on the coefficient of dispersion, to provide robust and computationally efficient analysis of both gene expression means and variability on RNA-seq counts. The MDSeq utilizes a novel reparametrization of the negative binomial to provide flexible generalized linear models (GLMs) on both the mean and dispersion. We address challenges of analyzing large-scale RNA-seq data via several new developments to provide a comprehensive toolset that models technical excess zeros, identifies outliers efficiently, and evaluates differential expressions at biologically interesting levels. We evaluated performances of the MDSeq using simulated data when the ground truths are known. Results suggest that the MDSeq often outperforms current methods for the analysis of gene expression mean and variability. Moreover, the MDSeq is applied in two real RNA-seq studies, in which we identified functionally relevant genes and gene pathways. Specifically, the analysis of gene expression variability with the MDSeq on the GTEx human brain tissue data has identified pathways associated with common neurodegenerative disorders when gene expression means were conserved.

Figure 1.

Type I errors in the absence of differential expression variability. There are n samples of cases and controls each and varying proportions of excess zeros s. The MDSeq, Levene's tests, and heteroscedastic regression have well controlled type I errors for moderate to large sample sizes, whereas Bartlett's and MAD tests have highly inflated type I errors. Results are based on 1,000 simulations without additional covariates. Reference lines (in red) are drawn at the 0.05 error rate.

Figure 2.

Powers of detecting differential expression variability. There are n samples of cases and controls each and varying proportions of excess zeros s and log₂ fold-changes log₂FC. The MDSeq often performs the best when sample sizes are moderate or large. Levene's tests and heteroscedastic regression tend to deteriorate in performance with increasing proportions of excess zeros s. Results are based on 1,000 simulations without additional covariates.

Figure 3.

Type I errors in the absence of differential expression means. There are n samples of cases and controls each and varying proportions of excess zeros s. The MDSeq controls type I errors well at moderate to large sample sizes. DESeq2 and edgeR methods may be conservative under the presence of excess zeros s > 0. Results are based on 1,000 simulations without additional covariates. Reference lines (in red) are drawn at the 0.05 error rate.

Figure 4.

Powers of detecting differential expression means. There are n samples of cases and controls each and varying proportions of excess zeros s and log₂ fold-change log₂FC. The MDSeq and ShrinkBayes often perform the best among methods compared under the presence of excess zeros s > 0. Results are based on 1,000 simulations without additional covariates.

Figure 5.

Hypothesis tests to evaluate absolute log fold-changes of expression means above given thresholds. Type I errors are shown at |log₂FC| less than or equal to a given threshold, and powers are presented when |log₂FC| is greater than the given threshold. The MDSeq and DESeq2 have well controlled type I errors, whereas edgeR methods have highly inflated type I errors when |log₂FC| is at or moderately less than the given thresholds. The MDSeq has greater power than DESeq2 when |log₂FC| is moderately above the given thresholds. There are 500 samples of cases and controls each. Results are based on 1,000 simulations generated from NB(2^log2FCμ₀, ϕ₀) with μ₀ = exp (5) and ϕ₀ = exp (4) for varying log₂FC. No excess zeros were generated with s = 0. Reference lines (in gray) are drawn at the corresponding threshold levels.

Figure 6.

Accuracy of the computationally efficient one-step estimator . The computationally efficient one-step estimator is compared with the leave-one-out influence measure I_i under scenarios when (A) there are no outliers and when (B) outliers are present. There are n = 250 samples each for cases and controls. Counts were generated from NB(μ_i, ϕ_i) for controls and NB(2^log2FCμ_i, ϕ_i) for cases with log₂FC = 2, where μ_i = exp (5 + (x_i1 + x_i2)/2) and ϕ_i = exp (4 + (x_i1 + x_i2)/2). Additional covariates x_i1 and x_i2 were simulated from binomial distributions Binom(2, prob = (0.5, 0.5)). In (B), five samples were randomly replaced by outliers simulated from Pois(exp (5)exp (4)). Non-outlying samples (in black) and outliers (in magenta) are plotted. Reference lines (in dashed blue) are drawn at the (α_out/2)th- and (1 − α_out/2)th-quantile of the variance-gamma distribution with α_out = 0.05. A diagonal reference line (in solid red) is drawn at equality of and I_i. In (B), all five outliers were identified by the one-step estimator .