Hi,

I am in the process of analyzing RNA seq data for an experiment in ants with three factors.

Factor 1: presence/absence of queens (levels:present, absent)

Factor 2: sample type (levels: forager, nurse, larva, etc.)

Factor 3: stage (1, 2, 3, 4, or 5)

I have been following along in the user manual, mostly following section 3.3.3 and 3.3.4. 

We want to be able to analyze the effect of all factors, as well as pairwise interactions. For questions using all the samples, the guide seems straightforward, but we also want to ask questions like how do forager samples vary according to stage WITHIN queen-absent colonies. So this is basically looking within just queen-absent colonies. 

There are two ways I can think of doing this.

What I've been doing: subset my treatments file and counts file so they only contain queenless samples, and proceed from here.

What I think may be better: build the full model, with all interactions, i.e. design=model.matrix(~queen_presence*type*stage). Then just extract the comparisons we want. I understand how to do this with pairwise comparisons. However, what I would like to be able to do is ask if there is an effect at any time, as in section 3.3.3. However, this requires requires reference to a control-type sample, coef=1, correct? Is it possible to ask the same question, but use a different coefficient as the control, i.e. compare something like coef=6:8 to coef=2 rather than 1?

The reason I would like to use the second method is because it seems to me it could have more power. This way, I fit a glm with all possible information and then just extract what I want, rather than neglecting information in the first approach.

Is there a good way to do the second approach, or should I stick with the first?

Thanks!

For such an experiment, I would find it easiest to use the group means parametrization:

group <- factor(make.names(factor1:factor2:factor3))
design <- model.matrix(~0+group)

And then use makeContrasts to test for differences between specific groups.