Can't change coefficient of dds object: error 'coef' should specify same coefficient as in results 'res'

Question

[macielrodriguez2]

Good evening :)

I'm completely new to RNA-seq and, because I'm going to do that analysis in the coming months, I was following the RNA-seq with Bioconductor course in Datacamp. I'm in the part where we going to run results() on my dds object.

My two conditions are normal and fibrosis and I want normal to be the reference or base level.

Following the course, I should specify this in the contrast argument of results() in the following way:

smoc2_res <- results(object = dds, contrast = c("condition of interest", "level to compare", "base or reference level"), alpha = 0.05, 
                 lfcThreshold = 0.32)

So, I wrote:

smoc2_res <- results(object = dds, contrast = c("condition", "fibrosis", "normal"), alpha = 0.05, 
                 lfcThreshold = 0.32)

And looking at the first line in the result with head(smoc2_res):

log2 fold change (MLE): condition fibrosis vs normal 

I got what I specified, the comparison is between fibrosis and normal, being normal the reference level.

As this guide says, to shrink the log fold changes, I have to "provide the dds object and the name or number of the coefficient I want to shrink, where the number refers to the order of the coefficient as it appears in resultsNames(dds)".

Running resultsNames(dds) to get the coefficient I got:

[1] "Intercept"                    "condition_normal_vs_fibrosis"

But, as normal is the reference level, I was expecting "condition_fibrosis_vs_normal" to be present in resultsNames(dds).

Following some posts like

• type='apeglm' shrinkage only for use with 'coef'
• Could not shrink after adjusting factor level reference in DESeq2
• http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#factorlevels

I tried to change the reference level of the dds object with:

dds$condition <- factor(x = dds$condition, levels = c("normal", "fibrosis"))

I was expecting "condition_fibrosis_vs_normal" but checking with resultsNames(dds) I got again:

[1] "Intercept"                    "condition_normal_vs_fibrosis"

Trying with relevel():

dds$condition <- relevel(dds$condition, ref = "normal")

And checking with resultsNames(dds) i get again:

[1] "Intercept"                    "condition_normal_vs_fibrosis"

I'm must be doing something wrong because the coefficient doesn't change...

And because I can get the coefficients right, when I run:

res_lfc <- lfcShrink(dds = dds, coef = "condition_normal_vs_fibrosis", res = smoc2_res)

As expected, I get an error message:

Error in lfcShrink(dds = dds_smoc2, coef = 2, res = smoc2_res) : 'coef' should specify same coefficient as in results 'res'

or

as normal is the reference level, was expecting condition_fibrosis_vs_normal to be present in 'resultsNames(object)'

Sorry if this is a basic mistake. I'm very new to this analysis so I'm sure there is something important that I'm missing.

I'll be very grateful for your kind help.

Thanks a lot in advance.

Answer 1

[ATpoint]

Don't worry, this is a common pitfalls. You will need to relevel the groups to be able to run apeglm for all contrasts you want, I made a short tutorial here:

type='apeglm' shrinkage only for use with 'coef'

[macielrodriguez2]

Thank you so much for your help!!!

I followed your tutorial and it worked!!! Max happiness!!!!!

I just have one doubt left. In this part

dds <- nbinomWaldTest(dds)
resultsNames(dds)
[1] "Intercept"        "condition_A_vs_C" "condition_B_vs_C"

We only the Wald test to get new MLE coefficients for the comparison of interest.

Why is that? the other elements stored in the dds object are not affected in any way by the change in level of the condition? in such a way that we only need to re run the Wald test?

Thanks a lot for your help!!!

[ATpoint]

No, the other relevant entries of dds are not affected if you only change the level of the group. Normalization stays the same as it is the exact same sample (and norm. is independent of the design), and dispersion stays the same as you do not change the covariates but only change the reference level. There is a post from the DESeq2 maintainer suggesting exactly what I wrote in the tutorial (this is how I learned how to run that code), I did not figure that out myself. In the vignette you can read about this https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html, search for the paragraph starting with:

Although apeglm cannot be used with contrast, we note that many designs can be easily rearranged such that what was a contrast becomes its own coefficient.

[macielrodriguez2]

No, the other relevant entries of dds are not affected if you only change the level of the group. Normalization stays the same as it is the exact same sample (and norm. is independent of the design), and dispersion stays the same as you do not change the covariates but only change the reference level.

Thank you so much for all your kind help and explanations ATpoint! I understand better now!

I'll thoroughly read the DESEq2 vignette, especially the part that you are pointing me.

Thank you so much for all your help and patience :)