How to Use edgeR, DESeq2, and t-test in RNA-Seq Differential Expression Analysis - Practical Rules by Data Type

  • Gene Expression
  • High-Throughput Sequencing

In RNA-Seq differential expression analysis, statistical methods designed for Gene Counts, such as edgeR and DESeq2, are commonly recommended. These methods are designed to take the characteristics of RNA-Seq Counts into account, and are often described as methods that can provide reliable analysis even with a small number of samples.

However, since the early days of microarray analysis, the omics field has had an unfortunate tradition of blindly following “the method everyone uses.” The efforts made to obtain more reliable results from small sample sizes are highly valuable, but there is still no magic cure that works for every dataset.

In real RNA-Seq data, the nature of the data differs greatly depending on many factors: whether the data come from an in vitro experiment or biopsy samples, whether the number of replicates is large or small, whether the input RNA amount is high or low, and whether the study uses a paired design. Naturally, when the nature of the data differs, a method that works well in one situation may not work well in another.

In this article, we first review the basic differences among edgeR, DESeq2, and the t-test. Then, by looking at four representative experimental data types, we examine what kinds of genes tend to be judged significant by each method, and consider how these methods can be used in practice.

Basic Differences Between the t-test and edgeR/DESeq2

A major difference between the t-test and edgeR/DESeq2 lies in how they handle variance.

However, there are three gaps between how edgeR and DESeq2 are commonly explained, the impressions created by those explanations, and what actually happens in the calculations and results.

The t-test evaluates whether the difference between groups is large relative to the variation actually observed among samples for each gene. In other words, the within-group variance observed for each gene is used directly in the test for that gene.

In contrast, edgeR and DESeq2 do not simply use the variance observed for each gene as it is. They use the variation observed across many genes with similar Gene Counts to estimate the relationship between the mean and variance across the RNA-Seq dataset, and then move the variance estimate for each gene toward that overall trend. This process is generally referred to as dispersion shrinkage.

In common explanations of RNA-Seq differential expression analysis, the use of a negative binomial distribution to model Gene Counts is often emphasized as a defining feature of edgeR and DESeq2. Certainly, in the low-Count range, random detection effects become more pronounced, so that a gene may be detected in one sample but not in another, and the variation tends to be large relative to the mean. Using a model suited to Count data is therefore meaningful.

However, in standard analysis workflows, many of these low-expression genes are removed by filtering. Therefore, after filtering, dispersion shrinkage performed by edgeR and DESeq2 may have a greater influence on the results than the use of the negative binomial distribution itself.

Although common explanations emphasize the negative binomial distribution, the influence of dispersion shrinkage, which can substantially affect the actual calculations, is not explained nearly as clearly. This is the first gap between common explanations and what happens in practice.

Under dispersion shrinkage, when the observed variance of an individual gene is larger than that of other genes with similar Gene Counts, the variance used in the test may be adjusted downward. As a result, some genes may be more likely to be identified as significant by edgeR or DESeq2 than by a t-test, which directly uses the variance actually observed for that gene.

This differs from the common impression that edgeR and DESeq2 are more conservative than the t-test. For genes with a small observed variance, the variance may be adjusted upward, making them less likely to be identified as significant. On the other hand, for genes with a large observed variance, the variance may be adjusted downward, making them more likely to be identified as significant than with a t-test.

Therefore, edgeR and DESeq2 are not necessarily more conservative than the t-test for every individual gene. This is the second gap.

Furthermore, dispersion shrinkage is not used only when the number of replicates is extremely small and the variance of each gene cannot be estimated reliably. It is still applied even when there is a moderate number of samples and the observed variance of each gene might be considered reasonably reliable.

As the number of samples increases, greater weight is given to the variance observed for each individual gene. However, the dispersion trend estimated from other genes with similar Gene Counts continues to influence the test results.

For this reason, simply describing this process as “stabilizing variance estimates by sharing information across genes” does not fully explain what is actually being done or how the results may be affected.

Even when the large variance observed for an individual gene reflects genuine biological heterogeneity, that variance may still be adjusted downward. In such cases, it is worth considering whether this process should simply be described as “stabilizing variance estimates.” This is the third gap.

However, apart from the characteristics of dispersion shrinkage discussed above, edgeR and DESeq2 also brought an important improvement to RNA-Seq differential expression analysis. They established the practice of evaluating differential expression using Gene Counts rather than TPM or FPKM. I consider this one of the most important contributions of edgeR and DESeq2.

