Source-linked AI summary
Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression
Belinda Phipson, Stanley Lee, Ian J. Majewski, Warren S. Alexander, Gordon K. Smyth
TL;DR
Differential-expression tests can be distorted by genes with unusually large or small variances, especially in small-replicate studies. This paper robustifies empirical-Bayes hyperparameter estimation, increasing power and decreasing false discovery rates when outliers are present.
Problem
Empirical-Bayes differential-expression tests can be unduly influenced by a small number of genes with extreme variances.
Method
The method assigns diffuse variance priors to hypervariable genes while retaining conjugate priors and exact small-sample null distributions.
Results
With outliers present, robust empirical Bayes simultaneously increased power and decreased the false discovery rate while matching the original method without outliers.
Takeaways & Limitations
The procedure protects non-outlier genes from hyperparameter distortion while downweighting hypervariable genes rather than removing them.
Takeaways & Limitations
The approach treats genes, rather than individual expression values, as outliers, so it is complementary to robust regression methods targeting erroneous observations.
Abstract
from arXiv · showhide
One of the most common analysis tasks in genomic research is to identify genes that are differentially expressed (DE) between experimental conditions. Empirical Bayes (EB) statistical tests using moderated genewise variances have been very effective for this purpose, especially when the number of biological replicate samples is small. The EB procedures can however be heavily influenced by a small number of genes with very large or very small variances. This article improves the differential expression tests by robustifying the hyperparameter estimation procedure. The robust procedure has the effect of decreasing the informativeness of the prior distribution for outlier genes while increasing its informativeness for other genes. This effect has the double benefit of reducing the chance that hypervariable genes will be spuriously identified as DE while increasing statistical power for the main body of genes. The robust EB algorithm is fast and numerically stable. The procedure allows exact small-sample null distributions for the test statistics and reduces exactly to the original EB procedure when no outlier genes are present. Simulations show that the robustified tests have similar performance to the original tests in the absence of outlier genes but have greater power and robustness when outliers are present. The article includes case studies for which the robust method correctly identifies and downweights genes associated with hidden covariates and detects more genes likely to be scientifically relevant to the experimental conditions. The new procedure is implemented in the limma software package freely available from the Bioconductor repository.
1. Introduction.
The paper robustifies empirical-Bayes differential-expression testing for high-dimensional gene-expression studies with few replicates by reducing the influence of genes with extreme residual variances. The resulting method preserves exact small-sample inference, is fast and numerically stable, and is implemented in limma.
- Motivation: High-dimensional gene-expression experiments often have few independent replicates, making gene-wise analyses imprecise and motivating information sharing across genes.Empirical Bayes methods borrow information across the dataset to improve genewise variance estimation.
- Contribution: The paper improves limma’s empirical-Bayes differential-expression tests by robustifying hyperparameter estimation for exceptionally large or small residual variances.Genes with extreme variances are treated as outliers, with large-variance outliers termed hypervariable genes.
- Method: The robust approach uses more diffuse conjugate priors for hypervariable genes while retaining exact small-sample null distributions for differential-expression test statistics.The method preserves conjugacy and reduces exactly to the original empirical-Bayes procedure when no outlier genes are present.
- Method: The robust empirical-Bayes procedure is fast, numerically stable, and avoids difficult convergence issues.These properties are presented as advantages of the robust estimation procedure while retaining key features of the original tests.
2. Linear models and moderated t-statistics.
The section models gene expression with weighted linear models and gene-specific variances, then moderates variance estimates using an empirical-Bayes scaled inverse chi-square prior. Under the null, moderated statistics follow exact t or F distributions with increased denominator degrees of freedom.
- Linear models: Weighted linear models represent experimental designs through a full-rank design matrix, with independent observations and known gene-specific weights.Least-squares coefficients are estimated using the diagonal weight matrix W_g.
- Linear models: Gene-wise residual variances use residual degrees of freedom, which are typically n−p but can decrease with missing values or zero weights.The residual variance is computed around the fitted values X β̂_g.
- Moderated variances: Empirical Bayes moderation combines gene-wise variance estimates with a scaled inverse chi-square prior parameterized by d_0 and s_0^2.The resulting moderated variance estimators replace ordinary residual variance estimates in downstream statistics.
- Moderated t-statistics: Under β_gj = 0, the moderated t-statistic follows a t-distribution with d_g+d_0 degrees of freedom, and ordinary genewise t or F-statistics can be moderated similarly.This increases the null denominator degrees of freedom from d_g to d_0+d_g.
3. Robust hyperparameter estimation.
This section replaces limma’s moment-based hyperparameter estimation with a robust procedure that Winsorizes extreme gene-wise variances before estimating the prior parameters. The approach reduces prior influence when variance heterogeneity is driven by outliers and supports identifying genes whose variances are inconsistent with the estimated prior.
- Motivation: As the observed variance of log sample variances increases, the estimated prior degrees of freedom decrease, so less prior information is borrowed for moderated t-statistics.The variance would be reduced if a small number of the most variable genes were excluded, motivating robust estimation.
- Winsorized estimation: The robust procedure applies moment estimation to Winsorized sample variances, resetting specified proportions of the most extreme upper- and lower-tail variances.Typical limits are p_l = 0.05 and p_u = 0.1, while any values strictly between 0 and 0.5 are permissible.
- Winsorized estimation: The Winsorized hyperparameters are estimated by matching the observed log-variance variance to the corresponding theoretical Winsorized F-distribution variance and solving for d0.When gene-wise residual degrees of freedom differ, variances are transformed to equivalent variables with equal degrees of freedom before applying the algorithm.
- Outlier handling: After robust estimation, genes with variances too large to have reasonably arisen from the estimated prior can be identified as outliers.The section raises handling such genes as a separate methodological question, including whether to set d0 = 0 and use ordinary t-tests.
4. Gene-specific prior degrees of freedom.
The method assigns gene-specific prior degrees of freedom by estimating each gene’s probability of being a hypervariable outlier. Most genes retain the common prior degrees of freedom, while clear outliers receive substantially less prior information.
- Gene-specific prior degrees of freedom: The procedure estimates doutlier by matching the largest observed residual variance to the median of the corresponding F distribution.A fast, stable numerical algorithm for finding doutlier is provided later in the paper.
- Gene-specific prior degrees of freedom: Each gene’s posterior probability of being non-outlier is estimated from its outlier-test p-value and its empirical residual-variance tail probability.Assuming most genes are not outliers, the method sets π0 = 1 and estimates the marginal tail probability from the gene’s rank.
- Gene-specific prior degrees of freedom: The estimated outlier probabilities are monotonized by ordering genes by p-value, applying cumulative mean and maximum filters, and approximating isotonic regression.This prevents non-monotonic estimates while allowing a group of cases to have small probabilities rather than requiring only the most extreme case to do so.
- Gene-specific prior degrees of freedom: Most genes retain d0, whereas genes with very small pg values receive much lower prior degrees of freedom.This reduces prior informativeness for clear hypervariable outliers while preserving the common prior for the main body of genes.
5. Covariate dependent priors.
The prior distribution for genewise variances is extended to depend smoothly on each gene’s average log-expression. Robust EB estimation uses a lowess variance trend, detrending before estimation and recombining the trend with the robust prior estimates.
- The variance prior is made gene-specific by allowing it to vary smoothly with each gene’s average log-expression, Ag.
- Robust EB with a variance trend fits a robust lowess curve to log s^2_g versus Ag, applies robust EB to detrended variances, and recombines the estimates with the unlogged trend.
- The robust hyperparameter procedure is implemented in limma’s fitFDist-Robustly, with user-settable Winsorizing tail defaults of p_l = 0.05 and p_u = 0.1.
6. Software implementation.
The robust EB strategy is integrated into edgeR pipelines for microarray and RNA-seq analyses, and glmQLFit also supports csaw and diffHic workflows for ChIP-seq and Hi-C data. Its performance was evaluated through simulations varying hyperparameter settings without differential expression genes or outliers.
- Simulation evaluation: Table 1 reports Type I error rates for standard and robust EB t-tests across p-value cutoffs, simulated d0 values, and datasets without differential-expression genes or outliers.Rates were averaged over all genes in 1,000 simulated datasets.
- Software implementation: The robust EB strategy is implemented in edgeR through estimateDisp and glmQLFit for gene-expression microarray and RNA-seq analysis pipelines.glmQLFit is also used in csaw and diffHic pipelines for ChIP-seq and Hi-C sequencing data.
- Simulation evaluation: Simulations generated expression values for 10,000 genes across six RNA samples in two groups of three, yielding 4 residual degrees of freedom.The hierarchical model used s0 = 0.2 and d0 values of 2, 4, or 10 to evaluate standard and robust hyperparameter estimators.
7. Evaluation using simulated data.
Simulations showed that robust and standard EB tests perform similarly without hypervariable genes, whereas robust tests reduce false discoveries and improve power when 250 outliers are present. The evaluation also examined hyperparameter estimates across datasets with and without outliers.
- 7. Evaluation using simulated data.: With 250 hypervariable genes, robust EB tests consistently produced fewer false discoveries and higher power than standard EB tests.Without hypervariable genes, the two tests were indistinguishable for false discoveries and power.
- 7. Evaluation using simulated data.: The simulations included 500 DE genes with log fold changes drawn from N(0, 4), and results were averaged over 1000 simulations.No genes were simultaneously DE and hypervariable in these simulations.
- 7. Evaluation using simulated data.: Hyperparameter estimates were compared across 1000 simulated datasets with no outliers or 250 hypervariable genes, using true s2 0 = 0.4 and d0 values of 2, 4, or 10.The figure reports standard and robust estimates under both outlier conditions.
8. Case studies.
Case studies show that robust empirical-Bayes estimation detects and downweights variance outliers caused by hidden covariates, improving differential-expression power while preserving biologically relevant findings. Across datasets, hypervariable genes were associated with sex, ribosomal, hemoglobin, and other sample-composition effects rather than the primary condition.
- Suz12 case study: 251 down-regulated and 35 up-regulated probes were identified by both methods, while the robust method uniquely found 22 down-regulated genes versus 16 unique to the standard method.The comparison used EB moderated t-statistics with false-discovery-rate control below 5%; the standard method’s unique genes tended to be sex-linked or hemoglobin-associated.
- Suz12 case study: The robust method identified and downweighted outlier variances associated with sex imbalance, ribosomal RNA retention, and blood or bone-marrow content rather than the Suz12 mutation.Two mutant embryos were female while the remaining mice were male, revealing an unwanted hidden covariate; outliers included sex-linked, ribosomal, and hemoglobin genes.
- HapMap case study: In HapMap Yoruba lymphoblastoid-cell RNA-seq profiles, the analysis compared males with females after voom transformation and examined genewise variability.The dataset contained 69 individuals, and Figure 1c displayed genewise standard deviations.
- Embryonic stem-cell case study: In embryonic-stem-cell RNA-seq data, hyper-variable genes were predominantly associated with ribosomal or hemoglobin genes, suggesting inconsistencies in cell purification and RNA processing.Other hyper-variable genes were located in the major histocompatibility region.
- Cross-study observations: Across routinely checked gene-expression studies, hypervariable genes commonly suggested mechanisms involving sex-linked, ribosomal, mitochondrial, or B-cell-related expression.The authors report that many studies contained subsets of outlier genes whose identities indicated potential sources of variability.
9. Discussion.
The robustified empirical-Bayes procedure protects differential-expression analysis from hyper- and hypovariable genes by fitting hyperparameters to non-outliers and reducing outlier influence. It preserves exact small-sample inference, improves power and false-discovery control when outliers occur, and may generalize beyond this setting.
- 9. Discussion.: Robustified EB protects against hyper- and hypovariable genes by letting non-outliers share information as if outliers were absent while retaining hypervariable genes in the analysis.Hypervariable genes borrow less information from the ensemble rather than being removed.
- 9. Discussion.: Conjugate Bayesian models provide closed-form posterior estimators and exact small-sample p-values; robustness comes from fitting priors to non-outliers and lowering outlier prior degrees of freedom.A practical algorithm assigns prior degrees of freedom to outlier genes, and the DE list is not sensitive to their exact values.
- 9. Discussion.: Default Winsorizing tail proportions work well in practice, but users may increase them when datasets are expected to contain high proportions of outlier genes.The recommendation is presented as an implementation choice rather than a change to the underlying robust EB procedure.
- 9. Discussion.: With outliers present, robustified EB simultaneously increased differential-expression power and decreased the false discovery rate, while matching the original method’s hyperparameter accuracy without outliers.The method therefore preserves the original procedure’s performance in clean settings while improving robustness when outliers occur.
- 9. Discussion.: The method treats genes, rather than individual expression values, as outliers; robust regression is presented as a complementary approach that could be used simultaneously.The limma package has included an M-estimation option for robust regression for over a decade.
- 9. Discussion.: The robust EB strategy may extend to other empirical-Bayes contexts by robustly estimating hyperparameters, testing for outlier cases, and assigning diffuse priors to those cases.The authors suggest this approach may be especially attractive for conjugate Bayesian models and differs from previous robust EB strategies.
10. Appendix: Details of computational algorithms.
The appendix details estimation of EB hyperparameters using logged residual variances, degree-of-freedom adjustments, and numerical computation of log-Winsorized F moments. It also describes abundance-dependent trends and an iterative procedure for determining the outlier degrees of freedom.
- Original hyperparameter estimation: The original non-robust method-of-moments algorithm estimates EB hyperparameters from logged gene-wise residual variances, whose moments are approximately normal and finite.Theoretical means and variances are matched to observed moments to solve for the hyperparameters.
- Abundance-dependent trend: An abundance-dependent trend can be incorporated by fitting a four-degree-of-freedom regression spline to transformed log variances against average log2-expression.The fitted spline values replace the corresponding non-trended transformed mean.
- Unequal degrees of freedom: When genes have unequal residual degrees of freedom, variances are transformed to a common maximum degree of freedom before robust estimation.The non-robust hyperparameters are estimated first, then the gene-wise variances are transformed using the F-distribution.
- Log-Winsorized moments: The mean and variance of the log-Winsorized F distribution are computed by separating tail masses from the central interval and evaluating conditional expectations numerically.After transforming to a uniform interval, Gauss-Legendre quadrature evaluates the required expectations.
- Outlier-degree estimation: The outlier degrees of freedom are found by iterating from the estimated prior degrees of freedom until the maximum variance has an F-distribution right-tail probability of 0.5.The iteration converges monotonically and usually reaches sufficient accuracy after two or three iterations.