[Bioc-devel] renameModelMatrixColumns mishap + patch for DESeq2

[Michael Love]

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
<michaelisaiahlove at gmail.com> 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" <lianoglou.steve at gene.com> 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