hi Steve, I have fixed this bug in version 1.0.18 and 1.1.22, and I added a unit test for such a design formula in the devel version.
best, Mike On Wed, Jul 10, 2013 at 8:52 AM, Michael Love <[email protected]> wrote: > Hi Steve, > > Thanks for pointing out this bug in the renaming code. I will fix this as > soon as possible. > > Best, > > Mike > > On Jul 10, 2013 1:00 AM, "Steve Lianoglou" <[email protected]> wrote: >> >> Greetings DESeq2ers, >> >> I've been playing around with different ways to encode a 2x2 factorial >> design I'm working with: 4 cell types, 4 treatments and exploring >> what's cooking in there. >> >> I thought that the nested interaction formula used in the limma user's >> guide (section 8.5.3) would be an easy way to explore some obvious >> questions since differential expression from different treatments >> within each cell type (among other things) shake out quite cleanly >> from the results of running the wald test since these end up as >> columns in the design matrix, eg: >> >> ~ cell + cell:treatment >> >> as long as the control treatment is the first level of the treatment >> factor, we got columns for all cell:treatment effects, eg. >> cellA:treatment1, cellA:treatment2, ..., cellD:treatment4) >> >> Yay! >> >> ... almost ... >> >> The problem is that the "renaming column" stuff in the fitNbinomGLMs >> function seems to assume that we'll have columns for all main effects >> in the design matrix, so the following line triggers an error: >> >> modelMatrixNames[match(convertNames$from, modelMatrixNames)] <- >> convertNames$to >> >> since you will find elements like "treatment2, ..., treatment4" among >> the elements in `covertNames$from`. The call to match then returns NA >> for these, and *boom*. >> >> An easy fix would be to add the following line after the call to >> `renameModelMatrixColumns`: >> >> convertNames <- subset(convertNames, from %in% modelMatrixNames) >> >> So the "if (renameCols)" block (line 1066 in core.R) now becomes: >> >> if (renameCols) { >> names(modelMatrixNames) <- modelMatrixNames >> convertNames <- renameModelMatrixColumns(modelMatrixNames, >> >> as.data.frame(colData(object)), >> modelFormula) >> convertNames <- subset(convertNames, from %in% modelMatrixNames) >> modelMatrixNames[match(convertNames$from, modelMatrixNames)] <- >> convertNames$to >> } >> >> >> >> -- >> Steve Lianoglou >> Computational Biologist >> Bioinformatics and Computational Biology >> Genentech _______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/bioc-devel
