Add testing against a threshold for ashr shrinkage

[frederikziebell]

Hi Mike,

I recently found out that it's easy to add testing against a fold-change in ashr and thought this could be added in lfcShrink. Here's an example:

library(airway)
data(airway)

keep <- rowSums(assay(airway) >= 5) >=3
dds <- DESeqDataSet(airway[keep,], design = ~ cell + dex)
dds <- DESeq(dds)

res <- results(dds, name = "dex_untrt_vs_trt")
betahat <- res$log2FoldChange
sebetahat <- res$lfcSE
fit <- ashr::ash(betahat, sebetahat, mixcompdist = "normal", method = "shrink")

lfcThreshold <- log2(1.5)
cdf_pos_lfc <- ashr::cdf_post(fit$fitted_g, lfcThreshold, ashr::set_data(betahat, sebetahat))
cdf_neg_lfc <- ashr::cdf_post(fit$fitted_g, -lfcThreshold, ashr::set_data(betahat, sebetahat))

res$log2FoldChange <- fit$result$PosteriorMean
res$lfcSE <- fit$result$PosteriorSD

res$lfsr <- ifelse(res$log2FoldChange > 0, cdf_pos_lfc, 1 - cdf_neg_lfc)
res <- res[order(res$lfsr),]
res$svalue <- dplyr::cummean(res$lfsr)
res$pvalue <- NULL
res$padj <- NULL

plotMA(lfcShrink(dds, coef="dex_untrt_vs_trt", type="ashr"))
plotMA(res)

Could you maybe add this functionality?

Best,
Frederik

[mikelove]

Thanks Frederik for the suggestion and code example.

I've started working on this code here, it will just take a week or so to cleanup the function and rewrite tests and docs.

https://github.com/mikelove/DESeq2/tree/ashr-lfc

[mikelove]

I just pushed a commit 0f2020b but I haven't tested it yet. My C++ appears to be a bit broken right now (I'm away from office) and I need to link it up to C++11 for Armadillo to work.

[mikelove]

Ok fixed my C++ (just had to add -std=gnu++11 to my .R/Makevars where it had been overwritten by a new install).

It appears to work on my end, if you could double check, then I'll fold it back in to the devel branch.

[frederikziebell]

Checked it and it works. But the line res$svalue <- apeglm::svalue(lfsr) requires to have apeglm installed when only ashr presence is checked in the lfcShrink call. Maybe copy over the code of apeglm::svalue to not force the user to install two packages when only one is really needed.

[mikelove]

Ok, this should work?

26fb6f4

[frederikziebell]

Works for me, perfect!

[mikelove]

Great, pushed to devel branch on Bioc. Thanks!