1. Support
2. Help - Theory & Case Study
3. Have We Focused Too Much on DEG Analysis in RNA-Seq? | Exploring Other Ways to Read the Data

What the Subtle Changes in 200 HuR Target Genes Suggest About How We Use Omics Data

In the previous case study, "Genes That Are Not Significant in DEG Analysis Can Still Have Experimentally Validated Functions," we examined a study of how HuR suppresses adipocyte differentiation and showed that genes that are not significant in DEG analysis can still have experimentally validated functions.

In this HuR study, candidate genes were not identified using RNA-Seq alone. The authors first used RIP-Seq to identify RNAs bound by HuR and selected 200 genes showing high enrichment in the HuR immunoprecipitation as HuR targets. They then used RNA-Seq to examine how the expression of these genes changed after HuR knockout. INSIG1, which the authors selected for further experimental validation, might not have been chosen as a candidate if they had relied on conventional DEG analysis alone.

Here, we take a closer look at the analysis used in this study.

The 200 genes were not evaluated one by one

In conventional RNA-Seq DEG analysis, Fold Change, P values, FDRs, and related statistics are calculated for each gene, and genes showing larger changes or stronger statistical significance are selected as candidates. The analysis performed on the 200 HuR target candidates in the original paper took a somewhat different approach.

Rather than asking whether each individual gene changed significantly, the authors compared the cumulative distributions of HuR knockout/control Fold Changes between the 200 HuR target candidates and the other genes. Using a Kolmogorov–Smirnov test, they showed that the HuR target group as a whole was shifted toward lower expression in the knockout.

In other words, the question was not:

"Which genes changed significantly?"

but rather:

"Do genes bound by HuR behave differently from other genes as a group?"

Could a small number of genes with large Fold Changes be driving the result?

This analysis raises a reasonable question.

Among the 200 HuR target candidates, most showed no clear expression difference between WT and KO, while only a small number showed relatively large changes. It is therefore reasonable to suspect that the difference between the cumulative distributions might simply be driven by a few genes with large Fold Changes.

To examine this possibility, in this case study we progressively removed HuR target candidates with the largest Fold Changes and compared the target and non-target distributions after each step.

The analysis below is not an exact reproduction of the analysis in the original paper. It is an approximate reanalysis using the RNA-Seq data and HuR target list available for this case study. The original paper treated 200 genes as HuR targets, whereas Fold Changes could be calculated for 199 genes in this reanalysis.

In the non-target group used for comparison, the mean KO/WT Fold Change was approximately 1.090 and the median was approximately 0.995. In the HuR target group, the mean was approximately 0.996 and the median approximately 0.987. The mean in the non-target group was influenced by a small number of large Fold Changes, whereas the medians of both groups were close to 1.

As genes with larger Fold Changes were removed, the median KO/WT ratio in the target group moved from 0.9872 to 0.9910, as expected, approaching 1. In other words, the individual expression changes among the remaining target genes became even smaller.

At the same time, the P value from the KS test comparing the target and non-target distributions decreased from 0.076 to 1.1×10⁻⁴. Because the KS test detects differences in the overall distributions rather than differences in the mean alone, the narrowing of the target distribution after genes with large Fold Changes were removed may itself have contributed to the decrease in P value. Therefore, the lower P value should not be interpreted directly as evidence that the biological difference became stronger.

However, even after the small number of genes with large Fold Changes were removed, the target group still contained slightly more genes with KO/WT below 1, and a distributional difference from the non-target group remained. The pattern observed in the original paper therefore does not appear to be explained simply by a few genes with large Fold Changes.

Does a small change mean there is no information?

In conventional DEG analysis, it is reasonable to examine genes with larger changes or smaller P values first. However, the results of this reanalysis suggest that some information may not be captured by that approach alone.

In the HuR target group, individual Fold Changes were small, and after the genes with the largest Fold Changes were removed, the median of the remaining genes moved even closer to 1. Even so, a distributional difference from the non-target group remained.

As discussed in the previous article, biological function is influenced not only by RNA abundance but also by multiple regulatory layers, including translation efficiency and protein stability.

Even if RNA abundance changes by only about 10% (1.1-fold), small changes in the same direction in translation efficiency or protein stability could combine multiplicatively and produce a larger effect on protein abundance or cellular function. Conversely, a change at one regulatory level may also be compensated for at another.

This suggests that a gene that does not stand out when candidates are ranked by RNA-Seq alone may become meaningful when other types of information are considered together. In this case, one such type of information was direct RNA binding.

There may be other important candidates among the 200 genes

Among the 200 HuR target candidates, INSIG1 did not show an especially large change in the RNA-Seq data. Nevertheless, in the original paper, its role was investigated by combining evidence of direct HuR binding, its known relationship with adipocyte differentiation, and subsequent qPCR, protein, mRNA stability, and functional assays.

