Analyzing expression time course with DESeq

[john_nl]

Hello,

I have 11 samples covering 10 time points. Time point 1 has a control sample and a treated sample. The following 9 time points are all treated samples.

I wish to be able to compute the ratios between the control sample and the treatment samples at each time point and then plot these ratios in order to observe the general expression trends over time for each gene. Ultimately I want to cluster these profiles in order to identify genes that roughly behave in a similar manner over time.

DESeq can compute fold changes between two samples and it has the getVarianceStablizedData() function to normalise the data such that it has similar variance (for clustering). In the vignette, the example is with two samples- is it OK to perform the getVarianceStablizedData on the entire set of samples and then compute moderated log fold changes (as in the vignette) at each timepoint (treatedTimepoint_x - control)? For each gene I can then generate a series of moderated fold changes- is it then appropriate to cluster genes based on these profiles?

[Simon Anders]

Hi John

Yes, your analysis problem is precisely the kind of tasks which we had in mind as potential applications for the variance stabilizing transformation. So, go ahead and give it a try.

Of course, this is not a perfect solution. The VST makes the data only _approximately_ homoskedastic, and it depends on circumstances whether this is good enough. So, inspect your results carefully and think about some sanity checks. The VST works worse, by the way, if your size factors are rather different; I'm still thinking about how to straighten this out.

Finally, it will make a difference whether you call 'estimateDispersions' with method 'blind' or 'pooled'. I would, in general, recommend the former if you want to cluster samples (because then you want to be unbiased about the samples' relations) and the latter if you cluster by genes (because then it is appropriate to use sample information). However, you cannot use "pooled" because you don't have real replicates, and "blind" will overestimate variance and so might skew the VST in some direction. In principle, a clustering analysis is proper (though not ideal) without replicates, but the VST will have issues handling it well. Hence, again, your mileage will vary.

In general, analysing RNA-Seq time courses seems to be an open question to me. I have not seen any papers yet suggesting good analysis methods, and there clearly is still work to do for us biostatisticians before we can recommend some standard analysis method.