I am doing a DGE analysis using RNAseq data to compare three conditions. I am using a standard pipeline (Create DGEList > Filter very lowly expressed genes > TMM normalize > DGE). Since there is a significant variation in library sizes (8-fold), I chose voom rather than limma-trend for DGE.

However, some of my top DEG are genes of low expression (e.g. Gene #1 : Voom-CPM = -5.0 in group A vs Voom-CPM = -1.0 in group B, hence contrast B versus A logFC =4 with adjusted-P = 10^-20). I am concerned that such fold-changes may be exaggerated. I am aware that voom calculates CPM using a prior.count of only 0.5, whereas the cpm() function uses a default prior.count of 2, leading to FC shrinkage.

Which approach would you recommend to estimate FC in these cases?

voom has already shrunk the logFC so they are not exaggerated. The raw fold-changes would be much larger.

If you are using voom on small counts, and if you have genes that are completely absent in one or more of the treatment groups, then I recommend edgeR::voomLmFit() instead of limma::voom().

Have you filtered low expressed genes using filterByExpr() as recommended in the limma workflows? Filtering is important for voom() although less necessary for voomLmFit().