Hello,

I am trying to perform DE gene analysis on RNA-Seq data. My groups consist of animals from two different populations (SD v. XE) and two different treatments (CON v. E2). I'm trying to look at nested interactions and am following section 3.3.2 in the edgeR user guide.

My metadata 'meta' is as follows:

   samples treat pop  group
1   THXL1A    E2  SD  E2.SD
2   THXL1B   CON  SD CON.SD
3   THXL1C   CON  XE CON.XE
4   THXL1D    E2  XE  E2.XE
5   THXL1E   CON  SD CON.SD
6   THXL1F    E2  SD  E2.SD
7    THXL7   CON  XE CON.XE
8    THXL8    E2  XE  E2.XE
9    THXL9   CON  SD CON.SD
10  THXL10   CON  SD CON.SD
11  THXL11    E2  SD  E2.SD
12  THXL12    E2  SD  E2.SD

My script looks like this:

meta$treat <- relevel(meta$treat, ref="CON") #sets reference level as 'CON'
design <- model.matrix(~treat + treat:pop, data=meta)                                                                                                                   
dispData <- estimateDisp(dgeData, design, robust=TRUE)
fit <- glmFit(dispData, design)

However, when I look at  'design,' I only get 4 coefficients (missing "treatCON:popSD" and "treatE2:popSD"):

[1] "(Intercept)"    "treatE2"        "treatCON:popXE" "treatE2:popXE"

As opposed to the 6 coefficients they list in the user guide.

Please help!

With 2 populations and 2 treatments, there are 4 possible combinations of population and treatment. Hence, you have 4 coefficients. If you're not sure how to interpret these 4 coefficients, consider the advice given in the very first sentence of the section of the user's guide that you cited: "We generally recommend the approach of the previous section [3.3.1], because it is so explicit and easy to understand."