Depower: An R Package for Simulation-Based Power Analysis of Differential Expression Studies

Abstract

Sample size calculations and power analyses are essential components of experimental design in modern biomedical research. Designs that account for sample correlation, multiple testing, and other sources of variability inherent to specific studies are routinely employed for identifying differential expressions. Despite recent advances in methodologies and software tools for power analysis, there remains a lack of statistical packages capable of accommodating these complex designs in differential expression studies. To fill this gap, we provide the R package depower, which implements the simulation-based framework presented in our recent publications. This unified framework covers both independent and dependent group comparisons and controls false positive rates by employing a simulation-based approach to calculate the empirical null distribution of test statistics.

Introduction

Using biological samples from different conditions to uncover changes in gene expression profiles has become the primary objective of many transcriptomic studies. Sample size calculations and power analyses are essential components of experimental design for these differential expression analyses. To address the need, many methods have been proposed specifically for commonplace RNA-seq experiments with count data. Most of these methods are based on Poisson or negative binomial distribution. However, these approaches often rely on asymptotic approximations, which can lead to inflated false positive rates, as demonstrated in prior studies.^1,2 To address this issue, Yu et al. ³ introduced a simulation-based framework that leverages the empirical null distribution to establish cutoff values for differential expression detection, thereby circumventing the reliance on asymptotic approximations. The simulation-based framework was originally intended for designs utilizing uncorrelated samples. Yu et al. ⁴ expanded its functionality to accommodate designs featuring diverse correlation structures. Other simulation-based approaches for power analysis of bulk and single-cell RNA-seq experiments, such as PROPER ⁵ and POWSC, ⁶ are available in literature. However, these approaches are also affected by inflated false positive rates due to their reliance on asymptotic tests. In addition, they need to simulate genome-wide expression data, whereas our proposed approach simulates expression data at gene-level. This gene-level simulation framework requires substantially fewer parameters to be specified for the simulations as a priori, thereby improving directness and robustness.

In this note, we present the R package depower, which is an implementation of the power analysis procedures developed by Yu et al.^3,4 Although originally developed for RNA-seq data, the methodologies and the package are equally applicable to other high-throughput technologies with count or lognormal data, such as ChIP-Seq and proteomics. depower provides a unified framework for power analyses with functions for simulating data, hypothesis testing, calculating power, and visualizing power for a variety of experimental designs.

Methods

depower supports 2 common outcome models: negative binomial (NB) for counts and lognormal for continuous intensities. It accommodates both independent and paired group comparisons. Throughout, the target fold ratio is a ratio of group means on the original scale (geometric means).

Negative Binomial Outcomes

Independent Groups

Gene count approximations $X$ is assumed to follow a negative binomial distribution NB $(\mu ,\theta )$ with mean $\mu$ and dispersion $\theta$ . The variance is given by $\mu +\frac{{\mu }^2}{\theta }$ , based on the indirect quadratic parameterization of the negative binomial distribution. Love et al. ⁷ assume a dependence relationship for the mean and the dispersion over all genes by assuming $\theta$ follows a lognormal distribution and its mean is a function of the mean of normalized read counts $\overline{\mu }$ :

$${\theta }_{tr}(\overline{\mu })=\frac{{\theta }_1}{\overline{\mu }}+{\theta }_0.$$

depower provides likelihood ratio test and Wald test for the ratio of group means. Both tests allow unequal dispersions across groups.

Paired Groups

For paired designs, the bivariate negative binomial (BNB) distribution is derived by compounding two conditionally independent Poisson variables $X_1|G=g~Poisson(\mu g)$ and $X_2|G=g~Poisson(\gamma \mu g)$ with a shared gamma-distributed random effect $G~\Gamma (\theta ,{\theta }^{-1})$ . $\gamma$ denotes a gene’s fold ratio between 2 conditions. To test hypotheses $H_0:\gamma ={\gamma }_0$ versus $H_1:\gamma \ne {\gamma }_0,$ likelihood ratio test and Wald test are available in the package.

Lognormal Outcomes

Independent Groups

Expression $X$ is assumed to follow a lognormal distribution. Thus, $Y=\log \left(X\right)$ has a normal distribution and the difference of arithmetic mean expressions for the log-transformed data can be back-transformed to obtain the ratio of geometric means of the lognormal data. Data is simulated from a normal distribution based on the parameters of the lognormal outcomes. To specify the variability of the lognormal data using its coefficient of variation (CV), we use the relationships $CV\left(X\right)=\sqrt{e^{{\sigma }_Y^2}-1}$ and ${\sigma }_Y=\sqrt{\log (CV{(X)}^2+1)}.$

Unequal CVs under different conditions are allowed for independent groups. Welch’s t-test is provided for hypothesis testing of the log-transformed data.

Paired Groups

Expression data $X_1$ and $X_2$ are assumed to follow a multivariate lognormal distribution, thus $Y_1=\log \left(X_1\right)$ and $Y_2=\log \left(X_2\right)$ follow a multivariate normal distribution. The relationship between correlation and CV is defined by $cor(Y_1,Y_2)=\frac{\log (cor\left(X_1,X_2\right)CV\left(X_1\right)CV\left(X_2\right)+1)}{{\sigma }_{Y_1}{\sigma }_{Y_2}}.$ Paired data may be simulated as 2 dependent samples or as 1 sample of differences. The 1-sample t-test is provided for hypothesis testing of the log-transformed data.

Type I Error Control

