Time-course for three-factor design
Entering edit mode
boczniak767 ▴ 720
Last seen 5 hours ago


*** Background ***

I have following design (currently test data):

4 inbred lines (A, B, C, D)

2 temperatures (x, c) "c" is a control

4 time points (possibly more in future)

Each combination of the above is repeated 3 times

  • possibly "c" and "x" could be joined somehow because plants grow first in "x" temperature and after say 24 h it is switched to "c" control temperature for regrowth. And each time three plants from each line (plants grows in blocks) are pooled for sample. So for say line "A" the full profile could look like:

Ax1 Ax2 ... AxN Ac1 Ac2 ... AcN

** Comparisons of interest **

I want to compare time course in "x" and "c" temperature for each line. I'd also like to see if lines reaction to "x" temperature is different.

*** Data *

Based on a real count data I've made matrix with 96 libraries (4 lines 2 temperatrures 4 time points * 3 biological replications).

I've also made colData object, here are first 48 lines of it

      linia temp czas
Ax15.1     A       x   15
Ax15.2     A       x   15
Ax15.3     A       x   15
Ax22.1     A       x   22
Ax22.2     A       x   22
Ax22.3     A       x   22
Ax29.1     A       x   29
Ax29.2     A       x   29
Ax29.3     A       x   29
Ax36.1     A       x   36
Ax36.2     A       x   36
Ax36.3     A       x   36
Ac15.1     A       c   15
Ac15.2     A       c   15
Ac15.3     A       c   15
Ac22.1     A       c   22
Ac22.2     A       c   22
Ac22.3     A       c   22
Ac29.1     A       c   29
Ac29.2     A       c   29
Ac29.3     A       c   29
Ac36.1     A       c   36
Ac36.2     A       c   36
Ac36.3     A       c   36
Bx15.1     B       x   15
Bx15.2     B       x   15
Bx15.3     B       x   15
Bx22.1     B       x   22
Bx22.2     B       x   22
Bx22.3     B       x   22
Bx29.1     B       x   29
Bx29.2     B       x   29
Bx29.3     B       x   29
Bx36.1     B       x   36
Bx36.2     B       x   36
Bx36.3     B       x   36
Bc15.1     B       c   15
Bc15.2     B       c   15
Bc15.3     B       c   15
Bc22.1     B       c   22
Bc22.2     B       c   22
Bc22.3     B       c   22
Bc29.1     B       c   29
Bc29.2     B       c   29
Bc29.3     B       c   29
Bc36.1     B       c   36
Bc36.2     B       c   36
Bc36.3     B       c   36

** Code and error **

So I proceeded to data set creation

dds=DESeqDataSetFromMatrix(countData=cts, colData=sym.coldata, design= ~ linia + temp + czas + linia:temp:czas)

but it gives an error

converting counts to integer mode
Error in checkFullRank(modelMatrix) : 
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':


I've read the vignette but I can see that my data are balanced. So I don't understand the problem. I've used similar designs in standard anovas and they was ok.

he same problem is for model (after all I'm not interested in the main "line" and "temp" effects)

dds=DESeqDataSetFromMatrix(countData=cts, colData=sym.coldata, design= ~ linia:temp + linia:temp:czas)

It gave the same error as above

To simplify I made subset of data to only have two lines (not four) but the problem remain.

Surely it is problem with my statistical ignorance. My plan "maximum" is to try also edgeR and maSigPro (at least the latter).

** Could someone help me with this design, please? *

I'd like to note, that I've read available documentation about time-course analysis in DESeq2 and successfully repeated it for design with two lines (mutant and wt) and four time points. So I'm able to analyse such design.

> sessionInfo()
R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

[1] C

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

other attached packages:
 [1] DESeq2_1.38.3               SummarizedExperiment_1.28.0
 [3] Biobase_2.58.0              MatrixGenerics_1.10.0      
 [5] matrixStats_0.63.0          GenomicRanges_1.50.2       
 [7] GenomeInfoDb_1.34.9         IRanges_2.32.0             
 [9] S4Vectors_0.36.2            BiocGenerics_0.44.0        

loaded via a namespace (and not attached):
 [1] KEGGREST_1.38.0        locfit_1.5-9.7         tidyselect_1.2.0      
 [4] lattice_0.20-45        generics_0.1.3         colorspace_2.1-0      
 [7] vctrs_0.5.2            utf8_1.2.3             blob_1.2.3            
[10] XML_3.99-0.13          rlang_1.0.6            pillar_1.8.1          
[13] glue_1.6.2             DBI_1.1.3              BiocParallel_1.32.5   
[16] bit64_4.0.5            RColorBrewer_1.1-3     GenomeInfoDbData_1.2.9
[19] lifecycle_1.0.3        zlibbioc_1.44.0        Biostrings_2.66.0     
[22] munsell_0.5.0          gtable_0.3.1           codetools_0.2-19      
[25] memoise_2.0.1          geneplotter_1.76.0     fastmap_1.1.1         
[28] parallel_4.2.2         AnnotationDbi_1.60.0   fansi_1.0.4           
[31] Rcpp_1.0.10            xtable_1.8-4           scales_1.2.1          
[34] cachem_1.0.7           DelayedArray_0.24.0    annotate_1.76.0       
[37] XVector_0.38.0         bit_4.0.5              ggplot2_3.4.1         
[40] png_0.1-8              dplyr_1.1.0            grid_4.2.2            
[43] cli_3.6.0              tools_4.2.2            bitops_1.0-7          
[46] magrittr_2.0.3         RCurl_1.98-1.10        RSQLite_2.3.0         
[49] tibble_3.1.8           pkgconfig_2.0.3        crayon_1.5.2          
[52] Matrix_1.5-3           httr_1.4.5             R6_2.5.1              
[55] compiler_4.2.2
DESeq2 • 85 views
Entering edit mode
Last seen 1 day ago
United States

Sorry I don't have sufficient time for answering non-software related questions on the support site. I recommend consulting with a local statistician.


Login before adding your answer.

Traffic: 321 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6