Analysis of RNA-seq time series with edgeR
Entering edit mode
Francesca • 0
Last seen 5.7 years ago

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,


edger rna-seq differential expression timecourse time course • 2.8k views
Entering edit mode
Aaron Lun ★ 27k
Last seen 8 hours ago
The city by the bay

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.

Entering edit mode

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).

Entering edit mode

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.


Login before adding your answer.

Traffic: 242 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6