Question: Analysis of RNA-seq time series with edgeR
gravatar for Francesca
20 months ago by
Francesca0 wrote:

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,


ADD COMMENTlink modified 20 months ago by Aaron Lun20k • written 20 months ago by Francesca0
gravatar for Aaron Lun
20 months ago by
Aaron Lun20k
Cambridge, United Kingdom
Aaron Lun20k wrote:

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.

ADD COMMENTlink modified 20 months ago • written 20 months ago by Aaron Lun20k

Hi Aaron, thanks a lot for your reply! I am indeed using this approach, but I wonder if there is a better solution than computing these three lists of DE genes and then merging them in with some logical rule (e.g. A1-list & A2-list & !B-list).

ADD REPLYlink modified 20 months ago • written 20 months ago by Francesca0

I don't think so. It's difficult to rigorously identify non-DE genes, because no detection may just be due to a lack of power. The approach I've described above is the easiest way to do it, and it's probably good enough for most applications.

ADD REPLYlink modified 20 months ago • written 20 months ago by Aaron Lun20k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 172 users visited in the last hour