I have successfully run DEXSeq using a simple design; two experimental conditions, each with two replicates. However, I am now trying to run a more complex model, and I am getting the following error:

> fullModel <- ~sample + exon + treatment:exon + time:exon + treatment:time:exon
> reducedModel <- ~sample + exon + treatment:exon + time:exon
> dxd <- testForDEU(dxd, BPPARAM = BPPARAM, fullModel = fullModel, reducedModel = reducedModel )
using supplied model matrix
using supplied model matrix
using supplied model matrix
using supplied model matrix
using supplied model matrix
using supplied model matrix
using supplied model matrix
using supplied model matrix
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘mcols’ for signature ‘"matrix"’

 I have nine samples in total; three treatment times (0h, 3h and 6h), each with three biological replicates. Ideally I would like to identify exons which are differential used in response to the length of treatment. I have arrived at the following design;

> fullModel <- ~sample + exon + treatment:exon + time:exon + treatment:time:exon
> reducedModel <- ~sample + exon + treatment:exon + time:exon

Prior to running the testForDEU command above (where the error occurs), I have used the following commands:

dxd <- DEXSeqDataSetFromHTSeq(  
  countFiles, 
  sampleData=sampleTable, 
  design= ~ sample + exon + treatment:exon + time:exon + treatment:time:exon, 
  flattenedfile=flattenedFile)

dxd <- estimateSizeFactors(dxd)
dxd <- estimateDispersions(dxd, formula = fullModel)

From reading other posts I thought the issue may be out of date dependencies, but I think this is not the case. Any help greatly appreciated.

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8   
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] RColorBrewer_1.1-2         biomaRt_2.26.1             DEXSeq_1.16.10             DESeq2_1.10.1              RcppArmadillo_0.6.600.4.0 
 [6] Rcpp_0.12.3                SummarizedExperiment_1.0.2 GenomicRanges_1.22.4       GenomeInfoDb_1.6.3         IRanges_2.4.8             
[11] S4Vectors_0.8.11           Biobase_2.30.0             BiocGenerics_0.16.1        BiocParallel_1.4.3        

loaded via a namespace (and not attached):
 [1] futile.logger_1.4.1  plyr_1.8.3           XVector_0.10.0       bitops_1.0-6         futile.options_1.0.0 tools_3.2.3         
 [7] zlibbioc_1.16.0      statmod_1.4.24       rpart_4.1-10         RSQLite_1.0.0        annotate_1.48.0      gtable_0.2.0        
[13] lattice_0.20-33      DBI_0.3.1            gridExtra_2.2.1      stringr_1.0.0        hwriter_1.3.2        genefilter_1.52.1   
[19] cluster_2.0.3        Biostrings_2.38.4    locfit_1.5-9.1       grid_3.2.3           nnet_7.3-12          AnnotationDbi_1.32.3
[25] XML_3.98-1.4         survival_2.38-3      foreign_0.8-66       latticeExtra_0.6-28  Formula_1.2-1        magrittr_1.5        
[31] geneplotter_1.48.0   ggplot2_2.1.0        lambda.r_1.1.7       Rsamtools_1.22.0     Hmisc_3.17-2         scales_0.4.0        
[37] splines_3.2.3        xtable_1.8-2         colorspace_1.2-6     stringi_1.0-1        acepack_1.3-3.3      RCurl_1.95-4.8      
[43] munsell_0.4.3   

Hi MedwayC, 

Thanks for reporting this! It was not a bug at the end, but the error message was very uninformative (it has now a nicer one in the latest development version). The design that you are specifying would work if you have a fully crossed design, i.e. where you would have samples for all possible combinations of all the levels across your two response variables (time and treatment). In your colData there are have two variables, but the design is not fully crossed (all the time 0's are only untreated, there are no treated samples with time 0, etc). When removing the dependent columns of your models, your full model and reduced model turn out to be identical and the test was failing.

If you want to test exons for whether exon usage increase or decrease over time you could convert the time column to a numerical value and test using the design below:

colData(dxd)$time <- as.numeric(as.character(colData(dxd)$time))

design(dxd) <- ~sample + exon + time:exon

 If you want to test exon usage effect differences in at least one of the time points, you would have to use the same 'design(dxd)' as above, but leave the 'time' column as a factor.

 Alejandro