hi Sridhar, On Wed, Jun 11, 2014 at 9:51 AM, Sridhar A Malkaram <smalkaram at="" wvstateu.edu=""> wrote: > Hi Michael, > > Thanks for explaining in detail. > Just for clarification, I think we should add the "~" before the > "genotype + time", right? Yes good catch. I should have written "~ genotype + time" > Does "~" stand for formula notation in R? Yes. > > dds <- DESeq(dds, test="LRT", reduced=genotype + time) > > > Also, when I do, > > dds <- nbinomWaldTest(dds, betaPrior=FALSE) > res<-results(dds, name="time2.gen2") > head(res) > > The first two lines of the output: > log2 fold change: time2.gen2 > Wald test p-value: time2.gen2 > ..... > > What two components is the log2 fold change between? I mean, how do I read time2.gen2 ? > How can I get the two expression values that are being compared. The older version of DESeq used to give these values in the results. > The time series model with interactions between time and a condition variable is a bit more complicated than the simple group comparisons. There are not simply two levels being compared. The way that interaction effects are encoded in model matrices in R is to write the two terms of which the interaction is composed next to each other. We use this convention as well. I'll try to explain the meaning of the model you have specified. The model fits a number of main effects: Firstly, we have the change in expression levels over time for the base level of the condition variable (so gen1). These effects are the fold change (or addition in the log2 space) of the time 2 over time1 (the intercept, or base level), time 3 over time1, time 4 over time 1. If you specify results(dds, name="time_2_vs_1") you would get a results table with a simple description of the comparison. But this is probably not of interest for you for testing, as it doesn't tell you anything about the effect of the experimental condition of interest, e.g. the genotype. Secondly, we have the fold change of genotype 2 over 1 for the first time period. This is also obtained with results(dds, name="gen_2_vs_1") Finally are the interaction terms. For example, "gen2.time2" is the fold change for genotype 2 at time 2 which is not explained by the genotype 2 main fold change and the time 2 main fold change multiplied together (or added in the log2 scale). If the log fold change of this intercept term equals 0, then we say there is no interaction. In other words, the effect of genotype 2 is the same at time 2 as it was at time 1, after accounting for the general effect of time on expression. Mike > -- With Kind Regards, > Sridhar Acharya > > > > > > > > On 6/10/2014 8:32 PM, Michael Love wrote: >> I left out one detail, below: >> >> On Tue, Jun 10, 2014 at 8:20 PM, Michael Love >> <michaelisaiahlove at="" gmail.com=""> wrote: >>> hi Sridhar, >>> >>> Let's keep the discussion on the mailing list in case the question is >>> relevant to others. >>> >>> After you have run: >>> >>> design(dds) <- ~ genotype + time + genotype:time >>> dds <- DESeq(dds, test="LRT", reduced=genotype + time) >>> res <- results(dds) >>> >>> The res object will contain the likelihood ratio test results, with >>> small p-values for genes which have a genotype effect which is >>> different than in the first time period. This tests all time periods >>> after the first time period. >>> >>> You also see: >>> >>> resultsNames(ddsLRT) >>> Intercept time2_vs_time1 time3_vs_time1 time4_vs_time1 >>> gen2_vs_gen1 time2.gen2 time3.gen2 time4.gen2 >>> >>> The three of these which might be interesting for your experiment are: >>> >> Before the results line below you should call: >> >> dds <- nbinomWaldTest(dds, betaPrior=FALSE) >> >>> results(dds, name="time2.gen2") >>> ...and same for time3.gen2, time4.gen2 >>> >>> which will return a results table with Wald tests of the additional >>> genotype effect in time 2 (additional beyond the genotype effect in >>> the first time period). This is similar to the first LRT results >>> above, except now we are asking for a different effect of genotype in >>> a specific time period, not in all time periods. >>> >>> The other coefficients are the main effect terms. Results tables for >>> these can also be built by using the 'name' argument to results(). >>> They are the intercept term, the effects of the different times over >>> the initial time, and the effect for genotype 2 over 1 in the first >>> time period. You don't want to use the contrast argument, which is for >>> other kinds of models. >>> >>> Mike >>> >>> On Tue, Jun 10, 2014 at 3:38 PM, Michael Love >>> <michaelisaiahlove at="" gmail.com=""> wrote: >>>> hi Sridhar, >>>> >>>> On Tue, Jun 10, 2014 at 3:05 PM, Sridhar A Malkaram >>>> <smalkaram at="" wvstateu.edu=""> wrote: >>>>> Hi, >>>>> >>>>> >>>>> I have been a user of DESeq and recently DESeq2 for my research work. >>>>> The latest DESeq2 seem to offer extensive differential testing options >>>>> suitable for various experimental designs. >>>>> >>>>> Recently I wanted to use DESeq for a differential gene expression >>>>> analysis between two plant genotypes across 4 different time points. >>>>> >>>>> I am basically a biologist and am finding hard to grasp the concepts of >>>>> testing results. I'd be very grateful if you could help me understand >>>>> some concepts (especially resultsNames) related to the DESeq2 package. >>>>> >>>>> >>>>> My experimental design is as below >>>>> >>>>> design<- ~ genotype + time + genotype:time >>>>> >>>>> There are two levels in genotype and 4 levels in time. >>>>> Basically I'd like to use binomLRT test to check if there is any >>>>> difference in gene expression between the genotypes across the time points. >>>>> >>>>> dds<-DESeq(dds) (dds is DESeq2 object obtained from, >>>>> dds<-DESeqDataSetFromMatrix(countData=counts, colData=coldata, >>>>> design=design) >>>>> >>>>> and I am using the reduced model for the liklihood test >>>>> >>>> Here is where things are getting confused. You have already run >>>> DESeq() using test="Wald". So it doesn't make sense at this point to >>>> instead perform a likelihood ratio test. In our vignette we explain >>>> this in the section on the LRT: "The likelihood ratio test can also be >>>> speci?ed using the test argument to DESeq, which substitutes >>>> nbinomWaldTest with nbinomLRT." >>>> >>>>> Is the model correct per my research question (is there a (time >>>>> influenced) difference between genotypes)? >>>>> >>>> Yes. If you want to find those genes which show a time influenced >>>> difference between genotypes, this is simply: >>>> >>>> dds <- DESeq(dds, test="LRT", reduced=genotype + time) >>>> res <- results(dds) >>>> >>>> You can then use heatmaps to inspect the patterns of gene expression >>>> for the differentially expressed genes. Visualization with heatmaps >>>> are covered in the vignette. >>>> >>>> If you have other more specific questions about how to generate >>>> results tables, I can answer them. With time series experiments, there >>>> are many possible combinations to test, but rather than going through >>>> all combinations, we recommend that users explore the results with >>>> heatmaps. >>>> >>>> Mike >>>> >>>>> Thanks, >>>>> Sridhar Acharya >>>>> >>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> _______________________________________________ >>>>> 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 >