resultsNames returning only "Intercept" after setting up paired multi-factor design in DESeq2
2
0
Entering edit mode
dh • 0
@dh-9632
Last seen 8.2 years ago

Hi - I have an experiment with paired samples (normal vs cancer) from mice in triplicate of two different genotypes (WT vs KO). It is basically the same situation as described in DESeq2's documentation in section 3.12.1. 

So I have used the edgeR trick of nesting individual mice in each genotype, as recommended. However, when I run the resultsNames, only "Intercept" is returned - so I am unable to test for any of the interactions I am interested in. Here is the code that I run:

> dds = DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=countFilesPath, design=~1)

> colData(dds)

DataFrame with 12 rows and 3 columns

                grp    mouse      cnd

           <factor> <factor> <factor>

D207-N2-KO       KO     D207   normal

D207-P1-KO       KO     D207   cancer

D521-N3-KO       KO     D521   normal

D521-P2-KO       KO     D521   cancer

D540-N2-KO       KO     D540   normal

...             ...      ...      ...

D627-P1          WT     D627   cancer

D686-N2          WT     D686   normal

D686-P5          WT     D686   cancer

D868-N3          WT     D868   normal

D868-P2          WT     D868   cancer

> dds$ind.n = factor(rep(rep(1:3,each=2),2))

> colData(dds)

DataFrame with 12 rows and 4 columns

                grp    mouse      cnd    ind.n

           <factor> <factor> <factor> <factor>

D207-N2-KO       KO     D207   normal        1

D207-P1-KO       KO     D207   cancer        1

D521-N3-KO       KO     D521   normal        2

D521-P2-KO       KO     D521   cancer        2

D540-N2-KO       KO     D540   normal        3

...             ...      ...      ...      ...

D627-P1          WT     D627   cancer        1

D686-N2          WT     D686   normal        2

D686-P5          WT     D686   cancer        2

D868-N3          WT     D868   normal        3

D868-P2          WT     D868   cancer        3

> dds$cnd = relevel(dds$cnd, ref="normal") # set reference comparison to normal

> dds$grp = relevel(dds$grp, ref="WT") # set reference comparison to WT

> design(dds) ~ grp + grp:ind.n + grp:cnd

design(dds) ~ grp + grp:ind.n + grp:cnd

> model.matrix(~ grp + grp:ind.n + grp:cnd, colData(dds)) # inspect model matrix

           (Intercept) grpKO grpWT:ind.n2 grpKO:ind.n2 grpWT:ind.n3

D207-N2-KO           1     1            0            0            0

D207-P1-KO           1     1            0            0            0

D521-N3-KO           1     1            0            1            0

D521-P2-KO           1     1            0            1            0

D540-N2-KO           1     1            0            0            0

D540-P1-KO           1     1            0            0            0

D627-N3              1     0            0            0            0

D627-P1              1     0            0            0            0

D686-N2              1     0            1            0            0

D686-P5              1     0            1            0            0

D868-N3              1     0            0            0            1

D868-P2              1     0            0            0            1

           grpKO:ind.n3 grpWT:cndcancer grpKO:cndcancer

D207-N2-KO            0               0               0

D207-P1-KO            0               0               1

D521-N3-KO            0               0               0

D521-P2-KO            0               0               1

D540-N2-KO            1               0               0

D540-P1-KO            1               0               1

D627-N3               0               0               0

D627-P1               0               1               0

D686-N2               0               0               0

D686-P5               0               1               0

D868-N3               0               0               0

D868-P2               0               1               0

attr(,"assign")

[1] 0 1 2 2 2 2 3 3

attr(,"contrasts")

attr(,"contrasts")$grp

[1] "contr.treatment"


attr(,"contrasts")$ind.n

[1] "contr.treatment"


attr(,"contrasts")$cnd

[1] "contr.treatment"


> dds = DESeq(dds, modelMatrixType="standard")

estimating size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

-- replacing outliers and refitting for 196 genes

-- DESeq argument 'minReplicatesForReplace' = 7 

-- original counts are preserved in counts(dds)

estimating dispersions

fitting model and testing

> resultsNames(dds)

[1] "Intercept"

> res = results(dds, name="grpWT.cndcancer")

Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) : 

  subscript contains invalid names

 

I'm baffled as to what I'm doing wrong - any help? 

 

Thanks!!

 

R version 3.2.1 (2015-06-18)

Platform: x86_64-unknown-linux-gnu (64-bit)

Running under: Ubuntu quantal (12.10)


locale:

 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              

 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    

 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   

 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 

 [9] LC_ADDRESS=C               LC_TELEPHONE=C            

[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       


attached base packages:

[1] parallel  stats4    stats     graphics  grDevices utils     datasets 

[8] methods   base     


other attached packages:

 [1] dplyr_0.4.3                DESeq2_1.10.1             

 [3] RcppArmadillo_0.6.500.4.0  Rcpp_0.12.3               

 [5] SummarizedExperiment_1.0.2 Biobase_2.30.0            

 [7] GenomicRanges_1.22.4       GenomeInfoDb_1.6.3        

 [9] IRanges_2.4.6              S4Vectors_0.8.11          

[11] BiocGenerics_0.16.1        biomaRt_2.26.1            

[13] BiocInstaller_1.20.1      


loaded via a namespace (and not attached):

 [1] RColorBrewer_1.1-2   futile.logger_1.4.1  plyr_1.8.3          

 [4] XVector_0.10.0       bitops_1.0-6         futile.options_1.0.0

 [7] tools_3.2.1          zlibbioc_1.16.0      rpart_4.1-10        

[10] annotate_1.48.0      RSQLite_1.0.0        gtable_0.1.2        

[13] lattice_0.20-33      DBI_0.3.1            gridExtra_2.0.0     

[16] genefilter_1.52.1    cluster_2.0.3        locfit_1.5-9.1      

[19] nnet_7.3-11          grid_3.2.1           R6_2.1.2            

[22] AnnotationDbi_1.32.3 XML_3.98-1.3         survival_2.38-3     

[25] BiocParallel_1.4.3   foreign_0.8-66       latticeExtra_0.6-26 

[28] Formula_1.2-1        magrittr_1.5         geneplotter_1.48.0  

[31] ggplot2_2.0.0        lambda.r_1.1.7       Hmisc_3.17-1        

[34] scales_0.3.0         splines_3.2.1        assertthat_0.1      

[37] xtable_1.8-0         colorspace_1.2-6     acepack_1.3-3.3     

[40] RCurl_1.95-4.7       munsell_0.4.2       

 

DESeq2 • 2.0k views
ADD COMMENT
1
Entering edit mode
@mikelove
Last seen 8 minutes ago
United States

hi,

This looks like a simple fix, you need to actually assign the design to the dds object, using <-

design(dds) <- ...

You can check to see if you successfully assigned by typing design(dds) and hitting Return. It should print the design formula.

ADD COMMENT
0
Entering edit mode
dh • 0
@dh-9632
Last seen 8.2 years ago

Thanks!! That fixes it - I knew I must have been missing something silly like that.. 

ADD COMMENT

Login before adding your answer.

Traffic: 798 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6