deseq2 model design : Different gene output

Question

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 ENSG00000173826significant ? 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

Answer 1

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:

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