Trying to run the following code to compare the interaction between three diagnosis groups (call them A, B, and C) and age (a continuous variable, not grouped into factors). The code works fine with DESeq2 1.6.3, but not with 1.10.1.

I checked and diagnosis is definitely a factor and age is definitely an integer.

design <- data.frame(Sample=paste(1:24,"D",sep=""),Diagnosis=factor(rep(c("A","B","C")),times=c(11,7,6)),Age=c(78,75,54,58,54,75,48,47,42,30,25,43,65,51,56,31,50,46,65,50,80,48,59,61))

dds <- DESeq(DESeqDataSetFromMatrix(countData=featureCounts[,paste(1:24,"D",sep="")],colData=design,design=~Diagnosis:Age))

First comparison is this, and is where I get the error. I would eventually also like to compare A vs. C and B vs. C.

results <- results(dds,contrast=c("Diagnosis","A.Age","B.Age"))

Error when running newer DESeq2:

Error in rowSums(cts.sub == 0) :

  'x' must be an array of at least two dimensions

SessionInfo:

R version 3.2.1 (2015-06-18)

Platform: x86_64-unknown-linux-gnu (64-bit)

Running under: CentOS Linux 7 (Core)


locale:

 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              

 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    

 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   

 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 

 [9] LC_ADDRESS=C               LC_TELEPHONE=C            

[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       


attached base packages:

[1] parallel  stats4    methods   stats     graphics  grDevices utils    

[8] datasets  base     


other attached packages:

 [1] DESeq2_1.10.1              RcppArmadillo_0.6.500.4.0 

 [3] Rcpp_0.12.3                SummarizedExperiment_1.0.2

 [5] Biobase_2.30.0             GenomicRanges_1.22.4      

 [7] GenomeInfoDb_1.6.3         IRanges_2.4.8             

 [9] S4Vectors_0.8.11           BiocGenerics_0.16.1       


loaded via a namespace (and not attached):

 [1] RColorBrewer_1.1-2   futile.logger_1.4.1  plyr_1.8.3          

 [4] XVector_0.10.0       futile.options_1.0.0 zlibbioc_1.16.0     

 [7] rpart_4.1-10         RSQLite_1.0.0        annotate_1.48.0     

[10] gtable_0.2.0         lattice_0.20-33      DBI_0.3.1           

[13] gridExtra_2.2.1      genefilter_1.52.1    cluster_2.0.3       

[16] locfit_1.5-9.1       grid_3.2.1           nnet_7.3-12         

[19] AnnotationDbi_1.32.3 XML_3.98-1.4         survival_2.38-3     

[22] BiocParallel_1.4.3   foreign_0.8-66       latticeExtra_0.6-28 

[25] Formula_1.2-1        geneplotter_1.48.0   ggplot2_2.1.0       

[28] lambda.r_1.1.7       Hmisc_3.17-2         scales_0.4.0        

[31] splines_3.2.1        colorspace_1.2-6     xtable_1.8-2        

[34] acepack_1.3-3.3      munsell_0.4.3       

You should include the main effects as well, I think you mean to use a design of ~Diagnosis*Age which expands to ~Diagnosis + Age + Diagnosis:Age.

Can you try this, and then see what you get with resultsNames(dds)?

I tried this, but still get the same error at the results step. 

resultsNames(dds) now returns the following:

[1] "Intercept"               "Diagnosis_Normal_vs_MDD"

[3] "Diagnosis_SCZ_vs_MDD"    "Age"                    

[5] "DiagnosisNormal.Age"     "DiagnosisSCZ.Age"  

Oh, sorry, I was previously using aliases for some of the variables.

Continuing with the aliases, this should read:

[1] "Intercept"               "Diagnosis_B_vs_A"

[3] "Diagnosis_C_vs_A"    "Age"                    

[5] "DiagnosisB.Age"     "DiagnosisC.Age"  

The investigator wants to look at the interaction between age and diagnosis. They wouldn't be opposed to grouping age into factors as appropriate I don't think, though I also wanted to just see if it could be done keeping age as a continuous variable.