By normalizing for gene length, TPM and FPKM alter the relationship between the mean and variance present in Gene Counts. This effectively discards much of the theory and practical knowledge accumulated through the analysis of experimental measurement data, not only in RNA-Seq, and has contributed to unreliable DEG analysis. The differences among TPM, FPKM, and Gene Counts are explained in more detail in this article.

With these characteristics in mind, we will now examine four representative data types to see how the results obtained with edgeR, DESeq2, and the t-test differ, and how these methods can be used appropriately.

Case 1: Low-Variance, Small-Sample in vitro Data

This case is examined in detail in the case study comparing edgeR, DESeq2, and the t-test using low-variance, small-sample in vitro data.

In in vitro experiments, there is no variation caused by individual biological differences, so within-group variance can easily become very small. In addition, when the number of replicates is extremely small, such as two replicates per group, the variance of a gene can be underestimated simply because the two observed values happen to be close.

Of course, such an experimental design is far from what statistical models ideally assume. However, it is also a condition that is often seen in real RNA-Seq data. For this reason, it is important to understand what kinds of genes tend to be judged significant by each method under such conditions.

In the actual analysis, because the observed within-group variance was extremely small, all three methods judged even genes with very small changes as significant.

Therefore, in this type of case, it is recommended not to extract DEGs based on p-values alone, but to combine the p-value condition with a fold-change condition. This removes genes that were judged significant despite only very small changes. When a fold-change condition was added, the differences among the DEG lists obtained by the three methods became much smaller.

However, even after adding a fold-change condition, some differences among the methods remain.

Genes significant only by the t-test tend to have extremely small observed within-group variance . For such genes, the p-value in the t-test can become small simply because the replicate values happened to be similar. However, this type of error can be handled relatively easily by combining the p-value condition with a fold-change condition. The same approach is also useful when using edgeR or DESeq2.

Next, genes judged significant by edgeR or DESeq2 tend to include genes with relatively large observed within-group variance. This is because variance is estimated not only from the observed variance of the gene itself, but also from the variation of other genes with similar Count levels. In low-variance in vitro data, many genes have very small observed variation, and as a result, some genes may be judged significant more easily by edgeR or DESeq2 than by the t-test.

Case 2: Medium-Sized Biopsy-Derived Data

This case is examined in detail in the case study comparing edgeR, DESeq2, and the t-test using medium-sized biopsy-derived data.

Next, let us consider medium-sized RNA-Seq data derived from biopsy samples. Compared with Case 1, this type of data has a larger number of samples, and in signal regions where expression levels are sufficiently high, variance estimation becomes relatively stable. Therefore, in regions where expression levels are sufficiently high, the differences among edgeR, DESeq2, and the t-test become small, and it is reasonable to expect that the major differentially expressed genes can be detected by any of these methods.

On the other hand, in this case, many genes in the low-expression region were judged significant only by edgeR or DESeq2. When these genes were examined in detail, many of them appeared to correspond not to clear biological expression differences, but rather to shifts in the lower end of normalized Gene Counts caused by differences in dynamic range among samples.

In RNA-Seq data, it is not unusual for total read counts to differ by several-fold among samples. Differences in total read counts and differences in dynamic range are not simply proportional, but differences in dynamic range can cause shifts in the lower end of normalized Gene Counts. This is a problem that is commonly seen in RNA-Seq experimental data of a certain scale. When this shift in the lower end is biased toward one group, it can appear as a group difference.

The statistical models of edgeR and DESeq2 are designed to account for the instability of Counts in low-expression regions. However, shifts in the lower end of normalized Gene Counts caused by differences in dynamic range among samples are not explicitly incorporated into these models. Therefore, for genes in the low-count region that are significant only by edgeR or DESeq2, it is risky to judge them based only on p-values.

In standard edgeR and DESeq2 workflows, low-expression genes are usually filtered before testing, so this problem tends to become apparent only when differences in dynamic range are extremely large. Even so, this property of RNA-Seq data should always be kept in mind. Otherwise, there is a risk of accepting DEG candidates that are actually caused by differences in dynamic range.

As with edgeR and DESeq2, appropriate preprocessing and filtering are also essential when using the t-test. When low-count measurements were handled properly , the differences among the three methods became smaller. In other words, for reliable DEG analysis, it is important not only to choose a statistical method, but also to perform appropriate preprocessing and filtering before statistical testing.

