Understanding DESeq2 design, contrast and results 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.
The read counts for the genes are summarized in a file that I load as follows in R:
&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;

I describe the experiments as follows:
&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

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 &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.
So I tried the following more symmetric formula:
&gt; dds &lt;- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
&gt; dds &lt;- DESeq(dds)

Now, can I for instance get the differential expression results when comparing treatment &quot;HS30&quot; against &quot;RT&quot; as a reference, in genotype &quot;prg1&quot;?
And how?
If I understand correctly, the above-mentioned &quot;list style&quot; uses names given by the resultsNames function. In my case, I have the following:
&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;

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.
I'm lost.
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.