Dear Dr Anders I have a few questions about applying deseq to a multifactor analysis that I am doing. I have three factors: three levels of genotype, three levels of "cell culture" and two time points. I have three replicates in each cell. The comparison(s) of most interest are changes in expression with genotype and days. I am also interested in knowing the effects of cell culture, but it's a technical issue. I have used HTSeq-count with the ENSEMBL v54 H.sapiens gtf to get counts from reads aligned with tophat. I have 54 libraries total. The first few rows of the conditions file look like: day culture genotype SL27450 7 33D6p63 NT SL27451 7 33D6p63 NT SL27452 7 33D6p63 NT SL27456 14 33D6p63 NT SL27457 14 33D6p63 NT SL27458 14 33D6p63 NT such that each set of factors is repeated three times. After normalization and dispersion estimation I did: fit0 <- fitNbinomGLMs(cds, count~genotype+day+culture) fit1 <- fitNbinomGLMs(cds, count~day+culture) then pvalsGLM = nbinomGLMTest( fit1, fit0 ) The result was entirely NaN or NAs. I'll put some output and the environment at the end of this email. I'm not sure where to start looking for the problem. At the moment, I am wondering how best to implement the pairwise comparison alluded to in the "fix me" in the manual I have. I am wondering: Does the formula in the fitNbinomGLMs() function work the same way as lm(), e.g., do I have to explicitly include interactions or will it calculate them for me? Is there a way to use relevel() to set the reference level? I have a wild-type genotype that would be natural for comparing two mutant genotypes. I don't see any references on the list, but has anybody tried to use normalized values from a CountDataSet outside of deseq? Lastly, I saw on the list while researching this issue that deseq2 is out. I'll give it try. Here is str() output from the fits (which seems OK to me) > str(fit0) 'data.frame': 37435 obs. of 8 variables: $ (Intercept) : num 9.01 1.8 7.69 7.06 6.63 ... $ genotypeNT : num -0.0666 0.0223 -0.08 -0.1554 -0.1967 ... $ genotypeWT : num -0.1665 -0.6294 -0.0716 -0.0908 -0.0754 ... $ day7 : num -0.0849 -0.4425 -0.0596 -0.1809 -0.2374 ... $ culture33D6p63: num -0.153 -1.917 0.119 0.715 0.518 ... $ culture34D6p62: num -0.2925 -0.3422 0.0181 0.242 0.0233 ... $ deviance : num 46.1 44.6 51.8 21.1 60.2 ... $ converged : logi TRUE TRUE TRUE TRUE TRUE TRUE ... - attr(*, "df.residual")= num 48 > str(fit1) 'data.frame': 37435 obs. of 6 variables: $ (Intercept) : num 8.94 1.6 7.64 6.99 6.54 ... $ day7 : num -0.0826 -0.4417 -0.0589 -0.1793 -0.2361 ... $ culture33D6p63: num -0.163 -1.899 0.114 0.708 0.507 ... $ culture34D6p62: num -0.3033 -0.3009 0.0161 0.2411 0.0191 ... $ deviance : num 49.7 46.9 52.7 23.9 64.1 ... $ converged : logi TRUE TRUE TRUE TRUE TRUE TRUE ... - attr(*, "df.residual")= num 50 > sessionInfo() R version 2.15.0 (2012-03-30) Platform: x86_64-unknown-linux-gnu (64-bit) 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=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid splines stats graphics grDevices utils datasets [8] methods base other attached packages: [1] gplots_2.11.0 MASS_7.3-23 KernSmooth_2.23-9 caTools_1.14 [5] gdata_2.12.0 gtools_2.7.0 RColorBrewer_1.0-5 multcomp_1.2-17 [9] survival_2.37-4 mvtnorm_0.9-9994 heatmap.plus_1.3 DESeq_1.10.1 [13] lattice_0.20-6 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.6 bitops_1.0-5 [4] DBI_0.2-5 genefilter_1.40.0 geneplotter_1.36.0 [7] IRanges_1.16.6 parallel_2.15.0 RSQLite_0.11.2 [10] stats4_2.15.0 tools_2.15.0 XML_3.95-0.2 [13] xtable_1.7-1 Michael Muratet, Ph.D. Senior Scientist HudsonAlpha Institute for Biotechnology mmuratet at hudsonalpha.org (256) 327-0473 (p) (256) 327-0966 (f) Room 4005 601 Genome Way Huntsville, Alabama 35806

Hi Michel On 22/03/13 20:52, Michael Muratet wrote: > fit0 <- fitNbinomGLMs(cds, count~genotype+day+culture) > fit1 <- fitNbinomGLMs(cds, count~day+culture) > then > pvalsGLM = nbinomGLMTest( fit1, fit0 ) > > The result was entirely NaN or NAs. This is easy. :-) It is simply because you supplied the fits in the wrong oder. The fit with 'genotype' is the full model and should come as first argument to nbinomGLMTest, and the reduced model without 'genotype' is second. So, just swap fit0 and fit1. > Does the formula in the fitNbinomGLMs() function work the same way as > lm(), e.g., do I have to explicitly include interactions or will it > calculate them for me? Yes, both fitNbinomGLMs and lm internally call model.matrix, which, as you know, adds an intercept automatically > Is there a way to use relevel() to set the reference level? I have a > wild-type genotype that would be natural for comparing two mutant > genotypes. Yes, you can use 'relevel' in the same way as you would with 'lm'. And to answer your question about the 'fixme'. There are two possibilities: You can compare a single genotype with wildtype by simply subsetting the CountDataSet to exclude the samples for the othe rgenotype before doing the fits. Or you can use our new DESeq2, which now offers a Wald test instead of a likelihood ratio test and can so give you p values for both genotype coefficient in one go. > I don't see any references on the list, but has anybody tried to use > normalized values from a CountDataSet outside of deseq? Depends on what you want to do with them. > Lastly, I saw on the list while researching this issue that deseq2 is > out. I'll give it try. Sure, do so. And tell us what you think. Simon

Hi [Please keep conversations on the list.] On 25/03/13 19:30, Michael Muratet wrote: > I thought I would use them in R to manually calculate the GLM where I > could have access to the model. Or, maybe there is a way to get that > in deseq2? What kind of "access" do you need? You can specify a formula (or even a model matrix, if you prefer), and you get all coefficients ande deviances back. Of course, you can call glm.fit or statmod::glmnb.fit directly. Then, however, you should _not_ use the normalized counts but rather specify the log size factors as offsets. > I tried to install deseq2 and discovered I need R 2.16. CRAN only > goes up to 15.3, but they have beta r003, is this what you used? No, I used R 2.15. See the explanations on the Bioc web page on how R and Bioc versions are tied together, and how devel packages are handled. Simon