hi, I cannot reproduce the bug with the data and code you sent me using the devel branch. DESeq() runs without error for me. It could be another library you have loaded other than phyloseq and DESeq2? My session info: R Under development (unstable) (2014-06-05 r65862) Platform: x86_64-apple-darwin12.5.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods [8] base other attached packages: [1] DESeq2_1.5.20 RcppArmadillo_0.4.300.0 Rcpp_0.11.1 [4] GenomicRanges_1.17.17 GenomeInfoDb_1.1.6 IRanges_1.99.15 [7] S4Vectors_0.0.7 BiocGenerics_0.11.2 phyloseq_1.9.9 [10] Defaults_1.1-1 devtools_1.5 knitr_1.6 [13] BiocInstaller_1.15.5 loaded via a namespace (and not attached): [1] ade4_1.6-2 annotate_1.43.4 AnnotationDbi_1.27.7 [4] ape_3.1-2 Biobase_2.25.0 biom_0.3.12 [7] Biostrings_2.33.9 cluster_1.15.2 codetools_0.2-8 [10] colorspace_1.2-4 compiler_3.2.0 data.table_1.9.2 [13] DBI_0.2-7 digest_0.6.4 evaluate_0.5.5 [16] foreach_1.4.2 formatR_0.10 genefilter_1.47.5 [19] geneplotter_1.43.0 ggplot2_1.0.0 grid_3.2.0 [22] gtable_0.1.2 httr_0.3 igraph_0.7.1 [25] iterators_1.0.7 lattice_0.20-29 locfit_1.5-9.1 [28] MASS_7.3-33 Matrix_1.1-3 memoise_0.2.1 [31] multtest_2.21.0 munsell_0.4.2 nlme_3.1-117 [34] permute_0.8-3 plyr_1.8.1 proto_0.3-10 [37] RColorBrewer_1.0-5 RCurl_1.95-4.1 reshape2_1.4 [40] RJSONIO_1.2-0.2 RSQLite_0.11.4 scales_0.2.4 [43] splines_3.2.0 stats4_3.2.0 stringr_0.6.2 [46] survival_2.37-7 tools_3.2.0 vegan_2.0-10 [49] whisker_0.3-2 XML_3.98-1.1 xtable_1.7-3 [52] XVector_0.5.6 zlibbioc_1.11.1 On Thu, Jun 26, 2014 at 4:44 PM, Michael Love <michaelisaiahlove@gmail.com> wrote: > hi, > > I'm not sure where that error is coming from. I don't know what 'NSBS > object' is referring to and couldn't find any relevant information from > google either. > > phyloseq has a DESeq() call in the vignette, and it's passing check in > devel, so not sure what's going on. > > If you can make a reproducible example dataset and email me privately, > I'll take a look and respond on the list. > > Mike > > > On Thu, Jun 26, 2014 at 4:32 PM, Shucong Li <lishucongnn@gmail.com> wrote: > >> Hello Michael, >> >> I tried to use the latest version of DESeq2 (1.5.20) and phyloseq (1.9.9) >> to analyze an illumina sequencing data set. However, it didn't work and I >> got the error: >> >> > library(DESeq2) >> >> > phyloseq.obj <- seq_all >> > dds <- phyloseq_to_deseq2(phyloseq.obj,design= ~ Treatment*Day) >> converting counts to integer mode >> > dds <- DESeq(dds, test = "Wald", fitType = "parametric",betaPrior=FALSE) >> estimating size factors >> estimating dispersions >> gene-wise dispersion estimates >> Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) : >> subscript is a NSBS object that is incompatible >> with the current subsetting operation >> >> >> I then tried use DESeq2(1.4.5) and phyloseq (1.8.2) to run the same >> codes, it worked without any issue. I also tried the combination of >> DESeq2 (1.4.5) + phyloseq (1.9.9), DESeq2 (1.5.20) and phyloseq (1.8.2), >> they didn't work either and I go the same message (Error in NSBS(i, x, >> exact = exact, upperBoundIsStrict = !allow.append) ). >> >> >> >> Here is the output of sessionInfo() >> >> R version 3.1.0 (2014-04-10) >> Platform: x86_64-apple-darwin13.1.0 (64-bit) >> >> locale: >> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8 >> >> attached base packages: >> [1] parallel grid splines stats graphics grDevices utils >> datasets methods base >> >> other attached packages: >> [1] DESeq2_1.5.20 ecodist_1.2.9 >> Biostrings_2.32.0 doParallel_1.0.8 foreach_1.4.2 >> [6] iterators_1.0.7 metagenomeSeq_1.6.0 gplots_2.14.0 >> limma_3.20.1 Biobase_2.24.0 >> [11] GenomicRanges_1.16.3 GenomeInfoDb_1.0.2 >> RcppArmadillo_0.4.300.8.0 Rcpp_0.11.2 XVector_0.4.0 >> [16] IRanges_1.22.6 BiocGenerics_0.11.2 locfit_1.5-9.1 >> phangorn_1.99-7 genefilter_1.46.1 >> [21] adephylo_1.1-6 scatterplot3d_0.3-35 analogue_0.12-0 >> rgl_0.93.996 princurve_1.1-12 >> [26] labdsv_1.6-1 mgcv_1.7-29 >> indicspecies_1.7.1 biom_0.3.13 ggplot2_1.0.0 >> [31] reshape2_1.4 plyr_1.8.1 phyloseq_1.9.9 >> pamr_1.54.1 cluster_1.15.2 >> [36] survival_2.37-7 vegan_2.0-10 lattice_0.20-29 >> permute_0.8-3 RColorBrewer_1.0-5 >> [41] matrixStats_0.10.0 MASS_7.3-33 ape_3.1-2 >> ade4_1.6-2 nlme_3.1-117 >> >> loaded via a namespace (and not attached): >> [1] adegenet_1.4-2 annotate_1.42.0 AnnotationDbi_1.26.0 >> bitops_1.0-6 brglm_0.5-9 caTools_1.17 >> codetools_0.2-8 >> [8] colorspace_1.2-4 data.table_1.9.2 DBI_0.2-7 >> digest_0.6.4 fastmatch_1.0-4 gdata_2.13.3 >> geneplotter_1.42.0 >> [15] gtable_0.1.2 gtools_3.4.1 htmltools_0.2.4 >> httpuv_1.3.0 igraph_0.7.0 KernSmooth_2.23-12 Matrix_1.1-4 >> [22] multtest_2.20.0 munsell_0.4.2 phylobase_0.6.8 >> proto_0.3-10 R.methodsS3_1.6.1 RJSONIO_1.2-0.2 >> RSQLite_0.11.4 >> [29] S4Vectors_0.0.8 scales_0.2.4 shiny_0.10.0 >> stats4_3.1.0 stringr_0.6.2 tools_3.1.0 XML_3.98-1.1 >> [36] xtable_1.7-3 zlibbioc_1.10.0 >> >> >> Shucong Li >> lishucongnn@gmail.com >> >> >> >> On May 17, 2014, at 11:32 PM, Shucong Li <lishucongnn@gmail.com> wrote: >> >> > Hi Mike, >> > >> > Thank you. With your help, I finally figured it out, sort of :-) >> > >> > Enjoy the rest of your weekend. >> > >> > Shawn >> > On May 17, 2014, at 11:07 AM, Michael Love <michaelisaiahlove@gmail.com> >> wrote: >> > >> >> hi Shawn, >> >> >> >> Let's keep all replies on the Bioconductor mailing list, so other >> >> users can follow the thread. >> >> >> >> >> >> On Sat, May 17, 2014 at 11:32 AM, Shucong Li <lishucongnn@gmail.com> >> wrote: >> >>> Hello Mike, >> >>> >> >>> Thank you for reply on weekend. >> >>> >> >>> Sorry for not make a part of my question properly. What I really >> meant is to compare : >> >>> >> >>> "Factor A level 1" vs "Factor A level 2" within "Fact B level 1" >> >>> "Factor A level 1" vs "Factor A level 2" within "Fact B level 2" >> >> >> >> I don't understand what you mean by these comparisons, specifically >> >> the 'within' part. >> >> >> >> There are four terms in this model. I will write in the log2 scale, so >> >> terms are additive (this is multiplication on the scale of counts): >> >> >> >> log2(counts) = intercept + treatmenttreat + dayB + dayB.treatmenttreat >> >> >> >> the treatmenttreat term is the difference between treatment and >> >> control for day A >> >> >> >> the dayB term is the difference between day B over day A for the >> control group >> >> >> >> the final term accounts for the case that the (day B and treatment) >> >> group is not merely the sum of the previous two terms. >> >> >> >> If you want to compare treatment vs control for day B, that term is >> >> the treatmenttreat term plus the dayB.treatmenttreat term. This can be >> >> pulled out easily with a numeric contrast (see the help for ?results >> >> for more details). >> >> >> >> results(dds, contrast=c(0,1,0,1)) >> >> >> >> where the numbers correspond to the order in resultsNames(dds) >> >> >> >> I'd recommend looking over some introductory material to interaction >> >> models (e.g. http://en.wikipedia.org/wiki/Interaction_(statistics) ). >> >> >> >> Mike >> >> >> >>> >> >>> Is it possible to do so? I checked all vignettes and manual of >> DESeq2 and still could not figure this out. >> >>> >> >>> Shawn >> >>> >> >>> On May 17, 2014, at 9:58 AM, Michael Love < >> michaelisaiahlove@gmail.com> wrote: >> >>> >> >>>> hi Shawn, >> >>>> >> >>>> >> >>>> On Sat, May 17, 2014 at 2:00 AM, Shawn [guest] < >> guest@bioconductor.org> wrote: >> >>>>> Hello Mike, >> >>>>> >> >>>>> Please allow me ask a basic question. What does 'log2FoldChange' >> in the results of DESeq2 analysis really mean for the interaction of a >> two-factor two-level design? >> >>>> >> >>>> The meaning is the same as for other linear models. The interaction >> >>>> term for the generalized linear model tests if the effect of both day >> >>>> and treatment is not simply multiplicative (or additive in the log2 >> >>>> space). If this term is significantly non-zero (the default test in >> >>>> DESeq2), then being in the group: (day=B and treatment=treat) is not >> >>>> simply the product of the day=B fold change and the treatment=treat >> >>>> fold change. >> >>>> >> >>>> >> >>>>> Is it possible to compare 'Factor A level 1' to ' Factor A level 2' >> or other similar comparison? >> >>>> >> >>>> For example, the day B over A effect: >> >>>> >> >>>> results(dds, contrast=c("day","B","A")) >> >>>> >> >>>> or equivalently >> >>>> >> >>>> results(dds, name="day_B_vs_A") >> >>>> >> >>>> Check the help for ?results for more examples. >> >>>> >> >>>> Mike >> >>>> >> >>>>> >> >>>>> Here are the part of codes I used: >> >>>>> >> >>>>> >> >>>>> dds <- phyloseq_to_deseq2(phyloseq.obj, design=~ Treatment*Day) >> >>>>> >> >>>>> colData(dds)$Treatment<- >> factor(colData(dds)$Treatment,levels=c("Control","Treat")); >> >>>>> colData(dds)$Day<- factor(colData(dds)$Day,levels=c("A","B")) >> >>>>> >> >>>>> dds$Treatment<- relevel(dds$Treatment, "Control") >> >>>>> dds$Day<- relevel(dds$Day, "A") >> >>>>> >> >>>>> dds <- DESeq(dds, fitType="local",betaPrior=FALSE) >> >>>>> >> >>>>> Shawn >> >>>>> >> >>>>> >> >>>>> -- output of sessionInfo(): >> >>>>> >> >>>>> sessionInfo() >> >>>>> R version 3.1.0 (2014-04-10) >> >>>>> Platform: x86_64-apple-darwin13.1.0 (64-bit) >> >>>>> >> >>>>> locale: >> >>>>> [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8 >> >>>>> >> >>>>> attached base packages: >> >>>>> [1] parallel grid splines stats graphics grDevices >> utils datasets methods base >> >>>>> >> >>>>> other attached packages: >> >>>>> [1] ecodist_1.2.9 Biostrings_2.32.0 >> doParallel_1.0.8 foreach_1.4.2 >> >>>>> [5] iterators_1.0.7 metagenomeSeq_1.6.0 gplots_2.13.0 >> limma_3.20.1 >> >>>>> [9] Biobase_2.24.0 DESeq2_1.4.5 >> GenomicRanges_1.16.3 GenomeInfoDb_1.0.2 >> >>>>> [13] RcppArmadillo_0.4.300.0 Rcpp_0.11.1 XVector_0.4.0 >> IRanges_1.22.6 >> >>>>> [17] BiocGenerics_0.10.0 locfit_1.5-9.1 >> phangorn_1.99-7 genefilter_1.46.1 >> >>>>> [21] adephylo_1.1-6 scatterplot3d_0.3-35 >> analogue_0.12-0 rgl_0.93.996 >> >>>>> [25] princurve_1.1-12 labdsv_1.6-1 mgcv_1.7-29 >> indicspecies_1.7.1 >> >>>>> [29] biom_0.3.13 ggplot2_0.9.3.1 plyr_1.8.1 >> phyloseq_1.9.2 >> >>>>> [33] pamr_1.54.1 cluster_1.15.2 >> survival_2.37-7 vegan_2.0-10 >> >>>>> [37] lattice_0.20-29 permute_0.8-3 >> RColorBrewer_1.0-5 matrixStats_0.8.14 >> >>>>> [41] MASS_7.3-33 ape_3.1-1 ade4_1.6-2 >> nlme_3.1-117 >> >>>>> >> >>>>> loaded via a namespace (and not attached): >> >>>>> [1] adegenet_1.4-1 annotate_1.42.0 AnnotationDbi_1.26.0 >> bitops_1.0-6 >> >>>>> [5] brglm_0.5-9 caTools_1.17 codetools_0.2-8 >> colorspace_1.2-4 >> >>>>> [9] data.table_1.9.2 DBI_0.2-7 digest_0.6.4 >> fastmatch_1.0-4 >> >>>>> [13] gdata_2.13.3 geneplotter_1.42.0 gtable_0.1.2 >> gtools_3.4.0 >> >>>>> [17] httpuv_1.3.0 igraph_0.7.0 KernSmooth_2.23-12 >> Matrix_1.1-3 >> >>>>> [21] multtest_2.20.0 munsell_0.4.2 phylobase_0.6.8 >> proto_0.3-10 >> >>>>> [25] R.methodsS3_1.6.1 reshape2_1.4 RJSONIO_1.2-0.2 >> RSQLite_0.11.4 >> >>>>> [29] scales_0.2.4 shiny_0.9.1 stats4_3.1.0 >> stringr_0.6.2 >> >>>>> [33] tools_3.1.0 XML_3.98-1.1 xtable_1.7-3 >> zlibbioc_1.10.0 >> >>>>> >> >>>>> >> >>>>> -- >> >>>>> Sent via the guest posting facility at bioconductor.org. >> >>> >> > >> >> > [[alternative HTML version deleted]]