Search
Question: DESeq2 LFC shrinkage is for designs with interactions
0
8 months ago by
yk.gatc0
yk.gatc0 wrote:

Hi,

It seems like LFC shrinkage is not implemented for designs with interactions. Is there a way to get shrinked LFC values for designs with interactions?

DESeq2_1.19.14

Thanks,

YK

modified 8 months ago • written 8 months ago by yk.gatc0

Thanks Michael for the quick response.

I did try "ashr" method, but it gave me the following warning -

> resultsNames(dds)
[1] "Intercept"    "grp_Y_vs_X"   "grpX.ind.n10" "grpY.ind.n10" "grpX.ind.n2"
[6] "grpY.ind.n2"  "grpX.ind.n3"  "grpY.ind.n3"  "grpX.ind.n4"  "grpY.ind.n4"
[11] "grpX.ind.n5"  "grpY.ind.n5"  "grpX.ind.n6"  "grpY.ind.n6"  "grpX.ind.n7"
[16] "grpY.ind.n7"  "grpX.ind.n8"  "grpY.ind.n8"  "grpX.ind.n9"  "grpY.ind.n9"
[21] "grpX.cndT"    "grpY.cndT"
> res <- results(dds, contrast=list("grpY.cndT","grpX.cndT"))

> resLFC <- lfcShrink(dds,type="ashr", res=res)
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
Due to absence of package REBayes, switching to EM algorithm
Warning message:
In estimate_mixprop(data, g, prior, optmethod = optmethod, control = control) :
Optimization failed to converge. Results may be unreliable. Try increasing maxiter and rerunning.

Any thoughts?

I would recommend to use apeglm in this case.

1
8 months ago by
Michael Love19k
United States
Michael Love19k wrote:

If you check ?lfcShrink in current release it says,

     For ‘type="normal"’, shrinkage cannot be applied to coefficients
in a model with interaction terms.

So for >= v1.18 with apeglm or ashr shrinkage estimators, you can perform this. We haven't written up a full explanation of these in the vignette, it has a section that I'm looking to expand, but haven't had the time yet. Here is some demo code though for an interaction term:

dds <- makeExampleDESeqDataSet()
dds$x <- factor(rep(1:2,6)) design(dds) <- ~x + condition + x:condition dds <- DESeq(dds, quiet=TRUE) resultsNames(dds) [1] "Intercept" "x_2_vs_1" "condition_B_vs_A" "x2.conditionB" ​res1 <- lfcShrink(dds, coef=4, type="apeglm") res2 <- lfcShrink(dds, coef=4, type="ashr") ADD COMMENTlink written 8 months ago by Michael Love19k Thanks for the suggestions. But with the design I am using, I cannot use coef argument with type="apeglm" to shrink LFC values. My objective here is to compare two treatment groups (with matched controls) which I supply as contrast=list(condA,condB). Do you have any suggestions to apply LFC shrinkage? I see LFC values way overestimated upto 30 folds without betaPrior? Please see the attached image. ADD REPLYlink modified 8 months ago • written 8 months ago by yk.gatc0 I see, you can use apeglm here by re-forming the design. I'm guessing your current design is something like ~group + group:individual + group:condition You can use: ~group + condition + group:individual + group:condiiton Which will give you one interaction term, "grpY.cndT", that now represents the difference between the groups in the condition effect. You can use apeglm to shrink this effect. You actually don't need to re-estimate the dispersion, because these designs are equivalent in terms of the estimated counts, only different in terms of the meaning of coefficients. ADD REPLYlink written 8 months ago by Michael Love19k Thank you Michael for your time and support. Btw, I am getting an error when I try to apply apeglm method - res2 <- lfcShrink(dds, coef = 22, type = "apeglm", res = res) using 'apeglm' for LFC shrinkage Error in apeglm::apeglm(Y = Y, x = design, log.lik = log.lik, param = disps, : unused argument (method = apeMethod) [1] apeglm_1.0.3 DESeq2_1.19.34 Any suggestions? ADD REPLYlink written 7 months ago by yk.gatc0 You are mixing versions. That is devel DESeq2 and release apeglm. Bioconductor only works when all packages are installed via biocLite. ADD REPLYlink written 7 months ago by Michael Love19k I reinstalled the packages from Bioconductor, now I get different errors - > resLFC <- lfcShrink(dds, coef=22, res=res, type="apeglm") using 'apeglm' for LFC shrinkage Error in solve.default(o$hessian) :
system is computationally singular: reciprocal condition number = 1.16041e-129
In sqrt(diag(sigma)) : NaNs produced

[1] DESeq2_1.18.1              SummarizedExperiment_1.8.1
[3] DelayedArray_0.4.1         matrixStats_0.53.1
[5] Biobase_2.38.0             GenomicRanges_1.30.1
[7] GenomeInfoDb_1.14.0        IRanges_2.12.0
[9] S4Vectors_0.16.0           BiocGenerics_0.24.0
[11] apeglm_1.0.3               BiocInstaller_1.28.0

Any suggestions? btw, I am running from a docker container.

Thanks for the report of this.

Hmm, I did a lot of work updating the fitting algorithm in apeglm in the devel branch, so that could help. If it's easy for you to switch to devel branch for both DESeq2 and apeglm, I'd be curious if the new code helps.

But also, you can just use type="ashr", and either ignore the warning above or try to install REBayes*. Because you have so many coefficients (at least 22), ashr is faster than DESeq2's or apeglm's shrinkage estimator (which runs a GLM).

* Installing REBayes is not easy. You have to install the mosek software outside of R (a tool for solving mathematical optimization problems), and for that you need to download an academic license, etc. Then you have to install Rmosek. It took me a few hours to get this working.

With the new code (developer version), lfcShrink estimation was done successfully with the following message  -

> resLFC <- lfcShrink(dds, coef=22, type="apeglm", res=res)
using 'apeglm' for LFC shrinkage
Some rows did not converge in finding the MAP

 [1] DESeq2_1.19.36             apeglm_1.1.28

Thank you for the quick responses and timely changes!

Long live Open Source and Open Science!

lfcShrink with "normal" and "apeglm" gives different volcano patterns (please see image attached). Is this expected behavior of different types of shrinkage methods? Or is it something odd here?