Question

DESeq2 lfcShrink using only dds and res arguments

0

Entering edit mode

Michael Steinbaugh ▴ 90

@mjsteinbaugh

Last seen 5 months ago

Cambridge, MA

Hi Mike,

According to the DESeq2 lfcShrink() documentation, it appears that the function should work when you provide a DESeqDataSet (dds) and corresponding DESeqResults (res) generated using the standard DESeq() %>% results() approach, correct? On BioC 3.8 (v1.22.2), I'm seeing this error return:

lfcShrink(dds = dds, res = res)
# using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
# additional priors are available via the 'type' argument, see ?lfcShrink for details
Error in results(dds.shr, name = coefAlpha, lfcThreshold = lfcThreshold) : 
  object 'coefAlpha' not found
Calls: lfcShrink -> results
Backtrace:
    █
 1. └─DESeq2::lfcShrink(dds = dds, res = res)
 2.   └─DESeq2::results(dds.shr, name = coefAlpha, lfcThreshold = lfcThreshold)

It seems like the coefAlpha value isn't passing through correctly to results() internally. If I run lfcShrink() using a coef or contrast argument, then the function works as expected. Also, it's a little confusing keeping coef here, since results() now only supports a name string (corresponding to coef inresultsNames()) or a contrast vector. This seems to still be the case in v1.23 on GitHub: lfcShrink.R. Should that be updated to use name instead to match results()?

deseq2 • 2.8k views

ADD COMMENT • link updated 5.2 years ago by Michael Love 41k • written 5.2 years ago by Michael Steinbaugh ▴ 90

score 1 · Accepted Answer · 2019-01-31

1

Entering edit mode

Michael Love 41k

@mikelove

Last seen 3 hours ago

United States

Thanks for catching this discrepancy Michael. lfcShrink was developed in a bit of a constrained way, because I wanted to provide 100% backwards compatibility, allowing users to re-generate previous shrinkage estimates, while also moving towards the new estimators which outperform the Normal prior.

I could in theory try to look at the results object and infer the name or contrast that was used upstream, but I think I'll just be more clear that it requires coef or contrast be specified when calling lfcShrink with type="normal", because I need that to perform shrinkage. The reason for allowing the res object to be passed is mostly so users can have full control of the generation of that object, e.g. usage of IHW / independentFiltering, alternative alpha target, etc. I've updated the documentation and checks to make it clear that coef or contrast needs to be specified for the Normal prior.

The reason I renamed name to coef, is that coef allows for a number or a name, while name is strictly a character. I found specifying a number was useful shorthand (and many are familiar with this argument coming from limma's topTable).

ADD COMMENT • link 5.2 years ago Michael Love 41k

0

Entering edit mode

HI Mike,

I get the same error message when I try to use lfcShrink with type="ashr" and previously calculated results if I used a more complex contrast in the call to "results". I've replicated the problem with the "pasilla" dataset used in the DESeq2 vignette. Is there a way to do this, or would I be better off just making a new model to test this particular contrast (i.e., ~ condition + type)?

load pasilla dataset:

library("pasilla")
pasCts <- system.file("extdata",
                      "pasilla_gene_counts.tsv",
                      package="pasilla", mustWork=TRUE)
pasAnno <- system.file("extdata",
                       "pasilla_sample_annotation.csv",
                       package="pasilla", mustWork=TRUE)
cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
coldata <- read.csv(pasAnno, row.names=1)
coldata <- coldata[,c("condition","type")]
rownames(coldata) <- sub("fb", "", rownames(coldata))
cts <- cts[, rownames(coldata)]
library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition + type)

make a grouping factor to allow comparisons between different combinations of condition and type:

dds$group <- factor(paste0(dds$condition, dds$type))
design(dds) <- ~ 0 + group
dds <- DESeq(dds)

test for overall effect of type:

z <- results(dds, 
    contrast = list(c("grouptreatedsingle.read", "grouptreatedpaired.end"), 
    c("groupuntreatedsingle.read", "groupuntreatedpaired.end")), 
    listValues = c(1/2, -1/2))

this works:

summary(z)

this does not:

zz <- lfcShrink(dds, z, type = "ashr")
summary(zz)

ADD REPLY • link 5.0 years ago keagy.jason • 0

1

Entering edit mode

The second argument of lfcShrink is coef. If you call a function in R without naming it, it assumes you are giving them in order of the function definition.

Try adding res= in front of z.

ADD REPLY • link 5.0 years ago Michael Love 41k

0

Entering edit mode

Awesome! I can't believe I overlooked that!

ADD REPLY • link 5.0 years ago keagy.jason • 0

0

Entering edit mode

Jason, try passing in a coef or contrast argument along with the res in your lfcShrink() call.

ADD REPLY • link 5.0 years ago Michael Steinbaugh ▴ 90

0

Entering edit mode

Hi Michael,

I know how to pass a contrast for a simple comparison between two groups.

zz <- lfcShrink(dds, 
    contrast = list(c("grouptreatedsingle.read"), 
    c("groupuntreatedsingle.read")), type = "ashr")

But if I want to do a contrast like I wrote above that requires the use of "listValues", I'm stuck because I don't know how to pass that parameter to "lfcShrink"

z <- lfcShrink(dds, 
    contrast = list(c("grouptreatedsingle.read", "grouptreatedpaired.end"), 
    c("groupuntreatedsingle.read", "groupuntreatedpaired.end")), 
    listValues = c(1/2, -1/2), type = "ashr")
using 'ashr' for LFC shrinkage. If used in published research, please cite:
    Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
    https://doi.org/10.1093/biostatistics/kxw041
Error in is.list(x) : unused argument (listValues = c(1/2, -1/2))

I should also add that this is a bit of a toy example. I'm actually trying to do this with 3 population x 3 treatments design.

ADD REPLY • link 5.0 years ago keagy.jason • 0