Hi,
I am trying to extract results from a DESeqDataSet object based on a customized contrast list. Depending on whether the design contains an interaction term, the effect/result names change, which makes it quite difficult to extract results in an automatic fashion (e.g. as part of a script or app, without checking resultsNames() manually). Here's an example based on ?results:
library(DESeq2) dds <- makeExampleDESeqDataSet(n=100,m=24) dds$genotype <- factor(rep(rep(c("I","II","III"),each=4),2)) design(dds) <- ~ genotype + condition dds <- DESeq(dds) resultsNames(dds) # [1] "Intercept" "genotypeI" "genotypeII" "genotypeIII" "conditionA" "conditionB"
# genotypeI compared to average of genotypeII and genotypeIII results(dds, contrast=list("genotypeI", c("genotypeII", "genotypeIII")), listValues=c(1, -0.5)) # works!
However, when adding an interaction term, pairwise effects (_vs_) instead of individual effects are extracted, which makes it difficult to extract results for individual levels or customized contrasts:
design(dds) <- ~ genotype + condition + genotype:condition dds <- DESeq(dds) resultsNames(dds) # [1] "Intercept" "genotype_II_vs_I" "genotype_III_vs_I" "condition_B_vs_A" "genotypeII.conditionB" "genotypeIII.conditionB"
results(dds, contrast=list("genotypeI", c("genotypeII", "genotypeIII")), listValues=c(1, -0.5)) # error
Additionally, there seems to be no way to extract results for e.g. genotypeI.conditionA vs ...
DESeq(dds, modelMatrixType="expanded") requires a beta prior, but betaPrior=TRUE results in: Error in designAndArgChecker(object, betaPrior) : betaPrior=FALSE should be used for designs with interactions
Is there a way to directly influence which effects are calculated (e.g. genotypeI and genotypeII instead of genotype_II_vs_I, or genotypeII.conditionA instead of genotypeII.conditionB)? (For cases where resultsNames(dds) are not checked manually)
Thanks.
I don't at this time have a convenience function to show the results names but you can run it on a few rows e.g. 10, if you want to see in advance. For coefficients that are built as comparisons to reference level (as in designs without the beta prior) we add that reference level in the result name for clarity to end users, for example
condition_B_vs_A
.Actually, the rules are pretty simple:
For designs with betaPrior, you get the same resultsNames as you would get with ~0 + condition, e.g.
conditionB
, except that in the designs with betaPrior we also have an intercept. Likewise, obviously for the designs which turn off beta prior but have ~0 + condition you get the standard terms you would get with model.matrix, e.gconditionB
.For the designs without beta prior, but which have an intercept term, you get the
condition_B_vs_A
style of resultsNames, to emphasize that the coefficient is built as a comparison against the reference level. Interaction terms force the beta prior off, so you get these style of terms for the main effects. The interaction terms use the default style from model.matrix, e.g.conditionB.genotypeII
.Pre-calculation using a few rows is a good idea. Thanks for the fast reply!