Hi all,

We have an RNA-seq time-series over three time points, measuring the response of a cell line to three treatments (A1, A2 and B). We also have the control measurements, over the same time points plus the time 0h (treatments were applied after 0h). We have three biological replicates for each time point and treatment.

     Treat Time
1   Ctrl   0h
2   Ctrl   0h
3   Ctrl   0h
4   Ctrl   1h
5   Ctrl   1h
6   Ctrl   1h
7   Ctrl   2h
8   Ctrl   2h
9   Ctrl   2h
10  Ctrl   3h
11  Ctrl   3h
12  Ctrl   3h
13     B   1h
14     B   1h
15     B   1h
16     B   2h
17     B   2h
18     B   2h
19     B   3h
20     B   3h
21     B   3h
22    A2   1h
23    A2   1h
24    A2   1h
25    A2   2h
26    A2   2h
27    A2   2h
28    A2   3h
29    A2   3h
30    A2   3h
31    A1   1h
32    A1   1h
33    A1   1h
34    A1   2h
35    A1   2h
36    A1   2h
37    A1   3h
38    A1   3h
39    A1   3h

We would like to identify the genes that respond to treatments A1 and A2, but not to B. Treatments A1 and A2 might have different dynamics, so we cannot merge their labels in a single factor. What would be the best way respond to this question with edgeR? For now, we are not interested in highlighting different response dynamics. 

Thank you in advance,

Francesca

The simplest solution seems to be to parameterise the experiment as a one-way layout:

groups <- paste0(Treat, ".", Time)
design <- model.matrix(~0 + groups)

... and go through the (quasi-)GLM framework in edgeR. For each treatment (say, B), you set up a contrast like:

con <- makeContrasts(B.1h - Ctrl.1h, B.2h - Ctrl.2h, 
    B.3h - Ctrl.3h, levels=design)

... which compares the treatment and control groups at each time point. You use this contrast to do an ANODEV

res <- glmLRT(fit, contrast=con)

... which looks for any effect of the treatment. If you do this for all treatments, you can identify the DE genes that respond (in any manner) to treatments A1 and A2, and intersect them with non-DE genes that do not respond to B (e.g., take genes with p-values > 0.2).

If you want to look at specific dynamics, then you'll have to look at comparisons between treatments at specific time points. This can also be easily done by just testing each of the contrasts above separately, rather than supplying the whole contrast matrix to glmLRT.