Question: DESeq2 LFC shrinkage is for designs with interactions
gravatar for yk.gatc
10 weeks ago by
yk.gatc0 wrote:


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?




ADD COMMENTlink modified 10 weeks ago • written 10 weeks 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.
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?

ADD REPLYlink modified 10 weeks ago • written 10 weeks ago by yk.gatc0

I moved your "Answer" to a "Comment".

I would recommend to use apeglm in this case.

ADD REPLYlink written 10 weeks ago by Michael Love16k
gravatar for Michael Love
10 weeks ago by
Michael Love16k
United States
Michael Love16k 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)
[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 10 weeks ago by Michael Love16k

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 weeks ago • written 8 weeks 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 weeks ago by Michael Love16k

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 5 weeks 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 5 weeks ago by Michael Love16k

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.



ADD REPLYlink written 5 weeks ago by yk.gatc0

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.

ADD REPLYlink written 5 weeks ago by Michael Love16k

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!

ADD REPLYlink modified 5 weeks ago • written 5 weeks ago by yk.gatc0

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?

ADD REPLYlink written 27 days ago by yk.gatc0

This is expected. apeglm will give back an LFC of near 0 if there is essentially no information to estimate the actual LFC. We'll have a manuscript up soon describing the methods.

ADD REPLYlink written 27 days ago by Michael Love16k
Please log in to add an answer.


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