[BioC] Using DESeq2: experimental design and extracting results

Sridhar A Malkaram smalkaram at wvstateu.edu
Wed Jun 11 15:51:02 CEST 2014


Hi Michael,

Thanks for explaining in detail.
Just for clarification, I think we should add the "~" before the
"genotype + time", right?
Does "~" stand for formula notation in R?

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.

-- 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
>>> specified 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



More information about the Bioconductor mailing list