However, this approach can also cause genes showing ON/OFF-like behavior in the low-count region to be missed. For such genes, it is safer to extract them separately based on the presence or absence of measured values, rather than relying on p-values that are affected by low-count instability and differences in dynamic range.

In contrast to Case 1, in medium-sized data with a certain number of samples, it is not always necessary to apply a fold-change condition uniformly. In signal regions where expression levels are sufficiently high, genes with small but consistent changes may reflect biologically meaningful subtle differences. However, low-change DEGs should still be interpreted more carefully than genes showing large expression changes.

Case 3: Paired Design Data

This case is examined in detail in the case study comparing edgeR, DESeq2, and the paired t-test using paired design data.

The third case is RNA-Seq data with a paired design. For example, when control and case samples are obtained from the same patient, the paired t-test evaluates the direction and magnitude of the change from control to case within each pair, after canceling out individual differences. For this reason, the paired t-test generally has higher detection power than an unpaired two-group comparison.

edgeR and DESeq2 can also account for individual differences by including patient information or other pairing information in the design. However, edgeR and DESeq2 do not simply evaluate the direction of change within each pair. They also consider the mean-variance relationship of Gene Counts, the uncertainty of low Counts, and variance estimation across the whole model.

Therefore, even when paired design is accounted for in edgeR or DESeq2, the number of significant genes does not necessarily increase. In some cases, the number of significant genes may even decrease compared with an unpaired two-group comparison.

Rather than arguing whether the paired t-test or edgeR / DESeq2 is correct, it is more important in practice to understand that these methods can produce such different results.

Case 4: Low-Input RNA-Seq Data

This case is examined in detail in the case study comparing edgeR, DESeq2, and the t-test using low-input RNA-Seq data.

The fourth case is low-input RNA-Seq. In this case, before choosing a statistical method, it is necessary to check whether the data themselves have enough reproducibility to support DEG analysis.

With extremely low input amounts, such as 100 cells, there was almost no reproducibility among replicate samples. Many of the factors that appeared as large expression differences were caused by stochastic amplification or stochastic detection failure.

In this type of data, the structure of the data differs greatly from the variance structure assumed for ordinary bulk RNA-Seq by the statistical models of edgeR and DESeq2. Therefore, even if genes are judged significant by edgeR or DESeq2, many of them may reflect apparent differences caused by amplification or detection rather than true biological expression differences.

This does not mean that edgeR or DESeq2 is wrong. Rather, it is a natural consequence of the input data structure deviating greatly from the variance model assumed for ordinary bulk RNA-Seq.

At 1000 cells, a slightly more RNA-Seq-like variance structure begins to appear on the scatter plot. In other words, high-expression regions show smaller variation, while variation spreads toward the low-expression region, approaching the mean-variance relationship assumed by edgeR and DESeq2.

However, the fact that this variance structure begins to appear does not mean that DEG analysis results from edgeR or DESeq2 can be trusted. Even at 1000 cells, the data remain unstable compared with ordinary bulk RNA-Seq, and stochastic amplification and detection still remain.

In this dataset, the 100k-cell result could be used as a reference. When the 1000-cell result was compared with the 100k-cell result, consistent changes were confirmed only for very highly expressed genes with Gene Counts above 5000.

In other words, even if 1000-cell data appear to approach the variance structure assumed by edgeR and DESeq2, the resulting DEG list should not be trusted without caution. DEG lists obtained from low-input RNA-Seq data need to be handled extremely carefully, for example by focusing only on very highly expressed genes.

Practical Rules by Data Type

