[BioC] Setting up model for time course analysis

[Heidi Dvinge]

Hi Casper,

>
> Dear Bioc,
>
>    I am trying to find DE genes in a time course experiment, and I am a
> bit confused at how to setup the matrix model. I wasn't able to follow
> the example in manual page 48-49 because I don't know what the variable
> "targets" refer to. But I have setup a model based on Section 7.2.
>
I guess you're referring to the limma users' guide. "targets" probably
refer to the Targets.txt file that has been read in. Just do a search
through the pdf and it should pop up.
>
>
>   My layout is: I have 10 time points, time 0,1,3,5,7,9...
>
>   Each time point has two replicates.
>
>   My aim is to find DE genes by comparing each time point to the reference
> point (i.e. at time 0).
>
Not directly related to your question, but when you have such a long time
series you might want to consider doing something else than pairwise
comparisons in limma. For example, the "timecourse" package uses an
approach similar to limma, but is specifically designed for time course
experiments, and takes into account that the individual time points are
not independent.
Comparing all time points directly to t0 does give you a list of genes
that respond during one or more time points, but the downstream
interpretation might become more difficult. For example, what if you get a
gene that is significant at 12, 18 and 30 hours but not at 24? Is there
actually a dip in expression level, or is the p-value at t24 just slightly
higher than whatever threshold you use? If you get a gene that is
significant only at a single time point, is this then a genuine spike in
expression level? By doing all these pairwise comparisons I suspect you'll
get a lot of genes with a "mixed" behaviour like this. At least if you
rely on some fixed p-value threshold and divide genes into differentially
expressed/non-differentially expressed at each time point, rather than
doing some clustering across the entire time series.

HTH
\Heidi
>
>
> Below are my codes related to setting up the model. It would help me a lot
> if someone could confirm if I'm doing this right, or whether some changes
> that need to be done.
>
> Thanks a lot in advance!
>
>
>
> design<-model.matrix(~0+factor(c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10)));
>  ##1=at time 0, 2=at time 1...and so on
> colnames(design)<-c("WT0hr","WT1hr","WT2hr","WT6hr","WT12hr","WT18hr","WT24hr","WT30hr","WT48hr","WT72hr");
> fit<-lmFit(allDatagcRMA, design); ##allDatagcRMA contains the expressions
> normalized by gcrma
> contrast.matrix <- makeContrasts(WT1hr-WT0hr, WT2hr-WT0hr, WT6hr-WT0hr,
> WT12hr-WT0hr, WT18hr-WT0hr, WT24hr-WT0hr, WT30hr-WT0hr, WT48hr-WT0hr,
> WT72hr-WT0hr, levels=design);
> fit2<-contrasts.fit(fit,contrast.matrix);
> fit2<-eBayes(fit2);
> topTableF(fit2, adjust="BH"); ##To get the top 10 genes that respond to at
> least one of the comparisons above
>
>
>
>
>
> Thanks again!
>
> Sincerely,
>
> Casper
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>