Sorry forgot the sessionInfo! Added at the end now

--


From: Natasha Sahgal
Sent: 14 May 2013 12:29
To: bioconductor@r-project.org
Subject: DESeq multifactor design paired analysis

Dear List,

I have a dataset with 3 conditions (control and 2 treatments), 2 replicates each. Initially I carried out a normal analysis in DESeq - each treatment vs the control, but later was told  it needed to be a paired analysis.

Thus, just to be consistent I decided to use DESeq itself. However, I am unsure how to extract my comparisons of interest.
------
Code:
>design
  condition pair
1        Cont    1
2        Cont    3
3        Trt1    1
4        Trt1    3
5        Trt2    1
6        Trt2    3

>cds2 <- newCountDataSet(gct, design)
>cds2 <- estimateSizeFactors(cds2)

>sizeFactors(cds2)
     Cont_1     Cont_3     Trt1_1     Trt1_3      Trt2_1     Trt2_3
0.9373964 1.1328686 1.0097990 1.0596419 0.8562104 1.0645546

.cds2 <- estimateDispersions(cds2,"pooled-CR",modelFormula=count ~ pair + condition)

>fit1 = fitNbinomGLMs(cds2, count ~ pair + condition)
>fit0 = fitNbinomGLMs(cds2, count ~ pair)

>str(fit1)
'data.frame':   4765 obs. of  6 variables:
$ (Intercept): num  10.19 9.99 6.89 6.48 4.46 ...
$ pair3      : num  -0.3018 -0.3222 0.068 0.0468 0.1504 ...
$ conditionTrt1: num  0.17 0.142 -0.495 0.125 0.319 ...
$ conditionTrt2: num  0.1882 0.00112 -0.2704 -0.09412 0.10651 ...
$ deviance   : num  0.0742 0.1569 0.9072 2.0612 0.3702 ...
$ converged  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
-  attr(*, "df.residual")= num 2

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



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.


sessionInfo()
R version 2.15.2 (2012-10-26)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8
 [7] LC_PAPER=C                 LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] gdata_2.12.0       WriteXLS_2.3.0     DESeq_1.10.1       lattice_0.20-6
[5] locfit_1.5-8       Biobase_2.18.0     BiocGenerics_0.4.0

loaded via a namespace (and not attached):
[1] annotate_1.36.0      AnnotationDbi_1.20.2 DBI_0.2-5
 [4] genefilter_1.40.0    geneplotter_1.36.0   grid_2.15.2
 [7] gtools_2.7.0         IRanges_1.16.4       MASS_7.3-22
[10] parallel_2.15.2      RColorBrewer_1.0-5   RSQLite_0.11.2
[13] splines_2.15.2       stats4_2.15.2        survival_2.36-14
[16] tcltk_2.15.2         tools_2.15.2         XML_3.95-0.1
[19] xtable_1.7-0

Many Thanks,
Natasha


	[[alternative HTML version deleted]]

