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...
Thank you for a very illuminating answer.
It seems that I am using two topTags, one for each experimental factor, so it seems that the issue is there (I just misunderstood where it was in the edgeR procedure). I think I will use coef=5:6 and use a single topTags() to avoid this. How does that sound?
It seems that one of my experimental factors (factor.2) has no effect whatsoever (zero genes are differentially expressed at 0.05 FDR) when tested on its own. Is it still useful to check for interaction since that is the way the experimental design is set up? In other words, should I include an interaction term in my design matrix because that is the way the experimental design was done?
If so, my design matrix would look like this:
Does that make sense? Can I do a single glmLRT(coef=5:6:7) for all of them (factor 1, 2 and interaction) or are those done separately glmLRT(coef=5:6), glmLRT(coef=7)? If so, how can topTags() correctly manage the number of tests? Or was that no problem because of your explanation that "FDR should be controlled under the nominal threshold in each DE list.."?
With regards to factor 2; it's important to stress that failure to detect any DE is not the same as no effect. For example, your original additive model might have been insufficient to describe the experimental set-up. This could lead to inflated dispersions and loss of power to detect genuine DE for factor 2. Now, you've got plenty of residual d.f. to play with, so you might as well specify the full model and see what you end up with. The approach I suggested with
group
is simpler to interpret, but if you insist on an interaction model, you should do:If you run
glmLRT
withcoef=5:7
, your null hypothesis is that there's no difference between any of the four groups. If you runcoef=5
, your null is that there's no difference between +/+ and -/+. If you runcoef=6
, your null is that there's no difference between +/+ and +/-. If you runcoef=7
, your null is that there is no interaction, i.e., effects are additive. Which one you should use depends on the question you want to answer, rather than any consideration of FDR control.The question that I would like to answer is either:
(1) what genes are differentially regulated by any of factor 1 and factor 2.
(2) what genes are differentially regulated by any of factor 1, factor 2 or interaction between factor 1 and factor 2.
If I understood your answer correctly, I should use coef=5:6 for (1) and coef:5:7 for (2). Is that correct?That sounds right. A more intiuitive interpretation of aim 1 is that you're looking for genes that are DE upon changing either of the factors individually. In contrast, aim 2 looks for DE genes between any of the four groups. So, any changes in expression in the double negative group will not cause rejections for aim 1 (as this is absorbed by the interaction term), but will be picked up by aim 2.
And doing coef=5 and coef=6 in two different
glmLRT / topTags
would be looking for those that differ for factor 1 specifically and those that differ for factor 2 specificaly?Yeah. If you intepret it from group-based perspective, the DE genes for
coef=5
are those that differ between +/+ and -/+ (i.e., after changing the first factor), whereas the DE genes forcoef=6
are those that differ between +/+ and +/- (i.e., after changing the second factor).