Question: DESeq2 lfcShrink using only dds and res arguments
0
gravatar for Michael Steinbaugh
6 months ago by
Harvard
Michael Steinbaugh30 wrote:

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 • 293 views
ADD COMMENTlink modified 6 months ago by Michael Love24k • written 6 months ago by Michael Steinbaugh30
Answer: DESeq2 lfcShrink using only dds and res arguments
1
gravatar for Michael Love
6 months ago by
Michael Love24k
United States
Michael Love24k wrote:

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 COMMENTlink written 6 months ago by Michael Love24k

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 REPLYlink modified 4 months ago • written 4 months ago by keagy.jason0
1

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 REPLYlink written 4 months ago by Michael Love24k

Awesome! I can't believe I overlooked that!

ADD REPLYlink written 4 months ago by keagy.jason0

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

ADD REPLYlink written 4 months ago by Michael Steinbaugh30

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 REPLYlink modified 4 months ago • written 4 months ago by keagy.jason0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 340 users visited in the last hour