Hello,

I am carrying out a DESeq2 analysis of a dataset with 16 groups and 4 replicates per group. I applied the variance stabilizing transformation, removed batch effects using limma::removeBatchEffect and plotted the PCA results. The PCA plot shows high variability within groups:

I ran a differential expression analysis with all samples together and identified up to 8 differentially expressed genes. However, when I split the data into pairs of groups (creating a separate DESeqDataSet object for each treatment), I found many more significant differences, with 6 groups having between 100 to 250 differentially expressed genes. I know that the DESeq2 vignette recommends to split the data into pairs of groups when a particular treatment has much higher within-group variability. What would be the recommended approach when most groups have high within-group variability?

I also performed other analyses using sleuth and limma, and it seems only limma::voomWithQualityWeights identifies more than 100 differentially expressed genes in some groups when running samples from all groups together.

Any help is greatly appreciated.

I don't have any particular advice here, but I'd warn that it seems like there aren't many differences and that by trying many options you may end up losing control of FPR simply by iterating across many combinations of methods.