Inferences based on likelihood ratio test or Wald test usually depend on asymptotic theory, which approximates their distributions as Chi-square for large sample sizes. However, this testing strategy may result in much smaller critical values than expected.^1,2 To address this issue, we implemented the simulation-based approach by Rettiganti and Nagaraja ² for proper false positive rate control in Yu et al.^3,4 For instance, the empirical null distribution of the test statistics (likelihood ratio test or Wald test) is obtained from simulated experimental data under the null hypothesis for a large number of iterations (eg, 100 000). Then the 100 $(1-\alpha )$ th percentile from the null test statistics is used as a significance cutoff value for testing under the alternative hypothesis by comparing the test statistics to this cutoff value. As demonstrated in Yu et al,^3,4 asymptotic tests substantially inflate false positive rates, whereas the simulation-based approach properly controls them.

Simulation Procedure for Power Analysis

1. Specify all input parameters: sample size per condition $n,$ mean expression $\mu ,$ dispersion $\theta$ , fold ratio $\gamma$ between conditions, nominal false positive rate $\alpha$ , and number of simulations $T$ .

2. Simulate count data T times from NB $(\mu ,\theta )$ or BNB $(\mu ,\gamma ,\theta )$ under both the null and alternative hypotheses using the input parameters listed in Step 1.

3. Fit the NB model or BNB model and obtain test statistics (likelihood ratio test or Wald test) under the null hypothesis for each simulation run.

4. Calculate the 100 $(1-\alpha )$ th percentile of the empirical null distribution of test statistics as the critical value.

5. Fit the NB model or BNB model and obtain test statistics under the alternative hypothesis for each simulation run.

6. Calculate power (percent of rejections under the alternative hypothesis) for the input parameters listed in Step 1.

Example

Consider an RNA-Seq study in which we are interested in the difference in expression of 2 independent groups. Often, thousands of genes are examined in RNA-Seq studies, however, power is often calculated using the average parameter values for a single representative gene of interest. Yu et al.^3,4 demonstrated that public resources such as TCGA and GEO are well suited for use as pilot datasets when planning new gene expression studies under similar experimental conditions (eg, tumor stages or normal tissues). Parameter estimates derived from these pilot datasets, such as gene-level mean expression and dispersion, can be incorporated into power calculations across a range of candidate sample sizes. In this example, using parameter estimates derived from the TCGA breast cancer dataset, we simulate data for 2 independent groups of negative binomial outcomes with a control group mean of 10, a dispersion parameter of 1 for both groups, and a minimum relevant fold change of 2. Suppose the total number of genes is 10 000, the proportion of truly non-differentially expressed genes is 0.8, and the number of acceptable type I errors is set to 50. Using the per-family error rate method, the type I error rate will be set to $\alpha =\frac{50}{10000(0.8)}$ . We use the approximate parametric likelihood ratio test to calculate expected power of detecting differential expression across varying sample sizes for both groups. Figure 1 illustrates a positive relationship between power and sample size in a 2-group experiment, based on user-specified parameters: mean expression of 10, dispersion of 1, and fold ratio of 2 between groups. According to this power analysis, a new study would require at least 60 samples per group to reach 80% power to detect a twofold differential expression at the specified gene-level parameters.

Figure 1.

Example: power versus sample sizes in a 2-group experiment with user-specified parameters: mean expression of 10, dispersion of 1, and fold ratio of 2 between groups.

library(depower)
set.seed(20251004)
sim_nb(
 n1 =c(3, 10, 20, 40, 60, 100),
 mean1 = 10,
 ratio = 2,
 dispersion1 = 1,
 dispersion2 = 1,
 nsims = 10000
) |>
 power (
  ”Simulated NB LRT” = lrt_nb(
   distribution = simulated(nsims = 20000)
 ),
 alpha = 50 / (10000 * 0.8)
) |>
plot(hline = 0.8)

Conclusion

The depower package addresses a critical gap in the design and analysis of differential expression studies by providing a flexible simulation-based framework for power analysis. Unlike traditional approaches that rely on asymptotic approximations, depower leverages empirical null distributions to ensure accurate control of type I error rates, even under dependent correlation structures and small sample sizes. This capability is particularly important for high-throughput experiments, where thousands of hypotheses are tested simultaneously and conventional methods may lead to inflated false positive rates.

Our proposed simulation-based framework is primarily designed for gene-level power analysis, but it can be easily extended to genome-wide power analysis by specifying a sequence of mean expression and dispersion levels representing different genes. For the proposed framework, the computation time is variable and depends mainly on 3 factors: the sample size of the simulated studies, the number of datasets simulated under the alternative hypotheses, and optionally, the number of datasets simulated under the null hypotheses. One limitation of the simulation-based approach is that the simulated data may become skewed under extreme parameter values (eg, mean expression, dispersion, or fold ratio), which can lead to larger variability of power estimates, particularly when sample sizes are relatively small.

By supporting both negative binomial and lognormal outcome distributions, as well as independent and paired designs, depower accommodates a wide range of experimental settings, including RNA-Seq, ChIP-Seq, and proteomics studies. The package offers integrated tools for data simulation, hypothesis testing, power estimation, and visualization, enabling researchers to make informed decisions about sample size and study design before data collection. Future extensions may include support for additional correlation structures, hierarchical models, and more complex experimental designs, further broadening the applicability of this framework. Overall, depower provides a robust and practical solution for researchers seeking reliable power calculations in complex high-dimensional biological studies.