Answers.Body: <p>It's not surprising to me that the adjusted fold change statistic changes in different ways for different genes when covariates are added or removed. That's sort of the point of using covariates, as long as it makes sense to include them (which is an experimental decision, not a mathematical or statistical decision).</p>
  <p>One thing I'd recommend checking is to make sure the results you're looking at are what you <em>think</em> you are looking at. You say that you're interested in testing <code>FAB</code>, but that is not included in the second (Age-only) design, which means that it won't be possible to extract any <code>FAB</code> comparisons from the results. It looks like your first design might be looking at <code>Sex</code> (based on the resultsNames ordering) or one of the <code>FAB</code> comparisons (based on reverse ordering of variables in the model), whereas the second might be looking at <code>Age</code>. This should be explicitly stated when displaying the results, but you haven't included that information in your output.</p>
  <p>If you want <code>FAB</code> comparisons in the results, that variable needs to be in the specified model, e.g. &quot;design = ~ Age + FAB&quot;. The Design 2 model that you have in your question only includes <code>Age</code>.</p>
  <p>A good way to do this without getting the wrong contrast is to fully specify the variables when fetching results, as explained in the <a href="https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs" rel="nofollow noreferrer">DESeq2 vignette</a>:</p>
  <blockquote>
  <p>In any case, the <code>contrast</code> argument of the function results takes a character vector of length three: the name of the variable, the name of the factor level for the numerator of the log2 ratio, and the name of the factor level for the denominator. The contrast argument can also take other forms, as described in the help page for results and below</p>
  </blockquote>
  <pre><code>res.age &lt;- results(dds, name=&quot;Age&quot;)
  res.FAB.M3.M0 &lt;- results(dds, contrast=c(&quot;FAB&quot;,&quot;M3&quot;,&quot;M0&quot;))
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-23T08:52:43.910Z
Answers.Text: `design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB) ` I used this one this resource https://github.com/hbctraining/DGE_workshop_salmon_online/blob/master/homework/DGE_assignment_2_answer_key.R  "which means that it won't be possible to extract any FAB comparisons from the results" i could extract the result for `FAB`
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-23T10:40:43.477Z
Answers.Text: this biostar post https://www.biostars.org/p/278684/ explains i hope im getting it right
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-24T08:54:35.887Z
Answers.Text: We can account for the different types of sequencing, and get a clearer picture of the differences attributable to the treatment. As condition is the variable of interest, we put it at the end of the formula. I followed this ..
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-24T23:57:04.470Z
Answers.Text: What gringer means is you can't get FAB results from design2, because it doesn't include FAB at all.  If you specify the contrast when getting results, **the order of the elements in the design does not matter at all**.   But since you won't say how you called results, no one really knows what you are doing.
Answers.UserDisplayName: swbarnes2
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-25T00:14:27.197Z
Answers.Text: Also, unless you have a huge number of samples, I am skeptical that you can really model so many different variables at once.  I would keep only those elements which correlate to the first few PCs.
Answers.UserDisplayName: swbarnes2
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-26T16:18:23.693Z
Answers.Text: @swbarnes2 yes I tried the second suggestion so if i want to achieve that 80% variance threshold it would be like 40 components may I need to filter or group the gene in better way...
Answers.UserDisplayName: kcm
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-06-23T04:24:27.627Z
Answers.LastActivityDate: 2022-06-25T03:37:25.517Z
Answers.LastEditDate: 2022-06-25T03:37:25.517Z
Answers.UserDisplayName: gringer
Body: <p>Here I have two designs:</p>
  <p><strong>Design 1</strong></p>
  <pre><code> dds &lt;- DESeqDataSetFromMatrix(countData=rpkm_ordered,
             colData=coldata, design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB) 
  </code></pre>
  <p>To see the results I did this</p>
  <pre><code>resultsNames(dds)
   [1] &quot;Intercept&quot;          &quot;Sex_Male_vs_Female&quot; &quot;Age&quot;                &quot;TMB&quot;                &quot;WBC&quot;                &quot;BM_percentage&quot;      &quot;FAB_M1_vs_M0&quot;      
   [8] &quot;FAB_M2_vs_M0&quot;       &quot;FAB_M3_vs_M0&quot;       &quot;FAB_M4_vs_M0&quot;       &quot;FAB_M5_vs_M0&quot;  
  </code></pre>
  <p>I extracted only the Age from the first design</p>
  <p><strong>Design 2</strong></p>
  <p>I again used the design using <code>Age</code> variable :</p>
  <pre><code>dds &lt;- DESeqDataSetFromMatrix(countData=rpkm_ordered,
            colData=coldata, design= ~ Age)
  </code></pre>
  <p>As an example I take this gene this output from the first design</p>
  <pre><code>  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
    &lt;chr&gt;              &lt;dbl&gt;          &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt;      &lt;dbl&gt;
  1 ENSG00000196188     97.4         -0.742  0.113 -6.58 4.57e-11 0.00000125
  </code></pre>
  <p>This from the second design:</p>
  <pre><code>  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
    &lt;chr&gt;              &lt;dbl&gt;          &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt;      &lt;dbl&gt;
    ENSG00000196188     97.4         -0.757  0.119 -6.34 2.28e-10   3.43e-07
          
  </code></pre>
  <p>Another example this gene from the first design; the result is not significant:</p>
  <pre><code>           gene  baseMean log2FoldChange  lfcSE    stat pvalue   padj
  ENSG00000173826  74.54467     0.01283834  0.589  0.0218  0.983  0.993
  </code></pre>
  <p>The same gene from <strong>Design 2</strong> <em>is</em> significant:</p>
  <pre><code>           gene  baseMean log2FoldChange  lfcSE    stat   pvalue     padj
  ENSG00000173826 74.54467       -2.004743  0.190   -10.6 4.76e-26 1.35e-21
  </code></pre>
  <p>Is this drastic outcome is due to the maths behind it or am I doing something wrong?</p>
  <p>So would my interpretation be as such when i consider Age as my only variable I can consider <code>ENSG00000173826</code>significant ? But it is not significant when my main effect which i want to test is FAB which is my design 1.</p>
  <p>I would like to know how to interpret and report the same.</p>
  <p>Any suggestion or help would be really really appreciated</p>
  
