Entering edit mode
dh
•
0
@dh-9632
Last seen 9.9 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
