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

<pre><code>dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)
</code></pre>

<p>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").</p>

<p>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 <code>~0+group</code> 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 <code>treatment HS30</code> and <code>geno prg1</code>.</p>

CC BY-SA 4.0
2018-06-16T22:08:10.823Z
"The simplest manner is to not use a wald test, but rather an LRT" -- could you explain why?
Unknown artist
CC BY-SA 4.0
2018-06-17T14:02:27.237Z
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.
Devon Ryan
CC BY-SA 3.0
2017-05-19T16:07:47.350Z
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
CC BY-SA 3.0
2017-05-19T18:05:01.773Z
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.
Devon Ryan
CC BY-SA 3.0
2017-05-19T14:51:01.017Z
2017-05-19T14:51:01.017Z
Devon Ryan
<p>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.</p>

<p>It seems possible do do it directly when building the <code>colData</code> and when calling <code>DESeqDataSetFromMatrix</code>:</p>

<p>Let's add a combined "geno" and "treat" factors to the future <code>colData</code> parameter:</p>

<pre><code>&gt; col_data$geno_treat &lt;- as.factor(paste(col_data$geno, col_data$treat, sep="_"))
&gt; col_data
DataFrame with 12 rows and 4 columns
                        geno       treat         rep     geno_treat
                 &lt;character&gt; &lt;character&gt; &lt;character&gt;       &lt;factor&gt;
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
</code></pre>

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

<pre><code>&gt; dds &lt;- DESeqDataSetFromMatrix(
    countData = counts_data,
    colData = col_data,
    design = ~ geno_treat)
</code></pre>

<p>We run the analysis:</p>

<pre><code>&gt; dds &lt;- DESeq(dds)
</code></pre>

<p>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":</p>

<pre><code>res &lt;- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))
</code></pre>

CC BY-SA 3.0
2017-05-22T12:05:18.343Z
2017-05-22T12:05:18.343Z
bli
<p>I have a set of high-troughput experiments with 2 genotypes (&quot;WT&quot; and &quot;prg1&quot;) and 3 treatments (&quot;RT&quot;, &quot;HS30&quot; and &quot;HS30RT120&quot;), and there are 2 replicates for each of the genotype x treatment combinations.</p>
<p>The read counts for the genes are summarized in a file that I load as follows in R:</p>
<pre><code>&gt; counts_data &lt;- read.table(&quot;path/to/my/file&quot;, header=TRUE, row.names=&quot;gene&quot;)
&gt; colnames(counts_data)
 [1] &quot;WT_RT_1&quot;          &quot;WT_HS30_1&quot;        &quot;WT_HS30RT120_1&quot;   &quot;prg1_RT_1&quot;       
 [5] &quot;prg1_HS30_1&quot;      &quot;prg1_HS30RT120_1&quot; &quot;WT_RT_2&quot;          &quot;WT_HS30_2&quot;       
 [9] &quot;WT_HS30RT120_2&quot;   &quot;prg1_RT_2&quot;        &quot;prg1_HS30_2&quot;      &quot;prg1_HS30RT120_2&quot;
</code></pre>
<p>I describe the experiments as follows:</p>
<pre><code>&gt; col_data &lt;- DataFrame(
    geno = c(rep(&quot;WT&quot;, times=3), rep(&quot;prg1&quot;, times=3), rep(&quot;WT&quot;, times=3), rep(&quot;prg1&quot;, times=3)),
    treat = rep(c(&quot;RT&quot;, &quot;HS30&quot;, &quot;HS30RT120&quot;), times=4),
    rep = c(rep(&quot;1&quot;, times=6), rep(&quot;2&quot;, times=6)),
    row.names = colnames(counts_data))
&gt; col_data
DataFrame with 12 rows and 3 columns
                        geno       treat         rep
                 &lt;character&gt; &lt;character&gt; &lt;character&gt;
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
</code></pre>
<p>I want to build a DESeq2 object that I could use to either:</p>
<ul>
<li><strong>find differentially expressed genes when the treatment varies for a given fixed genotype</strong></li>
</ul>
<p>or:</p>
<ul>
<li><strong>find differentially expressed genes when the genotype varies for a given fixed treatment</strong></li>
</ul>
<p>In the <a href="https://support.bioconductor.org/p/74101/#74163" rel="noreferrer">bioconductor help forum</a> I think I've found a somewhat similar situation, and I read the following:</p>
<blockquote>
<p>Try a design of ~ genotype + genotype:condition</p>
<p>Then you will have a condition effect for each level of genotype, including the reference level.</p>
<p>You can constrast pairs of them using the list style of the 'contrast' argument.</p>
</blockquote>
<p>However, this doesn't explain how to apply this &quot;list style&quot; to the &quot;contrast&quot; 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.</p>
<p>So I tried the following more symmetric formula:</p>
<pre><code>&gt; dds &lt;- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
&gt; dds &lt;- DESeq(dds)
</code></pre>
<p>Now, can I for instance <strong>get the differential expression results when comparing treatment &quot;HS30&quot; against &quot;RT&quot; as a reference, in genotype &quot;prg1&quot;</strong>?</p>
<p>And how?</p>
<p>If I understand correctly, the above-mentioned &quot;list style&quot; uses names given by the <code>resultsNames</code> function. In my case, I have the following:</p>
<pre><code>&gt; resultsNames(dds)
[1] &quot;Intercept&quot;               &quot;geno_WT_vs_prg1&quot;        
[3] &quot;treat_HS30RT120_vs_HS30&quot; &quot;treat_RT_vs_HS30&quot;       
[5] &quot;genoWT.treatHS30RT120&quot;   &quot;genoWT.treatRT&quot;
</code></pre>
<p>I guess I would need a contrast between &quot;genoprg1.treatRT&quot; and a &quot;genoprg1.treatHS30&quot;, but these are not in the above results names.</p>
<p>I'm lost.</p>

<r><bioconductor><deseq2><differential-expression>
Understanding DESeq2 design, contrast and results
# Understanding DESeq2 design, contrast and results

(Asked by: username_0 on 2017-05-19)

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
                   
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.

**Tags:** r, bioconductor, deseq2, differential-expression

---

## 2 answers

### Answer by username_0 on 2017-05-22. Score: 10 | Accepted answer

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


---

### Answer by username_1 on 2017-05-19. Score: 16

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.

> **Comments:**
> **username_0**: 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.  
> **username_1**: 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.  
> **username_2**: "The simplest manner is to not use a wald test, but rather an LRT" -- could you explain why?  
> **username_1**: 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.  

---

