I appreciate that these sorts of benchmarking papers exist to get an overview, but I will die on my hill saying that the specific combination of tools and these scores that are calculated there always depend on the dataset and situation, so I would not generally say something like "DESeq2+LTR+something" is best. I agree generally that also in my hands pseudobulking (if experimental design allows and sample size and depth is good) feels more robust that DE on pure single-cell level.
That having said, I personally for some time now almost exclusively use limma (voom or trend depending on situation) since it is the most generic one of common choices. I pair this with filtering by percent expression, meaning that I only test genes detected in >= x % of cells of at least one group in a pariwise comparion. Usually x is somewhere between 10-25%. That in my hands protects strongly from spurious calls stemming from outliers and reduces multiple testing burden a lot.
I would (if possible) always do paired tests, so e.g. within donor, within mouse etc, as this is most powerful. I typically use a test against a fold change (so not LRT-like tests) and I test against a fold change > 0 to avoid significant genes with small effect size. In limma/edgeR this is done via treat/glmTreat. Choice of minimum FC is arbitrary but something in the ballpart of 1.1-1.5 usually does the trick to emphasize what in my head are "meaningful" changes, avoiding fold change filtering which favour low-count genes (many threads where the limma/edgeR author recommends against fold change filtering).
Individual Qs:
Why is edgeR reporting so many DEG compared to DESeq2 in the B1vsB2 comparison? Is is more sensitive to small differences? Is it more prone to false positives?
Generally they perform at least similar, one would need to see code, data and prefiltering strategy to decide. After all, hard cutoffs can lead to overestimation of differences, you should plot plvalues and fold changes against each other to see how they correlate if you really want a direct comparison. I would check magnitude of fold changes and prioritize by fold change (as I discuss above) to avoid many DEGs hard to interpret.
Should we use all (22) pseudosamples in the data object, or only the 2+2 of the comparison of interest?
The usual "all together" vs "separate" discussion. For single-cell I usually run separate because we often assay tissues or heterogeneous populations so the mean-variance trend when using all samples is probably not too representative for each individual comparison.
Should we use a simple design ~cluster, or is it advisable to use ~cluster + sample for pseudobulk analyses trying to find differences between cell types? One of our samples always has more cells (of all cell types) than the other. Should we also include the number of cells per sample as covariates?
When paired analysis is possible I would do that, as it can eliminate donor variations etc and is more powerful, allowing to do more aggressive fold change testing as said above.
One of our samples always has more cells (of all cell types) than the other. Should we also include the number of cells per sample as covariates?
I am not sure how n of cells as covariate is behaving, but it's true that number of cells is important as fewer cells per pseudobulk gives lower counts, and this can not always be eliminated by the usual normalization methods. I usually run a subsampling function and compare DE numbers when using all cells or when subsampled to the smallest pseudobulk to see if absence if DEGs is due to low counts or "real".
Is it worth using LRT methods over Wald/Exact for simple pairwise comparisons? Is it more correct to use reduced=~1 or reduced=~sample in this type of analysis?
In my head testing against a fold change is more intuitive but this is probably because I have no formal stats education. Avoiding the garden of forked paths I never even explore LRT-like tests because after all, what is the strong and reproducible benchmark to decide for one over the other? So I simply stick to fold change testing which is the default for both DESeq2 (Wald) and edgeR (glmQLFTest as per their user guide).
I wouldn't mind running both DESeq2 and edgeR, and use only the DEG and pathways that are common in the results of both methods. However, how should I match the choice of parameters in both packages?
Ask yourself if inflating the burden of interpreting DEGs by running two tools is worth it. Do you plan to combine into a meta-analysis or just keep both lists? I would not do that. I assume you want to reveal interesting biology and not bother with benchmarking tools, no? Choose one and stick with it. Both edgeR and DESeq2 are well-established and there is not (to my knowledge) a reason to favour one over the other in general, as pseudobulk is basically bulk RNA-seq (with lower counts) and both tools work well with it, many benchmarks have shown it.
