Entering edit mode
Hi Mike,
I guess you are right. I use R 3.1.0. I notice that R3.1.1 will come
out in next a few days. I will then upgrade R and try it again. It the
issue still exists, I will let you know.
Thanks.
Shucong Li
lishucongnn@gmail.com
On Jun 26, 2014, at 4:45 PM, Michael Love
<michaelisaiahlove@gmail.com> wrote:
> 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]]
