Question: DESeq2 only runs interaction factor x continuous variable with old not new version
0
gravatar for hmgeiger
3.6 years ago by
hmgeiger10
hmgeiger10 wrote:

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       

 

deseq2 • 1.1k views
ADD COMMENTlink modified 9 months ago by Bioconductor Community ♦♦ 0 • written 3.6 years ago by hmgeiger10
Answer: DESeq2 only runs interaction factor x continuous variable with old not new versi
0
gravatar for Michael Love
3.6 years ago by
Michael Love25k
United States
Michael Love25k wrote:

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)?

ADD COMMENTlink written 3.6 years ago by Michael Love25k

Also see my note in the vignette on using continuous covariates like age. It is in the FAQ.

ADD REPLYlink written 3.6 years ago by Michael Love25k
Answer: DESeq2 only runs interaction factor x continuous variable with old not new versi
0
gravatar for hmgeiger
3.6 years ago by
hmgeiger10
hmgeiger10 wrote:

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"  

ADD COMMENTlink written 3.6 years ago by hmgeiger10

Ah, it should be giving a more helpful error but it looks like my error checking code is broken and I need to fix this in the development branch. Thanks for tipping me to this.

In order to use the style of argument: contrast=c(x,y,z), x needs to be a variable in the design, and y and z need to be levels of that factor. A.Age and B.Age are not levels of Diagnosis, which should throw a helpful error.

From ?results:

contrast: this argument specifies what comparison to extract from the
          ‘object’ to build a results table. one of either:

            • a character vector with exactly three elements: the name
              of a factor in the design formula, the name of the
              numerator level for the fold change, and the name of the
              denominator level for the fold change (simplest case)

But can you tell me what you are trying to extract? The interaction terms are:

results(dds, name="DiagnosisNormal.Age")

results(dds, name="DiagnosisSCZ.Age")

ADD REPLYlink written 3.6 years ago by Michael Love25k
Answer: DESeq2 only runs interaction factor x continuous variable with old not new versi
0
gravatar for hmgeiger
3.6 years ago by
hmgeiger10
hmgeiger10 wrote:

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"  

ADD COMMENTlink written 3.6 years ago by hmgeiger10
Answer: DESeq2 only runs interaction factor x continuous variable with old not new versi
0
gravatar for hmgeiger
3.6 years ago by
hmgeiger10
hmgeiger10 wrote:

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.

ADD COMMENTlink written 3.6 years ago by hmgeiger10

A small note about using the support site: you can use the ADD REPLY or ADD COMMENT buttons to reply to a previous answer. You are adding answers to your own question.

So, you should now be able to extract the interaction terms. If you have questions about their interpretation, check the interaction section of the vignette.

ADD REPLYlink written 3.6 years ago by Michael Love25k

Got it! Thanks.

This is how I ended up setting up my script. A is normal, B and C are my two test conditions.

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

B_vs_control_plus_interaction_age <- results(dds,contrast=list(c("Diagnosis_B_vs_control","Age.B")))

C_vs_control_plus_interaction_age <- results(dds,contrast=list(c("Diagnosis_C_vs_control","Age.C")))

 

ADD REPLYlink modified 3.6 years ago by Michael Love25k • written 3.6 years ago by hmgeiger10

You may want to speak to a local biostatistician to help with the interpretation here.

Note that, if you code Age as a numeric, the "Age" term in the results is the log2 fold change for the A group associated with an increase in Age of 1 unit. And the interaction "Age.B" is the difference between the Age LFC for group B and the Age LFC group A.

ADD REPLYlink written 3.6 years ago by Michael Love25k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 229 users visited in the last hour