Experimentally validated genes may not stand out in RNA-Seq
In RNA-Seq analysis, differentially expressed gene (DEG) analysis is commonly used to identify genes with small P values or FDRs as candidates for further investigation. P values and FDRs are useful for narrowing down candidates, but biologically important genes do not necessarily appear at the top of the list. What we usually want to know is not simply which gene is the most statistically significant, but which molecules are closer to the cause of the phenomenon and how their effects lead to the observed phenotype. DEG rankings alone cannot tell us where a gene lies in this chain of events.
Here, we will look at GSE124280, which contains data from a study of the RNA-binding protein HuR and adipocyte differentiation.
HuR suppresses adipocyte differentiation
The study showed that reducing HuR promotes adipocyte differentiation, whereas increasing HuR suppresses it. Mice with adipose-specific HuR knockout also showed increased fat mass and metabolic abnormalities.
This leads to the next question: through which genes does HuR suppress adipocyte differentiation?
RNA-Seq shows what happens after HuR loss
The authors performed RNA-Seq on BAT, eWAT, and iWAT from HuR knockout and control mice. They also used GSEA to examine which pathways were altered in each tissue. Different changes were observed across tissues, including adipogenesis, inflammation, browning, and myogenesis.
RNA-Seq and GSEA, however, mainly show what happens across the tissue after HuR is lost. Additional experiments are needed to determine which mRNAs are directly regulated by HuR. The pathways detected by GSEA are therefore more naturally interpreted as downstream changes following HuR loss rather than as direct HuR targets.
RIP-Seq was used to identify RNAs directly bound by HuR
Because HuR is an RNA-binding protein, the authors used RIP-Seq to identify RNAs associated with HuR. Genes were ranked according to their enrichment in the HuR immunoprecipitation, and the top 200 were treated as HuR targets.
These 200 genes were not selected because they had small P values in RNA-Seq. They were selected using a different type of information: whether their RNAs were enriched in the HuR immunoprecipitation.
The 200 HuR target genes do not stand out in RNA-Seq
When these 200 genes are examined in the RNA-Seq data, most do not show large differences between control and knockout. In a scatter plot, many lie close to the diagonal, and highlighting the HuR targets does not reveal an obvious distinct group.
When the 200 genes are considered as a group, however, there is a slight shift toward lower expression in the knockout. The authors generated cumulative distribution functions (CDFs) of the HuR-FKO/control fold changes for the 200 RIP-Seq targets and for the other genes, and compared the two distributions using a Kolmogorov–Smirnov test. The HuR targets were shifted toward the downregulated side.
Rather than asking whether each individual gene showed a large expression difference, this analysis asked whether the distribution of the 200 genes as a group was slightly shifted. The Kolmogorov–Smirnov test evaluates differences between two cumulative distributions. In this case, the genes were not simply classified one by one as significant or non-significant; instead, the analysis evaluated whether the HuR target group as a whole tended to shift in one direction.
RNA abundance of the 200 HuR target genes. Left: RNA-Seq data from WT and HuR knockout eWAT. Black dots indicate the 200 HuR targets selected by RIP-Seq. Most lie close to the diagonal, and no large difference is apparent in the scatter plot. Right: Expression distributions of the individual samples.
INSIG1 did not stand out in RNA-Seq
Among the HuR targets, the authors focused on INSIG1. INSIG1 was already known to be involved in adipocyte differentiation, making it an interesting functional candidate.
In the RNA-Seq data, however, the change in INSIG1 was not especially large. When we reanalyzed the data, the direction of change was consistent across BAT, eWAT, and iWAT, but the t-test did not yield P<0.05.
If candidate genes had been selected only from those with small P values in DEG analysis, INSIG1 might not have been included.
Volcano plot comparing WT and HuR knockout in eWAT. Black dots indicate the 200 HuR targets selected by RIP-Seq, and the reddish-brown dot indicates INSIG1. Most HuR targets do not occupy particularly prominent positions in the RNA-Seq DEG analysis, and INSIG1 was also not among the genes showing either a large expression change or a small P value.
INSIG1 was then validated experimentally
The authors went on to investigate INSIG1 in more detail. They examined HuR binding to INSIG1 mRNA, the involvement of the INSIG1 3′ UTR, the reduction in INSIG1 mRNA stability after HuR loss, and the effects on protein abundance and adipocyte differentiation.
These experiments supported a mechanism in which HuR stabilizes INSIG1 mRNA and thereby suppresses adipogenesis. INSIG1 was neither the gene with the largest RNA-Seq change nor the one with the smallest P value, yet it was experimentally validated as an important gene in the study.
Translation efficiency also shows no large change
The study also performed Ribo-Seq and examined translation efficiency. The authors reported that HuR knockout did not produce a significant change in translation efficiency.
In the scatter plot, control and knockout look very similar, and even when the 200 HuR targets are highlighted, there is no obvious separation. When we examine the same 200 genes as a distribution, however, a slight shift appears to be present.
This resembles the RNA-Seq result. RNA abundance also showed only a small difference in the scatter plot, but a shift became detectable when the distribution of the 200 genes was compared. For translation efficiency, however, the authors did not detect a significant change.
Translation efficiency of the 200 HuR target genes. Left: Comparison of translation efficiency between WT and HuR knockout. Black dots indicate the 200 HuR targets selected by RIP-Seq, and no clear separation is visible in the scatter plot. Right: Distributions in WT and KO. The authors reported no significant change in translation efficiency, although a slight difference in the distributions is visible.
A small shift in RNA abundance was detected statistically, whereas the change in translation efficiency was not considered significant. Still, it may be premature to conclude, based on the P<0.05 threshold alone, that a biological effect was present in one case while nothing happened in the other.
Small changes can accumulate across different regulatory levels
Protein abundance is not determined by mRNA abundance alone. It is influenced by multiple regulatory levels, including translation efficiency and protein stability. Even if the change at each level is small, changes in the same direction can combine to produce a larger overall effect. Changes in opposite directions can also compensate for one another.
If the affected protein is an upstream regulator, a small change can then propagate to many downstream molecules and eventually lead to a larger phenotypic effect.
From a systems biology perspective, it is not surprising for a relatively small change upstream to produce a larger effect downstream. Transcription factors, RNA-binding proteins, kinases, and receptors can all transmit relatively small changes to many downstream molecules.
It is therefore difficult to judge biological importance from an individual result alone simply because the RNA-Seq fold change is small or because the change in translation efficiency is small.
There may be other important candidates among the 200 genes
INSIG1 was investigated in detail in this study, but it was not the only HuR target among the 200 genes showing an RNA-Seq change of this magnitude. Its known relationship with adipocyte differentiation likely also contributed to the decision to investigate it further.
If the 200 genes as a whole show a weak shift in RNA abundance and a small shift in translation efficiency, the same list may contain other genes with effects comparable to INSIG1, or even stronger effects at another regulatory level.
For example, some genes may show a small decrease in RNA abundance with little change in translation efficiency. Others may show small changes in the same direction at both the RNA and translation levels. Still others may show an RNA-level change that is partly compensated for by translation efficiency.
The choice of which candidate to test next therefore cannot be determined from P-value ranking alone. Known function, direct binding, RNA abundance, translation, protein level, and other information need to be considered together.
What are we actually seeing in enrichment analysis?
GSEA can detect a coordinated trend even when the changes in individual genes are small, as long as genes belonging to the same pathway tend to move in the same direction.
In this study, GO enrichment analysis of the 200 HuR targets identified significant enrichment of processes related to RNA metabolism and RNA processing. In contrast, there was little overlap between the 200 HuR targets and the genes contained in pathways that were upregulated in the RNA-Seq GSEA. The authors therefore considered the inflammatory pathways and other signals detected by RNA-Seq GSEA more likely to reflect downstream changes following HuR loss than direct HuR targets.
The interpretation of enrichment analysis depends on which gene set is used as input. When enrichment analysis is applied to genes selected from standard RNA-Seq because they show differential expression, it may reveal downstream functional changes caused by the experimental condition. In contrast, when it is applied to candidate genes that are closer to the underlying regulatory mechanism, such as targets identified by RIP-Seq or ChIP-Seq, it can provide clues about the molecular mechanisms or functions of the regulator itself.
In a previous analysis using 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. In the HuR example here, however, a weak shift was detected by comparing the cumulative distribution of expression changes between HuR-bound genes and the other genes. It may therefore be worth revisiting the ChIP-Atlas analysis by visualizing the cumulative distributions of binding candidates and other genes and comparing them with a Kolmogorov–Smirnov test. Such an analysis might reveal trends that were not apparent when the genes were examined individually.
What matters is what we measure and how we connect the results
This study did not rely on RNA-Seq alone. RIP-Seq was used to examine which RNAs were bound by HuR, RNA-Seq to determine how expression changed after HuR loss, and Ribo-Seq to examine changes in translation efficiency. Protein-level measurements and functional assays were then used to investigate whether those changes were connected to molecular function and phenotype.
Each experiment measures something different. Research requires deciding which aspect of a phenomenon should be captured by which experiment, and then determining how the results should be connected. The same applies to data analysis. If the biological question has not been defined, simply running DEG analysis, GSEA, and clustering in sequence becomes little more than a routine workflow.
Recent advances in AI have made the practical work of data analysis easier than before. Writing R or Python code, running DEG analysis or GSEA, organizing candidate genes, and searching the literature have all become easier.
What still cannot be determined simply by automating the analysis is what should be measured, how the different results should be connected, which observations are closer to the cause and which are downstream, and what should be tested next. Deciding what to measure, connecting results from different experiments, and formulating the next question—designing the research itself will likely remain a role for researchers.
