My experimental design consists of four groups and two experimental factors (+/+, +/-, -/+ and -/-). In addition, each group has four biological replicates (the first biological replicate of each group was taken at the same time, the second replicates at another and so on) and a clear batch effect between replicates.

My design matrix looks like this:

design <- model.matrix(~batch+factor.1+factor.2)

             (Intercept) batch2 batch3 batch4 factor.1 factor.2
+/+.1           1      0      0      0               0      0
-/+.1           1      0      0      0               1      0
+/-.1           1      0      0      0               0      1
-/-.1           1      0      0      0               1      1
+/+.2           1      1      0      0               0      0
-/+.2           1      1      0      0               1      0
+/-.2           1      1      0      0               0      1
...
...
attr(,"assign")
[1] 0 1 1 1 2 3
attr(,"contrasts")
attr(,"contrasts")$batch
[1] "contr.treatment"

attr(,"contrasts")$factor.1
[1] "contr.treatment"

attr(,"contrasts")$factor.2
[1] "contr.treatment"

The factor.1 column represents presence/absence of factor 1 and factor.2 column represents presence/absence of factor 2.

So far, I have performed statistical significance testing between presence/absence of factor 1 and presence/absence of factor 2 in isolation:

lrt.factor1 <- glmLRT(fit, coef=5)
lrt.factor2 <- glmLRT(fit, coef=6)

However, this seems to mean that I am only correcting for x number of statistical significance tests (where x is the number of genes that survive filtering) in each test, but in total I am actually doing 2x tests spread over two separate applications of the glmLRT() function.

Is this really appropriate? It looks like the correction for multiple testing is artificially weak by this approach (since obviously the first glmLRT() doesn't "know" that I am running a second one after) and would thus produce more false positives / not be able to control FDR at the alleged level.

Is it more correct to do all the tests at once and correct for 2x number of tests instead?

lrt <- glmLRT(fit, coef=5:6)

Does the above command test to see if there are genes that are differentially expressed in any of the factor 1 (+/-) or factor 2 (+/-) comparisons and correct for 2x tests?

I have looked in the edgeR user manual and this kind of logic seems to be applied in the Arabidopsis case study when testing for differential expression between any of the batches (b1 vs b2, b2 vs. b3, b1 vs. b3) to see if there is a batch effect and then they use glmLRT(fit, coef=2:3), but I have only worked with R for a couple of weeks so a bit plagued by self-doubt...

Firstly, the Benjamini-Hochberg correction is applied in topTags, not glmLRT. Secondly, the FDR should be controlled under the nominal threshold in each DE list, when you run topTags on each object returned by glmLRT. There's no loss of control here, no matter how many DE lists you have. For example, if I get 100 DE genes for coef=5 and 500 DE genes for coef=6 at a FDR of 5% in each list, I'd expect no more than 5 + 25 = 30 false positives out of a total of 600 discoveries (any common genes between the two lists are considered as distinct discoveries). This gives an overall FDR of 5%, so no liberalness is observed.

Now, your coef=5:6 approach addresses a different question. You are correct in saying that it looks for genes that change between any of the experimental factors. However, this is done with a single test for each gene. Thus, the number of tests is still equal to the number of genes. Which approach you should use depends on what question you're trying to answer. If you want to find DE genes that change between states for factor 1, then it makes sense to set coef=5 (and similarly, coef=6 for factor 2). If you want to find genes that change between states for any factor, then your coef=5:6 approach is better.

Finally, your design matrix doesn't actually reflect a four-group setup. As it is now, it assumes that the log-fold change is additive between the two experimental factors. This may not be appropriate, e.g., for cases involving epistasis or synergy. In such cases, a better formulation might be:

group <- paste0(factor.1, ".", factor.2)
design <- model.matrix(~batch + group)