Thanks very much for the thorough explanation, I will add tm:trt1 to my model to see expression changes in two condition. Best wishes, On Tue, Mar 20, 2012 at 7:51 PM, James W. MacDonald <jmacdon@uw.edu> wrote: > Hi Xuan, > > > On 3/20/2012 2:13 PM, shao chunxuan wrote: > >> Hi Jim, >> >> Thanks for the reply, but I still have few more questions on the >> regression model of limma. >> >> Time points are not evenly spreaded, the real experiments like this: >> >> xx.1 <- c(0.5,1,3,6,12,24,72) >> time <- rep(xx.1,2) >> >> The design matrix is >> trt <- factor(rep(0:1, each = 7)) >> design=model.matrix(~time+trt) >> colnames(design) <- c("Intercept","time","treat") >> >> > design >> Intercept time treat >> 1 1 0.5 0 >> 2 1 1.0 0 >> 3 1 3.0 0 >> 4 1 6.0 0 >> 5 1 12.0 0 >> 6 1 24.0 0 >> 7 1 72.0 0 >> 8 1 0.5 1 >> 9 1 1.0 1 >> 10 1 3.0 1 >> 11 1 6.0 1 >> 12 1 12.0 1 >> 13 1 24.0 1 >> 14 1 72.0 1 >> attr(,"assign") >> [1] 0 1 2 >> attr(,"contrasts") >> attr(,"contrasts")$trt >> [1] "contr.treatment" >> >> Here is my understanding of what happens in limma >> >> In this case, Iimma is supposed to fit: >> y=p0 + p1*time + p2*treat, and y=p0 + p1*time, >> while time is quantitative variable, treat is group variable. >> A F-statistic will be calculated and the corresponding pvalue will be >> used to determine significance. >> >> Is it right to understand in this way? >> >> The implementation is : >> fit <- lmFit(ts.data, design) ## ts.data is the time-course data. >> fit <- eBayes(fit) >> xx.1 <- topTable(fit, coef="treat", adjust="BH") >> xx.2 <- topTable(fit, coef="time", adjust="BH") >> >> So xx.1 means the differentially expressed genes between the two >> treatments, after adjusting for the time effect, right? what does xx.2 mean? >> > > xx.2 gives you those genes where the slope of the fitted line is > significantly different from zero. In other words, these are the genes that > appear to go up or down as time progresses. > > You should note two things here. First, since you don't have replicates, > we are forced to use the time points as continuous covariates rather than > factor levels. This means that you are looking for a linear response to > time (e.g., genes that go up or down more as time progresses). This is not > normally what one would want to do, as it is more likely that a gene will > go up for a bit and then go back down (or vice versa). To account for that > curvature in the expression, you might want to add a time^2 covariate as > well. > > Second, linear regression is in some ways more complicated than analysis > of variance, and if you were just fitting a single model you would be able > to assess how well your model was conforming to the assumptions of > regression, test to see if you need a quadratic time term, look for data > with high leverage that might be skewing your results, etc. For ANOVA, you > are in general just testing to see if the mean of the groups are different. > > When doing microarray analysis, you just fit what you hope is a generally > reasonable model and let it rip on thousands of genes at once. There is no > way to check all the little details of a linear regression on thousands of > models, so you may (will?) have numerous genes that are clearly > mis-specified by the model you fit, or have glaring inconsistencies that > cause them to appear significant when it is just that the model is really > wrong. > > So that is my long winded way of saying that you are not IMO analyzing > your data in an optimal way, but you are forced to do so by the lack of > replication. So for each 'significant' gene you find, you should look at > plots of expression vs time prior to attempting validation. > > > > I am more interested in finding genes changed in time course in >> treatment, after adjusting for the control. >> Is tm*trt in design2 what I need? >> > > Not really. The model you are fitting here (xx.2) gives the genes that > change linearly with time, after adjusting for control. The assumptions you > are making are that > > a) Genes go up or down in a linear fashion (e.g., you can plot expression > values on the vertical axis and time on the horizontal axis and then draw a > line over the data, and it looks 'reasonable'). > b) The only difference between control and treatment is a difference in > the intercept, so you assume the same slope for control and treated samples. > > adding a tm:trt1 coefficient allows the slopes to be different as well. > This coefficient tests for a significant difference between slopes (for > example, it might be that geneX goes up in treated, but down in control). > > Best, > > Jim > > > >> Best, >> >> >> >> On Tue, Mar 20, 2012 at 2:48 PM, James W. MacDonald <jmacdon@uw.edu<mailto:>> jmacdon@uw.edu>> wrote: >> >> Hi Xuan, >> >> >> On 3/20/2012 6:50 AM, Chunxuan Shao wrote: >> >> Hi everyone: >> >> I have one microarray data set considering differentiation in >> a cell line. About 7 time points for both control and >> treatment, no replicates >> I would like to use limma to find the differentially expressed >> genes between time points for control and treatment, and want >> to compare the gene expression between control and treatment >> for the same time point. But I don't know how to make the >> right design and contrast matrix. >> >> After searching mail archive, the closest related answer is >> "https://stat.ethz.ch/**pipermail/bioconductor/2010-** >> June/033849.html<https: stat.ethz.ch="" pipermail="" bioconductor="" 2010-j="" une="" 033849.html=""> >> ", >> which suggest: >> >> " >> time=1:10 >> design=model.matrix(~time) >> " >> >> In my case, is it correct to set this? >> >> time=rep(1:7,2) >> design=model.matrix(~time) >> >> >> Probably not. This implies that your time points are equally >> spaced, with one time period between them (e.g., you collected >> samples after 1,2,3,4,5,6,7 hours or days or some other period). >> If you didn't use equally spaced time points, you need to change >> to reflect that. >> >> In addition, you need to set the design matrix up to include the >> control/treatment comparison. If I assume you are in fact using >> seven equally spaced time intervals, then you would want: >> >> tm <- rep(1:7,2) >> trt <- factor(rep(0:1, each = 7)) >> >> design <- model.matrix(~tm+trt) >> >> And the trt1 coefficient estimates the difference in the >> intercepts, which is what you are looking for. In other words, >> this is fitting a model where you are allowing the two sample >> types to have different intercepts, but assuming the same slope. >> If the intercepts are different, then one sample type has overall >> higher expression than the other. >> >> You could also allow for different slopes by adding a >> time/treatment interaction term: >> >> design2 <- model.matrix(~tm*trt) >> >> Here the main coefficient of interest would be tm:trt1, which >> measures the difference in slope between the treatment and control. >> >> Best, >> >> Jim >> >> >> >> design >> >> (Intercept) time >> 1 1 1 >> 2 1 2 >> 3 1 3 >> 4 1 4 >> 5 1 5 >> 6 1 6 >> 7 1 7 >> 8 1 1 >> 9 1 2 >> 10 1 3 >> 11 1 4 >> 12 1 5 >> 13 1 6 >> 14 1 7 >> attr(,"assign") >> [1] 0 1 >> >> >> Then how to set the contrast matrix? >> >> >> Thanks! >> >> -- >> xuan >> >> >> >> [[alternative HTML version deleted]] >> >> ______________________________**_________________ >> Bioconductor mailing list >> Bioconductor@r-project.org <mailto:bioconductor@r-**project.org<bioconductor@r-project.org> >> > >> >> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https: stat.ethz.ch="" mailman="" listinfo="" bioconductor=""> >> Search the archives: >> http://news.gmane.org/gmane.**science.biology.informatics.** >> conductor<http: news.gmane.org="" gmane.science.biology.informatics.c="" onductor=""> >> >> >> -- James W. MacDonald, M.S. >> Biostatistician >> University of Washington >> Environmental and Occupational Health Sciences >> 4225 Roosevelt Way NE, # 100 >> Seattle WA 98105-6099 >> >> >> >> >> -- >> xuan >> > > -- > James W. MacDonald, M.S. > Biostatistician > University of Washington > Environmental and Occupational Health Sciences > 4225 Roosevelt Way NE, # 100 > Seattle WA 98105-6099 > > -- Chunxuan Shao [[alternative HTML version deleted]]