Pairwise differential expression comparisons can work, up to a point. ANOVA may be another option. However, I have used with much success a multinomial lasso-penalised regression for this, previously. If done correctly, it should identify genes that are more highly or lowly expressed in just one stage.

You can get more advanced by fitting trend models to the data - this is an entire area in its own right, and there are different ways to interpret and use these.

A much simpler way is to transform the entire dataset into Z-scores and take Z>3 as being statistically significantly expressed above the mean (Z of 3 is equivalent to p=0.05 on a 2-tailed distribution).

Kevin

Dear Dr. Kevin I am thankful for your exact answer. that was great and useful. I would be appreciated if you attach an example of what you have performed before with more details. for instance, the R codes and packages you have used would be suitable. For example, this post was fascinating to me "Tutorial: Clustering of DAVID gene enrichment results from gene expression studies".

with best regards.