# Understanding DESeq2 design, contrast and results

I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.

The read counts for the genes are summarized in a file that I load as follows in R:

> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
[1] "WT_RT_1"          "WT_HS30_1"        "WT_HS30RT120_1"   "prg1_RT_1"
[5] "prg1_HS30_1"      "prg1_HS30RT120_1" "WT_RT_2"          "WT_HS30_2"
[9] "WT_HS30RT120_2"   "prg1_RT_2"        "prg1_HS30_2"      "prg1_HS30RT120_2"


I describe the experiments as follows:

> col_data <- DataFrame(
geno = c(rep("WT", times=3), rep("prg1", times=3), rep("WT", times=3), rep("prg1", times=3)),
treat = rep(c("RT", "HS30", "HS30RT120"), times=4),
rep = c(rep("1", times=6), rep("2", times=6)),
row.names = colnames(counts_data))
> col_data
DataFrame with 12 rows and 3 columns
geno       treat         rep
<character> <character> <character>
WT_RT_1                   WT          RT           1
WT_HS30_1                 WT        HS30           1
WT_HS30RT120_1            WT   HS30RT120           1
prg1_RT_1               prg1          RT           1
prg1_HS30_1             prg1        HS30           1
...                      ...         ...         ...
WT_HS30_2                 WT        HS30           2
WT_HS30RT120_2            WT   HS30RT120           2
prg1_RT_2               prg1          RT           2
prg1_HS30_2             prg1        HS30           2
prg1_HS30RT120_2        prg1   HS30RT120           2


I want to build a DESeq2 object that I could use to either:

• find differentially expressed genes when the treatment varies for a given fixed genotype

or:

• find differentially expressed genes when the genotype varies for a given fixed treatment

In the bioconductor help forum I think I've found a somewhat similar situation, and I read the following:

Try a design of ~ genotype + genotype:condition

Then you will have a condition effect for each level of genotype, including the reference level.

You can constrast pairs of them using the list style of the 'contrast' argument.

However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.

So I tried the following more symmetric formula:

> dds <- DESeqDataSetFromMatrix(
countData = counts_data,
colData = col_data,
design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)


Now, can I for instance get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"?

And how?

If I understand correctly, the above-mentioned "list style" uses names given by the resultsNames function. In my case, I have the following:

> resultsNames(dds)
[1] "Intercept"               "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120"   "genoWT.treatRT"


I guess I would need a contrast between "genoprg1.treatRT" and a "genoprg1.treatHS30", but these are not in the above results names.

I'm lost.

• this is a more similar example (I think): support.bioconductor.org/p/63201. Conclusion is similar to Devon's answer - use LRT. May 19 '17 at 14:59

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:

dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)


The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").

As an aside, this is probably a case where the edgeR-preferred way of creating groups of genotype-treatment combinations and then using a model of ~0+group might make your life a bit easier. You'll get the same results (more or less) regardless, but it'll probably be easier for you to think in those terms rather than remembering that the base level will be treatment HS30 and geno prg1.

• I tried the following: dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno + treat + geno:treat) then dds = DESeq(dds, test="LRT", reduced=~geno+geno:treat) but this fails with Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior = betaPrior, : less than one degree of freedom, perhaps full and reduced models are not in the correct order. Besides, what I want is not "regardless of the type of treatment". I want results for a given fixed treatment or for a given fixed genotype. Maybe my question was not clear.
– bli
May 19 '17 at 16:07
• I'd have to go over the model matrix to see why it's producing the error. Anyway, I would suggest going the route the edgeR authors prefer and making the treatment-genotype combinations into individual groups. Then the contrasts are more obvious. May 19 '17 at 18:05
• "The simplest manner is to not use a wald test, but rather an LRT" -- could you explain why? Jun 16 '18 at 22:08
• Because treatment has multiple levels, so with a wald test one would need to handle the contrasts appropriately, whereas with an LRT one can simultaneously omit all the levels at once. Jun 17 '18 at 14:02

It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts.

It seems possible do do it directly when building the colData and when calling DESeqDataSetFromMatrix:

Let's add a combined "geno" and "treat" factors to the future colData parameter:

> col_data$geno_treat <- as.factor(paste(col_data$geno, col_data\$treat, sep="_"))
> col_data
DataFrame with 12 rows and 4 columns
geno       treat         rep     geno_treat
<character> <character> <character>       <factor>
WT_RT_1                   WT          RT           1          WT_RT
WT_HS30_1                 WT        HS30           1        WT_HS30
WT_HS30RT120_1            WT   HS30RT120           1   WT_HS30RT120
prg1_RT_1               prg1          RT           1        prg1_RT
prg1_HS30_1             prg1        HS30           1      prg1_HS30
...                      ...         ...         ...            ...
WT_HS30_2                 WT        HS30           2        WT_HS30
WT_HS30RT120_2            WT   HS30RT120           2   WT_HS30RT120
prg1_RT_2               prg1          RT           2        prg1_RT
prg1_HS30_2             prg1        HS30           2      prg1_HS30
prg1_HS30RT120_2        prg1   HS30RT120           2 prg1_HS30RT120


We can now use a design where differential expression will be explained by these combined factors:

> dds <- DESeqDataSetFromMatrix(
countData = counts_data,
colData = col_data,
design = ~ geno_treat)


We run the analysis:

> dds <- DESeq(dds)


Then we can query results for a particular contrast between such factor combinations. For instance, to have the results for the effect of treatment "HS30" against the reference state "RT" in genotype "prg1":

res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))