We have no plans to offer confidence intervals for logFCs from edgeR.

There is actually a standard formula for the standard errors for the coefficients in any generalized linear model and, in principle, this could easily be computed in edgeR. We don't plan to do this however for two reasons:

1. Standard errors for count-based glms (binomial, Poisson, negative binomial) are notorious unreliable. We showed in the earliest edgeR paper (Robinson and Smyth, Biostatistics, 2008) that the t-tests obtained by dividing coefficients by their SEs are not even consistent (they do not necessarily become more significant as the effect size increases). We don't wish to include misleading information in the edgeR output. Likelihood ratio tests and exact tests are far, far more reliable.

2. The logFCs returned by edgeR are actually shrunk towards zero and have smaller SEs than ordinary glm coefficients would have. How to compute SEs for shrunk logFCs is not clear.

To compute bootstrap intervals, you would need a reasonably large number of samples. There is also a theoretical problem arising with the fact that the different RNA samples are not exchangeable in that they are associated with different sequence depths.

If you really do need confidence intervals for some reason, then I would suggest that you use voom and limma instead of edgeR. SEs are reliable for Gaussian models, so limma is able to return reliable confidence intervals for coefficients. Just set confint=TRUE when running topTable().

Thank you so much for this answer! I am building a random effects model to compare log fold change across several studies. I had used DESEq2 initially since they provide a log fold change SE and have been previously documented to be used with microbiome data (which is what I am working with).

However, in another set of data I used voom + limma since I had a paired design for the samples. Is there any way we could extract the SE from the confidence intervals in a voom + limma model? By extension, does it also imply that the log fold change standard error (lfcse) reported by DESEq2 are not stable since they follow an underlying negative binomial distribution?

Could you kindly elaborate why the logFC are shrunk towards 0 in edgeR / limma even for the non zero counts? I noticed this as my log fold changes in DESEq2 are somewhat higher in DESeq2 than in edgeR/ voom + limma for the same features. Why is that the case?

Thanks a lot for your time,

Manasi (PhD Epidemiology Candidate, UT SPH)