I am currently learning to do differential expression analysis on bulk and pseudobulk RNAseq data with limma::voom. In the function documentation it says
Note that edgeR::voomLmFit is now recommended over voom for sparse counts with a medium to high proportion of zeros.
However, edgeR::voomLmFit does not seem to be very often used by the community (yet?), is not included in the standard limma tutorials, and is not an option in the pseudo-bulk differential state analysis function of the muscat package (muscat::pbDS), which makes me a bit hesitant.
Is there a specific reason to not use voomLmFit at least if the data is sparse, if not always, instead of limma::voom + limma::lmFit? The paper and the documentation don't mention any downsides.