Package interactions in R

Question

[mikael.lenz.strube]

Hi,

having a weird problem. I wrote a package which me and my team use to analyse NGS data, in particular to extract and sanitize it from a complicated workflow. I then thought i would give DESeq2 a go, which i could not get to work for the life of me. For some reason I figured my package might be the reason, so i saved the sanitized data as an .Rdata file and tried to do the DESeq2 on the data without my own package loaded, which worked.

Example

library(MyPackage) 

countsTable <- aBunchOfMyFunctions
Group <- factor(SomeGroupNames)
coldata1 <- data.frame(Group=Group, row.names = someRowNames )

dds <- DESeqDataSetFromMatrix(countData = countsTable, colData = coldata1, design = ~ Group  )
dds <- DESeq(dds)

using pre-existing size factors
estimating dispersions
found already estimated dispersions, replacing these
gene-wise dispersion estimates
Error in model.matrix.formula(design(object), data = colData(object)) : 
  data must be a data.frame

If i save the data, restart, reload the data and start over without loading my package it works like a charm.

Whats going on? I'm guessing I might be overriding some function(s), but how do I tell which?

[mikael.lenz.strube]

So I went snooping around in the DESeq2 source code and found the offending call in the function estimateDispersionsGeneEst: modelMatrix <- model.matrix(design(object), data=colData(object))

so ran this

for( i in listOfPackagesMyPackageDependsOn)
{
  library(i, character.only = T)
  print(i)
  model.matrix( object  = design(dds),  data = data.frame( colData(dds) ) )
}

to see where it fails, and it turns out that its "agricolae": DESeq2 is not compatible with agricolae. It does not help to detach agricolae.

Funny thing is that model.matrix needs three arguments when agricolae is loaded and 2 when it is not. Weird.

Any suggestions?

[Devon Ryan]

Wow, redefining model.matrix will certainly cause problems. Maybe reloading the stats package will get around this? agricolae should really be fixed.

[Devon Ryan]

If you load your package first and DESeq2 second then you should get warning messages like "DESeq2 is overriding function foo()". Having said that, what's the output of colData(dds) when you make the DESeqDataSet? You're not doing something crazy like redefining data.frame() I hope ;)

Btw, you might have better luck with this on the bioconductor site. The Bioconductor folks are likely able to track this down a bit faster than us.

[mikael.lenz.strube]

No, no crazy redefinitions, and it doesn't complain at all.

output of colData(dds) is

DataFrame with 47 rows and 3 columns
           Group sizeFactor replaceable
        <factor>  <numeric>   <logical>
41_IL_A        A  1.3980176        TRUE
42_IL_A        A  3.2374793        TRUE
43_IL_A        A  3.7708212        TRUE
44_IL_A        A  0.9493888        TRUE
45_IL_A        A  0.6235690        TRUE
...          ...        ...         ...
84_IL_F        F  0.2930707        TRUE
85_IL_F        F  2.4115678        TRUE
86_IL_F        F  0.3769485        TRUE
87_IL_F        F  0.3268223        TRUE
88_IL_F        F  0.3881433        TRUE

Which I think looks like it should, right?

Maybe ill head over to BionConductor, good idea

[igor]

Post an update if you figure this out. Weird error.

[igor]

When you load a library, it should tell you if there are any conflicts.

[mikael.lenz.strube]

Yeah, I thought so to, but it doesnt say anything. Must be something other than redefined functions then.

[franniehwu]

I have encountered the same problem. It seems that the DEseq2 is really incompatible with the "agricolae." Whenever I run agricolae before the DEseq2, the error comes out. It would be nice if someone can solve this problem, otherwise restart the R will be the only solution I have now.

[igor]

That usually happens if you have the same function name in multiple packages. If you load both packages, you may not be calling the function that you think you are. Instead of calling functions with function(), you can specify the exact package you would like to use with package::function(). That will avoid conflicts.

Answer 1

[mikael.lenz.strube]

I have migrated the question to BioConductor for anyone interested: https://support.bioconductor.org/p/85845/