Issue/Error with DEXSeq design formulae in testForDEU function
1
0
Entering edit mode
acciardo • 0
@acciardo-19125
Last seen 3.3 years ago

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
deseq2 DEXSeq BioCParallel • 560 views
ADD COMMENT
0
Entering edit mode
Alejandro Reyes ★ 1.8k
@alejandro-reyes-5124
Last seen 6 weeks ago
Novartis Institutes for BioMedical Reseā€¦

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.

ADD COMMENT
0
Entering edit mode

Hi Alejandro,

Thank you for the help - unfortunately I tried changing the formulae to what you suggested and received the same 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

Then, I also removed the KG5.4 sample and once again received the same error. I am really unsure of what could be causing the issue now that the KG5.4 sample is gone.

ADD REPLY

Login before adding your answer.

Traffic: 514 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6