Data type Main issue Points to note when using edgeR / DESeq2 Points to note when using the t-test Practical rule
Low-variance, small-sample in vitro data Within-group variance can appear extremely small. Because n is small, gene-wise variance tends to be underestimated. Because variance is estimated using genes with similar Count levels, genes with large observed variance can still be judged significant. In such cases, significance calls may be more lenient than those from the t-test. Genes with extremely small observed within-group variance are more likely to be judged significant. Use a fold-change condition together with the p-value. Adding a fold-change condition reduces the differences among the three methods.
When using edgeR / DESeq2, consider excluding genes with large observed variation from the DEG list.
Medium-sized biopsy-derived data Because there is a certain number of samples, variance estimation becomes more stable. However, individual differences, tissue composition, total read counts, dynamic range, and other factors can still make DEG reliability unstable. In the boundary region near noise, apparent differences caused by dynamic range differences may be judged significant.
Genes with large observed variance can also be judged significant.
In the boundary region near noise, apparent differences caused by dynamic range differences may be judged significant.
Genes with extremely small observed within-group variance are more likely to be judged significant.
Check the distribution of normalized Gene Counts and set the low signal cutoff according to the sample group with the narrowest dynamic range. Perform ordinary DEG analysis in the signal region, and separately extract ON/OFF-type genes in the low-count region based on the presence or absence of measured values, rather than relying on p-values.
Paired design If individual differences are not accounted for, baseline differences among patients are mixed into the evaluation of group differences. Patient information can be included in the design to account for individual differences, but because the overall uncertainty of the Count model is also evaluated, the number of significant genes does not necessarily increase. Because within-pair differences are evaluated directly, genes with small but consistent changes within patients are more easily detected. However, these genes cannot simply be dismissed as biologically meaningless noise. Always account for the paired structure. Genes significant only by the paired t-test should also be examined as genes showing small but consistent within-patient changes, using clustering and functional analysis.
Low-input RNA-Seq Stochastic amplification and detection have a large effect, and the data can easily deviate from the variance structure of ordinary bulk RNA-Seq. Genes judged significant may reflect apparent differences caused by amplification or detection rather than true biological differences. The instability of observed values is reflected directly, so even large p-value or fold-change signals may not be reliable. Do not trust DEG lists based on p-values alone. Interpret the results even more conservatively than they may appear from the plots.

Detailed Case Studies

Summary: Choosing Methods Based on the Data and Research Objective

The t-test is not an old or immature method that should simply be replaced by edgeR or DESeq2. More precisely, the t-test places greater emphasis on the variance actually observed for each gene, whereas edgeR and DESeq2 place greater emphasis on stabilizing variance estimates by sharing information across genes. Which approach is more appropriate depends on the characteristics of the data and the objective of the analysis.

When the sample size is small, many low-Count genes are present, and uncertainty in gene-specific variance estimates has a major influence on the results, dispersion shrinkage is a reasonable approach. However, as shown in Cases 1 and 2, even when the methods produce noticeably different results when compared using p-values alone, the final DEG lists become much more similar once a fold change criterion is also applied.

On the other hand, when there is a sufficient number of samples, the data quality is high, low-Count genes have been appropriately removed, and effect size and gene-specific variation are also considered important information, the benefit of sharing information across genes to stabilize variance estimates becomes smaller, while the disadvantage of moving the variance observed for each gene toward an overall trend becomes relatively greater. Under such conditions, using the variance observed for each gene directly is also well justified.

Therefore, the question we should really ask is not, “Should the t-test be avoided in RNA-Seq analysis?”
The more important question is: under what conditions, and to what extent, is it appropriate to share variance information across genes?
Dispersion shrinkage can be viewed as a trade-off between uncertainty in variance estimation and bias. In other words, it reduces the sampling variability of variance estimates at the cost of moving the variance observed for each gene toward an overall trend estimated from many genes with similar Gene Counts. This trade-off can be advantageous in some situations, but disadvantageous in others, particularly when gene-specific variation is itself considered important information.

In practice, neither method should be selected automatically in every situation. The choice should be based on the sample size, the proportion of low-Count genes, data quality, the observed variance and effect size, and the objective of the analysis. The t-test, edgeR, and DESeq2 are all simply tools that analysts select according to the purpose of the analysis.

edgeR/DESeq2 and the t-test differ in the expression patterns for which they tend to produce overly lenient significance calls. With the t-test, genes with small observed variance are more likely to be judged significant, but this can be addressed by combining the test with a fold change criterion. In contrast, edgeR and DESeq2 may identify genes as significant even when their observed variance is large. In situations where highly variable genes are not considered important, it may be worth considering measures such as removing genes with a high CV from the DEG list. In addition, comparing DEG lists obtained with multiple methods and visualizing both their overlap and the characteristics specific to each method can provide more useful practical information than relying on a single list of p-values.

It is also important to take a broader view and to make decisions by considering not only the statistical test itself, but also the steps that come before and after it. Furthermore, regardless of which method is used, visualizing the results and understanding their characteristics is fundamental to DEG analysis.

Related Topics

edgeR & DESeq2 Analysis with ChatGPT (RNA-Seq Tutorial using R)

Chef Choosing the Right Knife