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
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")
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.
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.
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?
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 addition: Warning message: 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!
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
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 moved your "Answer" to a "Comment".
I would recommend to use apeglm in this case.