Hello,

I am trying to run DEXSeq analysis for RNA-Seq data from mice. I want to see if there is differential exon usage between the experimental group and control group. I am including blocking for individuals since there are duplicate mice from different sequencing runs in my dataset. I also have Library Kit in the formulae since a different kit was used for the two different runs.

Here is the information for each sample:

"Sample.ID" "Library.kit" "Group" "IndSample"
"GV1"  "Nextflex" "CTRL" "GV1"
"GV2" "Nextflex" "EXP" "GV2"
"GV3" "Nextflex" "EXP" "GV3"
"GV4" "Nextflex" "EXP" "GV4"
"GV5" "Nextflex" "CTRL" "GV5"
"GV6" "Nextflex" "CTRL" "GV6"
"KG5.1" "Illumina_truseq" "EXP" "GV2"
"KG5.2" "Illumina_truseq" "EXP" "GV3"
"KG5.3" "Illumina_truseq" "EXP" "GV4"
"KG5.4" "Illumina_truseq" "EXP" "KG5.4"
"KG5.5" "Illumina_truseq" "CTRL" "GV5"
"KG5.6" "Illumina_truseq" "CTRL" "GV6"
"KG5.7" "Illumina_truseq" "CTRL" "GV1"

So for example, Sample.ID "GV1" comes from the same individual as Sample.ID "KG5.7" as indicated by the IndSample column.

My full model is: ~Library.kit + Group + exon + IndSample:exon + Group:exon My reduced model (redMod) is: ~Library.kit + Group + exon + IndSample:exon

Creating the DEXSeq object, estimating size factors, and estimating dispersions all runs normally.

But when I run:

dxd <- testForDEU(dxd, fullModel = design(dxd), reducedModel = redMod, BiocParallel::MulticoreParam(workers = 15))

I receive this error

Error: BiocParallel errors
  element index: 1, 2, 3, 4, 5, 6, ...
  first error: less than one degree of freedom, perhaps full and reduced models are not in the correct order

Running the code without BiocParallel still outputs an error

 dxd <- testForDEU(dxd, fullModel = design(dxd), reducedModel = redMod)
Error in nbinomLRT(x, reduced = reducedModelMatrix, full = fullModelMatrix) : 
  less than one degree of freedom, perhaps full and reduced models are not in the correct order

When I remove the IndSample:exon term from the formulae, the testForDEU() function runs fine, so I think the error lies in that term somehow, but I'm not sure how to fix it.

Thank you for any help!

Session Info:

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8    
 [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_GB.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] DEXSeq_1.30.0               DESeq2_1.24.0               argparser_0.4               data.table_1.12.2          
 [5] dplyr_0.8.3                 RColorBrewer_1.1-2          AnnotationDbi_1.46.0        SummarizedExperiment_1.14.1
 [9] DelayedArray_0.10.0         matrixStats_0.54.0          GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[13] IRanges_2.18.1              S4Vectors_0.22.0            Biobase_2.44.0              BiocGenerics_0.30.0        
[17] BiocParallel_1.18.1        

loaded via a namespace (and not attached):
 [1] httr_1.4.1             bit64_0.9-7            splines_3.6.1          Formula_1.2-3          assertthat_0.2.1      
 [6] statmod_1.4.32         latticeExtra_0.6-28    blob_1.2.0             Rsamtools_2.0.0        GenomeInfoDbData_1.2.1
[11] yaml_2.2.0             progress_1.2.2         pillar_1.4.2           RSQLite_2.1.2          backports_1.1.4       
[16] lattice_0.20-38        glue_1.3.1             digest_0.6.20          XVector_0.24.0         checkmate_1.9.4       
[21] colorspace_1.4-1       htmltools_0.3.6        Matrix_1.2-17          XML_3.98-1.20          pkgconfig_2.0.2       
[26] biomaRt_2.40.3         genefilter_1.66.0      zlibbioc_1.30.0        purrr_0.3.2            xtable_1.8-4          
[31] scales_1.0.0           htmlTable_1.13.1       tibble_2.1.3           annotate_1.62.0        ggplot2_3.2.1         
[36] nnet_7.3-12            lazyeval_0.2.2         survival_2.44-1.1      magrittr_1.5           crayon_1.3.4          
[41] memoise_1.1.0          hwriter_1.3.2          foreign_0.8-72         prettyunits_1.0.2      tools_3.6.1           
[46] hms_0.5.0              stringr_1.4.0          locfit_1.5-9.1         munsell_0.5.0          cluster_2.1.0         
[51] Biostrings_2.52.0      compiler_3.6.1         rlang_0.4.0            grid_3.6.1             RCurl_1.95-4.12       
[56] rstudioapi_0.10        htmlwidgets_1.3        bitops_1.0-6           base64enc_0.1-3        gtable_0.3.0          
[61] DBI_1.0.0              R6_2.4.0               gridExtra_2.3          knitr_1.24             bit_1.1-14            
[66] zeallot_0.1.0          Hmisc_4.2-0            stringi_1.4.3          Rcpp_1.0.2             vctrs_0.2.0           
[71] geneplotter_1.62.0     rpart_4.1-15           acepack_1.4.1          tidyselect_0.2.5       xfun_0.8

Hi @acciardo,

I think your formulae are not well set up. If you want to add Library.kit and IndSample as blocking factors, your formulae should look like this:

formulaFullModel    =  ~ sample + exon + IndSample:exon + Library.kit:exon + Group:exon
formulaReducedModel =  ~ sample + exon + IndSample:exon + Library.kit:exon

As an additional comment, the sample KG5.4 does not have a corresponding pair in the column IndSample. I think this can result in a matrix that is not of full rank and you might get an error. If you want to add IndSample as blocking factor, you might have to remove this sample.