Dear Michael,

Thank you for your quick response, I shall try what you say after upgrading R and relevant packages.

Many Thanks,
Natasha

--


From: Michael Love [mailto:michaelisaiahlove@gmail.com]
Sent: 14 May 2013 14:39
To: Natasha Sahgal
Cc: bioconductor@r-project.org
Subject: Re: [BioC] DESeq multifactor design paired analysis

Hi Natasha,

On Tue, May 14, 2013 at 1:30 PM, Natasha Sahgal <nsahgal@well.ox.ac.uk<mailto:nsahgal@well.ox.ac.uk>> wrote:

>head(fit1)
         (Intercept)       pair3 conditionTrt1  conditionTrt2  deviance converged
gene0      10.188336 -0.30184662   0.1698471  0.188197344 0.0741656      TRUE
gene1       9.992372 -0.32221027   0.1415406  0.001119304 0.1568869      TRUE
gene10      6.893693  0.06795915  -0.4954365 -0.270402107 0.9072335      TRUE
gene100     6.477847  0.04684995   0.1250743 -0.094116036 2.0611685      TRUE
gene1002    4.463446  0.15038277   0.3186187  0.106512332 0.3701619      TRUE
gene1003    4.090079 -0.05247162  -0.3560189 -0.054765642 1.9631768      TRUE

>pvalsGLM = nbinomGLMTest( fit1, fit0 )
>padjGLM = p.adjust(pvalsGLM, method="BH" )
---------

Based on the above how do I extract:

1)       Trt1 vs Cont? - is this the 3rd column of fit1, depicted as log2 FoldChange?



2)       Trt2 vs Cont? - is this the 4th  column of fit1, depicted as log2 FoldChange?

Yes, you are right. The 3rd and 4th columns of fit1 are the log2 fold changes for trt1 vs control and trt2 vs control.



3)      What would the p-values for each comparison be? As pvalsGLM and padjGLM gives only one set, and I don't think it will be the same for both comparisons.



The p-values and adjusted p-values here are from a likelihood ratio test of a model including the condition and pair variables vs a model which only includes the pair variable. However, we can provide tests of significance of each treatment from a model including both now in the DESeq2 package. We tried to simplify the interface, so it shouldn't take too much time to set up the same analysis. Note that the methods have been updated from DESeq, and these changes are described in the vignette in Appendix A.

To use DESeq2, you would create a DESeqDataSet and extract the results (log2 fold changes, p-values and adjusted p-values) for each treatment, with something like:

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ pair + condition)
dds <- DESeq(dds)
resTrt1 <- results(dds, "condition_Trt1_vs_Cont")
resTrt2 <- results(dds, "condition_Trt2_vs_Cont")

Mike

	[[alternative HTML version deleted]]

