Hi everyone! This is a bit long but I hope to state my problem more clearly

I've been trying to perform an analysis on the effects of a treatment over different time points using RNASeq with both limma and edgeR but I'm a bit confused with how to approach and interpret the results.

The experiment consists of samples taken at 5 different time points, with both treatment and control (this gives a total of 10 groups). Additionally, I have a baseline group (control at time 0). In total there are 11 groups. For each sample we have 4 replicates, so in total it's 44 libraries. Only one lane was sequenced (all libraries were multiplexed and sequenced in one single lane).

We wish to know what are the effects of the treatment at each time point (i.e., which genes are differentially expressed at time 1 in treatment vs control, which at time 2, and so on.)

First I tried using edgeR using a model that takes into account the time, treatment, interaction and an experimental effect of the tube strips where the libraries where made:

> design <- model.matrix(~Time+Treatment+Strip+Time:Treatment) 

I coded the times tm1-5. The control is a and the treatment is coded as b.

Here comes one of my first confusions. The design matrix only takes displays the interaction between the time and the treatment but show's no time:control. Also, I had to manually get rid of a term in the matrix (b) to make it full rank:

         (Intercept) tm1 tm2 tm3 tm4 tm5 b Strips2 Strips3 ... tm1b tm2b tm3b tm4b tm5b

For what I understand the last part of the design matrix accounts in part for what I'm looking for but I don't know how to interpret in relation to a treatment effect relative to a control in a particular time point.  

I tried two different approaches:

1) using an ANOVA for the coefficientes corresponding to the terms I wanted resulted only in 15 significant DE genes (FDR < 0.05)

> glmLRT(Fitglm, coef=12:16)

2) using matrix of contrasts I found ~700 DE genes.

> con <- makeContrasts('tm1'=tm1b - tm1, tm2'=tm2b - tm2, tm3'=tm3b - tm3, 'tm4'=tm4b - tm4,'tm5'=tm5b - tm5                    levels=design)
> glmLRT(Fitglm, contrast = con)

I next tried using limma converting the normalized read counts with voom and following section 9.6 of the limma vignette.

Using the same model as in edgeR gave no DE genes using coefficients and only 10 genes using contrasts: 

> design <- model.matrix(~0 + Time + Treatment + Strip + Time:Treatment, data = targets)
> fit <- lmFit(y, design)
> fit <- eBayes(fit)
#Contrasts:
> topTable(fit, coef=12:16,n=Inf)
#Coefficients:
> topTable(fit, coef=12:16,n=Inf)
>cont.dif <- makeContrasts(  Dif1 =Time1.Treatmentb-Time1, Dif2 =Time2.Treatmentb-Time2, (...), levels=design)
>fit2 <- lmFit(y, design)
>fit2 <- contrasts.fit(fit2, cont.dif)
>fit2 <- eBayes(fit2)

3) Following the Many timepoints section of the limma guide (9.6.2) and an advice from a post in http://permalink.gmane.org/gmane.science.biology.informatics.conductor/53752 also gave no DE genes using the contrast as adviced in the post nor using the advice from the limma guide.

X <- ns(as.numeric(targets$Time), df=4)

Treat <- factor(targets$Treatment)

>design3 <- model.matrix(~Treat+Treat:X)
>colnames(design3) <- make.names(colnames(design3))
>fit3 <- lmFit(y,design3)
>fit3 <- eBayes(fit3)

## different responses between the  different time courses.

cont.mat <- makeContrasts(
  Treatb.X1-Treata.X1,
  Treatb.X2-Treata.X2,
  Treatb.X3-Treata.X3,
  Treatb.X4-Treata.X4,
  levels=design3)
fit4 <- contrasts.fit(fit3,cont.mat)
fit4 <- eBayes(fit4)
topTable(fit4)

Has anyone had any situations like this? Are there significative differences between using coefficients in limma and edgeR? What would be the best approach for this case (testing the effects of treatment vs control at different times). Is the baseline group (control@ time0) of any use.

Any insights and suggestions are welcome!

### Edited to add sessionInfo() output:

> sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-apple-darwin15.3.0 (64-bit)
Running under: OS X 10.11.4 (El Capitan)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] edgeR_3.12.1 limma_3.26.9

loaded via a namespace (and not attached):
[1] tools_3.2.4

If you don't understand the interpretation of coefficients in an interaction model, I strongly recommend you avoid such a design matrix in favor of the equivalent group means design (ensure that Time and Treatment are both factors):

targets$Group <- interaction(targets$Time, targets$Treatment, drop=TRUE)
design <- model.matrix(~0 + Group + Strip, data=targets)

Now the first 11 coefficients will represent the means of the 11 combinations of Time and Treatment (after subtracting the batch effects of Strip). Knowing this, you can now phrase any contrasts you want as differences between group means.