Hi David, On 2/26/2014 11:04 AM, David Mori?a Soler wrote: > Hi Jim, > > thank you very much for your help! If I'm right, the solution you propose > > > contrast <- makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, > Treat2diff = Treat2.6h-Treat2.0h, levels = design) > > gives the difference in expression between times inside each treatment > group, while I'm more interested in the 'conventional' treatment main > effect, as you suggest, a main effect for treatment, after > controlling for time. Is there a way to achieve that using limma? Of course! You just use model.matrix(): > treat <- factor(rep(paste0("Treat", 1:2), each = 2, times = 6)) > treat [1] Treat1 Treat1 Treat2 Treat2 Treat1 Treat1 Treat2 Treat2 Treat1 Treat1 [11] Treat2 Treat2 Treat1 Treat1 Treat2 Treat2 Treat1 Treat1 Treat2 Treat2 [21] Treat1 Treat1 Treat2 Treat2 Levels: Treat1 Treat2 > time <- factor(rep(paste0("Hr", c(0,6)), times = 12)) > time [1] Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 Hr6 Hr0 [20] Hr6 Hr0 Hr6 Hr0 Hr6 Levels: Hr0 Hr6 > design <- model.matrix(~treat + time) > colnames(design) <- gsub("treat|time", "", colnames(design)) ## this step just necessary because R likes to append things to the column names... > head(design) (Intercept) Treat2 Hr6 1 1 0 0 2 1 0 1 3 1 1 0 4 1 1 1 5 1 0 0 6 1 0 1 > In this parameterization Treat2 is implicitly the difference between Treat2 and Treat1, and Hr6 is implicitly the difference between the 6 hour and 0 hour times so you don't use a contrasts.fit() step: fit <- lmFit(eset, design) fit <- eBayes(fit) topTable(fit, 2) ## treatment main effect topTable(fit, 3) ## time main effect But note that these main effects are nonsensical in the presence of a significant interaction term. So you could hypothetically do > design <- model.matrix(~treat*time) > colnames(design) <- gsub("treat|time", "", colnames(design)) > head(design) (Intercept) Treat2 Hr6 Treat2:Hr6 1 1 0 0 0 2 1 0 1 0 3 1 1 0 0 4 1 1 1 1 5 1 0 0 0 6 1 0 1 0 And the Treat2:Hr6 coefficient is the interaction term. You would then want to first find the genes with significant interaction, and then ignore those genes when testing for significant main effects. Best, Jim > > Thanks! > > David > > > On 26/02/14 16:55, James W. MacDonald wrote: >> Hi David, >> >> On Wednesday, February 26, 2014 8:09:21 AM, David Mori?a Soler wrote: >>> Dear Bioconductor users, >>> >>> I just began to use limma package for the analysis of microarray data. >>> I want to include in the analysis two factors (treatment and time) and >>> the interaction between them. My design matrix looks like >>> > head(design) >>> Treat1.0h Treat1.6h Treat2.0h Treat2.6h >>> 1 0 0 1 0 >>> 2 0 0 0 1 >>> 3 1 0 0 0 >>> 4 0 1 0 0 >>> 5 0 0 1 0 >>> 6 0 0 0 1 >>> >>> Then, if I run >>> > fit <- >>> lmFit(data.norm,design,block=pData(data)$Subject,correlation=corfi t$consensus), >>> >>> >>> >>> I have this output for the coefficients: >>> > head(fit$coefficients) >>> Treat1.0h Treat1.6h Treat2.0h Treat2.6h >>> KCNE4 5.020670 4.981786 5.038805 4.924326 >>> IRG1 6.119265 6.015868 6.095171 6.027943 >>> SNAR-G2 8.242385 8.186429 8.144942 8.230391 >>> MBNL3 10.438644 10.417312 10.417042 10.358303 >>> HOXC4 8.985834 8.854698 8.969801 8.939682 >>> ENST00000319813 3.913602 4.102653 4.000681 3.960431 >>> >>> How can I obtain for example the effect of the treatment, using >>> contrasts and eBayes? >> >> It depends on what you mean by 'the effect of the treatment'. Are you >> looking for a main effect for treatment, after controlling for time >> (e.g., a conventional main effect)? Or are you looking for the >> difference in treatment at each time separately? >> >> The easiest way to create a contrasts matrix is to use >> makeContrasts(), which is as simple as deciding what comparisons you >> want to make. So for instance, let's say you want to know the >> difference in expression between the treatments at each time: >> >> contrast <- makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, >> Treat2diff = Treat2.6h-Treat2.0h, levels = design) >> >> If you want the interaction too, it's easy to do that as well >> >> contrast <- makeContrasts(Treat1diff = Treat1.6h-Treat1.0h, >> Treat2diff = Treat2.6h-Treat2.0h, Interaction = >> Treat1.6h-Treat1.0h-Treat2.6h+Treat2.0h, levels = design) >> >> Best, >> >> Jim >> >> >>> >>> Thank you very much, >>> >>> David >>> >>> _______________________________________________ >>> Bioconductor mailing list >>> Bioconductor at r-project.org >>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>> Search the archives: >>> http://news.gmane.org/gmane.science.biology.informatics.conductor >> >> -- >> James W. MacDonald, M.S. >> Biostatistician >> University of Washington >> Environmental and Occupational Health Sciences >> 4225 Roosevelt Way NE, # 100 >> Seattle WA 98105-6099 > > -- James W. MacDonald, M.S. Biostatistician University of Washington Environmental and Occupational Health Sciences 4225 Roosevelt Way NE, # 100 Seattle WA 98105-6099