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. 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). 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,1 0,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]]

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 >