Setting Contrast for DESeq2 results This question was also asked on BioStars
I am trying to recreate this heatmap. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.
I need to compare both Ly49+ versus Ly49− and MOG versus MOG plus SP. Column data in at the bottom. The function I used for DESeq2 is,
dds = DESeqDataSetFromMatrix(countData=countData,

colData=colData,

design= ~ cell_treatment)

I am not sure if the design is right. For the results, how should I write the contrasts? I am very confused and would really appreciate any help!
DataFrame with 11 rows and 2 columns

cell treatment

&lt;character&gt; &lt;character&gt;

MOGSP_2 SPtetramer+CD8+TCell MOGSP

MOGSP_3 SPtetramer+CD8+TCell MOGSP

MOG_3 SPtetramer+CD8+TCell MOG

MOG_4 SPtetramer+CD8+TCell MOG

MOG_5 SPtetramer+CD8+TCell MOG

Ly49N_1 Ly49-CD8+TCell MOGSP

Ly49P_1 Ly49+CD8+TCell MOGSP

Ly49N_2 Ly49-CD8+TCell MOGSP

Ly49P_2 Ly49+CD8+TCell MOGSP

Ly49N_3 Ly49-CD8+TCell MOGSP

Ly49P_3 Ly49+CD8+TCell MOGSP

Thank you!!
Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

I generate an example dataset:
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

You can compare them like this:
head(results(dds,c(&quot;Ly&quot;,&quot;Ly49plus&quot;,&quot;Ly49minus&quot;)))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus 
Wald test p-value: Ly Ly49plus vs Ly49minus 
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         &lt;numeric&gt;          &lt;numeric&gt;        &lt;numeric&gt;          &lt;numeric&gt;
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          &lt;numeric&gt;         &lt;numeric&gt;
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

To compare groups within another factor:
results(dds,c(&quot;treatment&quot;,&quot;MOGSP&quot;,&quot;MOG&quot;))
log2 fold change (MLE): treatment MOGSP vs MOG 
Wald test p-value: treatment MOGSP vs MOG 
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         &lt;numeric&gt;          &lt;numeric&gt;        &lt;numeric&gt;          &lt;numeric&gt;
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           &lt;numeric&gt;         &lt;numeric&gt;
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

The colData i used:
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L, 
1L, 2L, 1L, 2L), .Label = c(&quot;Ly49-CD8+TCell&quot;, &quot;Ly49+CD8+TCell&quot;, 
&quot;SPtetramer+CD8+TCell&quot;), class = &quot;factor&quot;), Ly = c(&quot;SPtetramerplus&quot;, 
&quot;SPtetramerplus&quot;, &quot;SPtetramerplus&quot;, &quot;SPtetramerplus&quot;, &quot;SPtetramerplus&quot;, 
&quot;Ly49minus&quot;, &quot;Ly49plus&quot;, &quot;Ly49minus&quot;, &quot;Ly49plus&quot;, &quot;Ly49minus&quot;, 
&quot;Ly49plus&quot;), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L), .Label = c(&quot;MOG&quot;, &quot;MOGSP&quot;), class = &quot;factor&quot;)), class = &quot;data.frame&quot;, row.names = c(NA, 
-11L))