Tags: <deseq2>
Title: deseq2 model design : Different gene output
ThreadText: # deseq2 model design : Different gene output
  
  (Asked by: username_0 on 2022-06-21)
  
  Here I have two designs:
  
  **Design 1**
  
  ```
   dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
             colData=coldata, design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB) 
  
  ```
  
  To see the results I did this
  
  ```
  resultsNames(dds)
   [1] "Intercept"          "Sex_Male_vs_Female" "Age"                "TMB"                "WBC"                "BM_percentage"      "FAB_M1_vs_M0"      
   [8] "FAB_M2_vs_M0"       "FAB_M3_vs_M0"       "FAB_M4_vs_M0"       "FAB_M5_vs_M0"  
  
  ```
  
  I extracted only the Age from the first design
  
  **Design 2**
  
  I again used the design using `Age` variable :
  
  ```
  dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
            colData=coldata, design= ~ Age)
  
  ```
  
  As an example I take this gene this output from the first design
  
  ```
    gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
    <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
  1 ENSG00000196188     97.4         -0.742  0.113 -6.58 4.57e-11 0.00000125
  
  ```
  
  This from the second design:
  
  ```
    gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
    <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
    ENSG00000196188     97.4         -0.757  0.119 -6.34 2.28e-10   3.43e-07
          
  
  ```
  
  Another example this gene from the first design; the result is not significant:
  
  ```
             gene  baseMean log2FoldChange  lfcSE    stat pvalue   padj
  ENSG00000173826  74.54467     0.01283834  0.589  0.0218  0.983  0.993
  
  ```
  
  The same gene from **Design 2** *is* significant:
  
  ```
             gene  baseMean log2FoldChange  lfcSE    stat   pvalue     padj
  ENSG00000173826 74.54467       -2.004743  0.190   -10.6 4.76e-26 1.35e-21
  
  ```
  
  Is this drastic outcome is due to the maths behind it or am I doing something wrong?
  
  So would my interpretation be as such when i consider Age as my only variable I can consider `ENSG00000173826`significant ? But it is not significant when my main effect which i want to test is FAB which is my design 1.
  
  I would like to know how to interpret and report the same.
  
  Any suggestion or help would be really really appreciated
  
  **Tags:** deseq2
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2022-06-23. Score: 2 | Accepted answer
  
  It's not surprising to me that the adjusted fold change statistic changes in different ways for different genes when covariates are added or removed. That's sort of the point of using covariates, as long as it makes sense to include them (which is an experimental decision, not a mathematical or statistical decision).
  
  One thing I'd recommend checking is to make sure the results you're looking at are what you *think* you are looking at. You say that you're interested in testing `FAB`, but that is not included in the second (Age-only) design, which means that it won't be possible to extract any `FAB` comparisons from the results. It looks like your first design might be looking at `Sex` (based on the resultsNames ordering) or one of the `FAB` comparisons (based on reverse ordering of variables in the model), whereas the second might be looking at `Age`. This should be explicitly stated when displaying the results, but you haven't included that information in your output.
  
  If you want `FAB` comparisons in the results, that variable needs to be in the specified model, e.g. "design = ~ Age + FAB". The Design 2 model that you have in your question only includes `Age`.
  
  A good way to do this without getting the wrong contrast is to fully specify the variables when fetching results, as explained in the [DESeq2 vignette](https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs):
  
  > In any case, the `contrast` argument of the function results takes a character vector of length three: the name of the variable, the name of the factor level for the numerator of the log2 ratio, and the name of the factor level for the denominator. The contrast argument can also take other forms, as described in the help page for results and below
  
  ```
  res.age <- results(dds, name="Age")
  res.FAB.M3.M0 <- results(dds, contrast=c("FAB","M3","M0"))
  
  ```
  
  > **Comments:**
  > **username_0**: `design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB) ` I used this one this resource https://github.com/hbctraining/DGE_workshop_salmon_online/blob/master/homework/DGE_assignment_2_answer_key.R  "which means that it won't be possible to extract any FAB comparisons from the results" i could extract the result for `FAB`  
  > **username_0**: this biostar post https://www.biostars.org/p/278684/ explains i hope im getting it right  
  > **username_0**: We can account for the different types of sequencing, and get a clearer picture of the differences attributable to the treatment. As condition is the variable of interest, we put it at the end of the formula. I followed this ..  
  > **username_2**: What username_1 means is you can't get FAB results from design2, because it doesn't include FAB at all.  If you specify the contrast when getting results, **the order of the elements in the design does not matter at all**.   But since you won't say how you called results, no one really knows what you are doing.  
  > **username_2**: Also, unless you have a huge number of samples, I am skeptical that you can really model so many different variables at once.  I would keep only those elements which correlate to the first few PCs.  
  > **username_0**: @username_2 yes I tried the second suggestion so if i want to achieve that 80% variance threshold it would be like 40 components may I need to filter or group the gene in better way...  
  
  ---
  
  