Dear all,

I have a dataset with several time points and after working on the whole dataset, I was interested on analyzing a specific question related to the time effect. I used the dds object and applied a subsetting like explained in the "Beginner's guide to using the DESeq2 package". But then when I droped some levels that didn't have anymore samples (these levels originally appeared in the design formula) and run DESeq 2 with the the subsets left, it doesn't work and I get the following message :

>  dds0 <- DESeq(dds0)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
  contrasts can be applied only to factors with 2 or more levels

Here is my code for the beginning (only for one time point):

##Assign conditions (first four are c1, second four are c2)
conditions <- factor(c(rep("c1", 12), rep("c2", 12)))

time <- factor(c(rep("0", 3), rep("3", 3),rep("5", 3),rep("7", 3)
                 , rep("0", 3), rep("3", 3),rep("5", 3),rep("7", 3)))

rep <- factor(c(rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1),
                rep("A", 1), rep("B", 1),rep("C", 1)))

batch <- factor(c(rep("b1",17), rep("b2",1), rep("b1",6)))

##Fake metadata
coldata <- data.frame(row.names=colnames(cts), conditions, time, rep, batch)

#---Run DESeq2 for batch, rep, conditions and time---
##Preparation of counts table for all libraries
dds <- DESeqDataSetFromMatrix(countData = cts, 
                                colData = coldata, 
                                design = ~ batch + rep + time + conditions)
dds$conditions <- relevel(dds$conditions, "c1") #make sure control is first, here c1
dds

dds <- DESeq(dds)

##--Per time point c1 vs c2
##-->First interested on the DEG at each time point between c1 vs c2 with c1 = control/base_level

###For 0
#Substet and data preparation
dds0 <- dds[ , dds$time == "0" ] #Subset for analysis of interest
dds0$time <- droplevels(dds0$time) #drop other levels of time
dds0$batch <- droplevels(dds0$batch) #drop other levels of batch
dds0$conditions <- relevel(dds0$conditions, "c1") #put control first here c1
dds0
as.data.frame(colData(dds0)) #To check if we applied the right subset

###Run DESeq on subset data for 0
dds0 <- DESeq(dds0)

It was working a month ago but now drop levels don't work anymore and I don't understand why.

Thank you in advance for your help !

Here the session info :

sessionInfo() R version 4.0.2 (2020-06-22) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale: [1] LCCOLLATE=FrenchFrance.1252 LCCTYPE=FrenchFrance.1252 LCMONETARY=FrenchFrance.1252 [4] LCNUMERIC=C LCTIME=French_France.1252

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

other attached packages: [1] ggrepel0.8.2 ggplot23.3.2 pheatmap1.0.12
[4] ashr2.2-47 DESeq21.28.1 SummarizedExperiment1.18.2 [7] DelayedArray0.14.0 matrixStats0.56.0 Biobase2.48.0
[10] GenomicRanges1.40.0 GenomeInfoDb1.24.2 IRanges2.22.2
[13] S4Vectors0.26.1 BiocGenerics0.34.0 RColorBrewer1.1-2
[16] gplots3.0.4

loaded via a namespace (and not attached): [1] bitops1.0-6 fs1.4.1 usethis1.6.1 devtools2.3.1
[5] bit640.9-7 rprojroot1.3-2 tools4.0.2 backports1.1.7
[9] R62.4.1 irlba2.3.3 KernSmooth2.23-17 DBI1.1.0
[13] colorspace1.4-1 withr2.2.0 tidyselect1.1.0 prettyunits1.1.1
[17] processx3.4.2 bit1.1-15.2 compiler4.0.2 cli2.0.2
[21] desc1.2.0 caTools1.18.0 scales1.1.1 SQUAREM2020.3
[25] genefilter1.70.0 callr3.4.3 mixsqp0.3-43 digest0.6.25
[29] XVector0.28.0 pkgconfig2.0.3 sessioninfo1.1.1 invgamma1.1
[33] rlang0.4.6 rstudioapi0.11 RSQLite2.2.0 generics0.0.2
[37] BiocParallel1.22.0 gtools3.8.2 dplyr1.0.0 RCurl1.98-1.2
[41] magrittr1.5 GenomeInfoDbData1.2.3 Matrix1.2-18 Rcpp1.0.4.6
[45] munsell0.5.0 fansi0.4.1 lifecycle0.2.0 zlibbioc1.34.0
[49] pkgbuild1.1.0 grid4.0.2 blob1.2.1 gdata2.18.0
[53] crayon1.3.4 lattice0.20-41 splines4.0.2 annotate1.66.0
[57] locfit1.5-9.4 ps1.3.3 pillar1.4.6 geneplotter1.66.0
[61] pkgload1.1.0 XML3.99-0.3 glue1.4.1 remotes2.2.0
[65] BiocManager1.30.10 vctrs0.3.1 testthat2.3.2 gtable0.3.0
[69] purrr0.3.4 assertthat0.2.1 xtable1.8-4 survival3.1-12
[73] truncnorm1.0-8 tibble3.0.1 AnnotationDbi1.50.3 memoise1.1.0
[77] ellipsis_0.3.1

Can you try to build a DESeqDataSet with the subset data, it may give you some indication as to why DESeq() cannot run.

The constructor function runs a number of checks to make sure you don't end up with an in-estimable object. This may point you to the issue.