Hi there, I am trying to perform a differential expression analysis on large RNA-seq data (761 samples x 26000 genes), where the samples are divided into 5 classes, e.g. group_1 (n=275); group_2 (n=76); group_3 (n=151); group_4 (n=163) and group_5 (n=96). The use of DESeq2 and edgeR gives similar results, with a large number of differentially expressed genes (adjusted P-value <0.05) when factor levels are compared pairwise, representing about 40-50% of the genes in the whole dataset (which is too many for the phenotype in question). Their number drops to about 3-5% (expected proportion of differentially expressed) when I apply a threshold of |log2 fold-change| >0.5. I have read that limma-voom is recommended for large samples due to speed problems, but this is not a big issue at the moment with either DESeq2 or edgeR on this data. I was wondering whether there might be specific reasons to use alternative differential expression analysis methods to handle large samples and a variable of interest with many levels, or whether DESeq2 and edgeR also work effectively in this case. Thank you for any suggestions.

Just one note, in addition to the above pertinent pieces from ATpoint :

You should really consider to include RUV or SVA factors in such a large experiment. Even low correlations between batch variation and comparisons of interest will produce spurious rejections at such sample sizes. More here under Experiments with many samples

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs

And there is code for doing this in our workflow:

https://bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html#removing-hidden-batch-effects

I routinely work with dataset in the 100s and every time we model unwanted variation.

When you have large sample sizes and/or a lot of power you may want to test against a different Null hypothesis than "a gene has a fold change of 0", which is the default. You can use the lfc argument in DESeq2::lfcShrink or in edgeR::glmTreat for this, thereby enriching for genes with good evidence to have fold changes beyond the value you choose, for example 1.5. Other than that, you could use limma-voom which scales better at large sample size. The function for a fold change testing would be limma::treat here. As for any analysis, be sure that you explore the data first, e.g. by PCA or MDS plots to ensure that there is no unwanted confounding of any sort (batch effects for example) that confound the analysis.