Answers.Body: <p>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:</p>
  <pre><code>                   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
  </code></pre>
  <p>I generate an example dataset:</p>
  <pre><code>counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)
  </code></pre>
  <p>And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:</p>
  <pre><code>dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
  dds = DESeq(dds)
  </code></pre>
  <p>You can compare them like this:</p>
  <pre><code>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
  </code></pre>
  <p>To compare groups within another factor:</p>
  <pre><code>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
  </code></pre>
  <p>The colData i used:</p>
  <pre><code>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))
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-16T05:27:46.427Z
Answers.Text: so it can used for multiple variable or cell type comparisons instead of many pairwise comparison which i do normally
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-16T16:15:23.957Z
Answers.Text: it depends on the design but in general yes. In this example, you can see it's necessary because the groups have some overlap
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-21T06:36:38.923Z
Answers.Text: conceptual doubt so if i have to perform multiple pairwise comparison do i need to set reference level for every comparison or i can simply use "results(dds,c("treatment","MOGSP","MOG"))" or the contrast function?
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-22T01:09:10.543Z
Answers.Text: using ```results(dds,c("treatment","MOGSP","MOG"))``` takes care of everything. The reference level of the factor does not matter here
Answers.UserDisplayName: StupidWolf
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-22T11:12:18.613Z
Answers.Text: will try this..
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-13T18:34:53.213Z
Answers.LastActivityDate: 2020-09-13T18:34:53.213Z
Answers.UserDisplayName: StupidWolf
Body: <p><em>This question was also asked on <a href="https://www.biostars.org/p/461091/" rel="nofollow noreferrer">BioStars</a></em></p>
  <p>I am trying to recreate <a href="https://www.nature.com/articles/s41586-019-1467-x/figures/6" rel="nofollow noreferrer">this heatmap</a>. How would I compare the four variables: Ly49+/- and MOG/MOGSP to get a single heatmap? Thank you for helping me through this.</p>
  <p>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,</p>
  <pre class="lang-r prettyprint-override"><code>dds = DESeqDataSetFromMatrix(countData=countData,
                                  colData=colData,
                                  design= ~ cell_treatment)
  </code></pre>
  <p>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!</p>
  <pre class="lang-r prettyprint-override"><code>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
  </code></pre>
  <p>Thank you!</p>
  
Tags: <rna-seq><deseq2><design>
Title: Setting Contrast for DESeq2 results
ThreadText: # Setting Contrast for DESeq2 results
  
  (Asked by: username_0 on 2020-09-13)
  
  *This question was also asked on [BioStars](https://www.biostars.org/p/461091/)*
  
  I am trying to recreate [this heatmap](https://www.nature.com/articles/s41586-019-1467-x/figures/6). 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,
  
  ```r
  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!
  
  ```r
  DataFrame with 11 rows and 2 columns
  cell treatment
  <character> <character>
  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!
  
  **Tags:** rna-seq, deseq2, design
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-09-13. Score: 3 | Accepted answer
  
  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("Ly","Ly49plus","Ly49minus")))
  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
           <numeric>          <numeric>        <numeric>          <numeric>
  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
            <numeric>         <numeric>
  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("treatment","MOGSP","MOG"))
  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
           <numeric>          <numeric>        <numeric>          <numeric>
  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
             <numeric>         <numeric>
  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("Ly49-CD8+TCell", "Ly49+CD8+TCell", 
  "SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus", 
  "SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus", 
  "Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus", 
  "Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L, 
  2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA, 
  -11L))
  
  ```
  
  > **Comments:**
  > **username_2**: so it can used for multiple variable or cell type comparisons instead of many pairwise comparison which i do normally  
  > **username_1**: it depends on the design but in general yes. In this example, you can see it's necessary because the groups have some overlap  
  > **username_2**: conceptual doubt so if i have to perform multiple pairwise comparison do i need to set reference level for every comparison or i can simply use "results(dds,c("treatment","MOGSP","MOG"))" or the contrast function?  
  > **username_1**: using ```results(dds,c("treatment","MOGSP","MOG"))``` takes care of everything. The reference level of the factor does not matter here  
  > **username_2**: will try this..  
  
  ---
  
  