If so, there may be other genes among the 200 that do not stand out in RNA-Seq alone but become important candidates when additional experimental information is considered.

In this case study, in addition to the FPKM-based t-test used in the original paper, we analyzed the RNA-Seq data using a t-test on Gene Counts recalculated from the FASTQ files, as well as DESeq2 and edgeR. Low-count genes were also handled according to the respective DESeq2 and edgeR workflows.

The analysis method used in the original paper was only one of several possible analytical paths. Another researcher could have chosen a different method. We therefore applied multiple methods and compared how many of them produced P<0.05 for each gene, as an indicator of whether a gene would tend to emerge as a significant candidate regardless of the analytical path taken.

For the gene groups classified in this way, we used AI-assisted literature searches to identify genes for which previous studies had reported functional validation directly related to adipocyte differentiation. The results are summarized below.

* The number of methods with P<0.05 was counted across four analyses: Student's t-test (FPKM), Student's t-test (counts), DESeq2, and edgeR.

What is interesting is that genes with previous functional validation were not found only among those that reached P<0.05 in all four analyses—that is, the genes most consistently identified as significant candidates across different analysis methods. Functional evidence was also found for genes that reached P<0.05 in three, two, or one method, and even for genes that did not reach P<0.05 in any of the four analyses.

At least among these 200 genes, this suggests that there is no simple relationship in which genes that are more consistently statistically significant are necessarily more biologically important.

The ChIP-Atlas analysis may also be worth revisiting from another perspective

In a previous analysis, we used ChIP-Atlas to compare transcription-factor binding candidates with RNA-Seq expression changes. The candidate target genes showed a mixture of increases and decreases, and it was difficult to identify a clear trend from scatter plots or individual Fold Changes alone.

The HuR study examined in this case study suggests another way to ask the question. What would happen if we compared the cumulative distributions of expression changes between transcription-factor binding candidates and non-binding genes?

Even when changes in individual candidate genes are small, the binding-candidate group as a whole may show some directional tendency. Visualizing the CDFs and comparing them with a Kolmogorov–Smirnov test might reveal patterns that were not apparent from scatter plots or individual Fold Changes.

This is not simply a matter of using a different statistical test. It means asking a different question of the same RNA-Seq data.

Does this mean conventional analysis methods were wrong?

The purpose of this article is not to reject conventional DEG analysis. Because it is impossible to experimentally validate every gene, it is entirely reasonable to narrow down candidates using P values and Fold Changes and to begin with genes showing larger or more reproducible changes. DESeq2 and edgeR were also developed to address clearly defined statistical problems.

Underlying this approach is the expectation that biologically important genes are more likely to be found among the top-ranked candidates, or that genes selected by P-value or Fold Change thresholds are enriched for biologically meaningful hits. This is a natural and reasonable idea.

However, the fact that these approaches have been widely adopted does not mean that the broader question of how RNA-Seq should be used in research has been settled. Methods now regarded as "standard" began as ideas proposed, tested, and published by individual researchers. Because they proved useful, they became widely adopted, implemented in software, and extended into many related methods.

At the same time, their success may have helped establish the assumptions that "analyzing RNA-Seq means identifying DEGs" and that "candidate discovery fundamentally means filtering genes by P values or FDRs." In the HuR study examined in this case study, a biologically meaningful group of genes was first defined using RIP-Seq, and RNA-Seq was then used to examine the behavior of that group as a whole. This is quite different from conventional DEG analysis, but biologically it is a very natural question.

We still do not know how omics data should be analyzed

RNA-Seq has many well-established analytical methods, including normalization, DEG analysis, PCA, clustering, GSEA, and pathway analysis. However, these methods do not themselves answer the broader question of how to extract the information that matters most for research from omics data.

Should we examine individual genes or groups of genes? Should we look only at RNA abundance, or combine it with translation, protein-level measurements, or other regulatory layers? Should we focus on genes with large changes, or on subtle shared changes within gene groups defined by other experiments? What we can see in the same dataset depends on the question we choose to ask.

A great deal of discussion in omics analysis has focused on statistical models, normalization methods, and DEG methods. By contrast, there is still much more room to discuss ideas such as, "Given this biological question, what if we looked at this group of genes in this way?" It may be increasingly important to bring the perspective of wet-lab researchers—who routinely formulate hypotheses, perform experiments, and interpret unexpected results—more deeply into discussions of data-analysis methods.

The fundamental discussion of how omics data should be used in research has barely begun. We need to recognize that there is still no general answer and encourage wet-lab and computational researchers to discuss how the data themselves should be viewed, rather than remaining confined to established analysis workflows.