I am trying to run edgeR with the GLM on some data that I've already used with the 'classic' edgeR analysis, as I would like to try it with RUVSeq. However, I am getting the following errors when I attempt to use the edgeR estimateGLMCommonDisp function:
If I use estimateCommonDisp instead, it works just fine, so I don't think it's a data issue. I have reinstalled the LAPACK libraries several times and reinstalled edgeR after each time, but the error remains. Error code -10 doesn't tell me much(!) so I'm about at my wit's end.
As Davide says, can you give us the design matrix that causes the problem, and a set of counts for one gene that can reproduce the error? DSYRK is called by DPOTRF for the Cholesky decomposition, so I guess that's where the problem is occurring within mglmLevenberg.
I've roughly narrowed it to the following code that I copied/pasted from the RUVSeq user's guide. I should emphasize that this is not a problem with RUVSeq, as I have another example with a similar issue that doesn't use it at all. However, since you're the one who responded...
# Run RUV adjustment on spike-ins
spikeRUVSet <- RUVg(set, spikes, k=1)
pData(spikeRUVSet)
plotRLE(spikeRUVSet, outline=FALSE, ylim=c(-4, 4), col=colors[x], main='New version')
# Now use the revised matrix for DE analysis (either one causes error below)
#design <- model.matrix(~x + W_1, data=pData(uqRUVSet))
design <- model.matrix(~x + W_1, data=pData(spikeRUVSet))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
# ERROR:
y = estimateGLMCommonDisp(y, design)
I've copied/pasted the tail of my output below. The W_1 values are slightly different, but I don't imagine it is the source of the problem. I didn't do anything with LAPACK until I got this error the first time; then I thought it must be a library version issue, so I reinstalled the latest LAPACK before I reinstalled R and then the R packages. However, I was unaware of these other LAPACK nuances (didn't even know about the --with-lapack flag) so I will read that documentation -- many thanks for the link.
We have been unable to reproduce the problem with estimateGLMCommonDisp on our end. Just to be sure we're talking about the same thing, is this what the design matrix looks like on your machine?
On my system (CentOS 6.6, edgeR 3.10.2, LAPACK 3.5.0), the design matrix constructed with spikeRUVSet gives a common dispersion of 0.74, while that with uqRUVSet gives a common dispersion of 0.76.
Also, your original post suggests that you're installing LAPACK separately from R. Be careful with how you do this - perhaps a reading of the R documentation here may be helpful. In particular:
Please do bear in mind that using --with-lapack is ‘definitely not recommended’: it is provided only because it is necessary on some platforms and because some users want to experiment with claimed performance improvements. Reporting problems where it is used unnecessarily will simply irritate the R helpers.
Hi Aaron,
As Gordon requested, I'm sending you the sessionInfo() output and the
design matrix as attachments. I still haven't tried to address possible
LAPACK install issues as I've got a meeting this morning; I hope to return
to it this afternoon.
Many thanks to you all!
Mark
--
Dr. Mark Rogers
Research Associate
Intelligent Systems Laboratory
Department of Engineering Mathematics
0117 33 15328
On 25 July 2015 at 10:49, Aaron Lun [bioc] <noreply@bioconductor.org> wrote:
> Activity on a post you are following on support.bioconductor.org
>
> User Aaron Lun <https: support.bioconductor.org="" u="" 6732=""/> wrote Comment:
> LAPACK/BLAS error with edgeR GLM
> <https: support.bioconductor.org="" p="" 70249="" #70385="">:
>
> We have been unable to reproduce the problem with on our end. Just to be
> sure we're talking about the same thing, is this what the design matrix
> looks like on your machine?
>
> > design
> (Intercept) xTrt W_1
> Ctl1 1 0 -0.04539413
> Ctl3 1 0 0.50347642
> Ctl5 1 0 0.40575319
> Trt9 1 1 -0.30773479
> Trt11 1 1 -0.68455406
> Trt13 1 1 0.12845337
> attr(,"assign")
> [1] 0 1 2
> attr(,"contrasts")
> attr(,"contrasts")$x
> [1] "contr.treatment"
>
> On my system, the design matrix constructed with spikeRUVSet gives a
> common dispersion of 0.74, while that with uqRUVSet gives a common
> dispersion of 0.76.
>
> Also, your original post suggests that you're installing LAPACK separately
> from R. Be careful with how you do this - perhaps a reading of the R
> documentation here
> <https: cran.r-project.org="" doc="" manuals="" r-admin.html#lapack=""> may be
> helpful. In particular:
>
> *Please do bear in mind that using --with-lapack is ‘definitely not
> recommended’: it is provided only because it is necessary on some platforms
> and because some users want to experiment with claimed performance
> improvements. Reporting problems where it is used unnecessarily will simply
> irritate the R helpers.*
>
> ------------------------------
>
> Post tags: glm, edger, lapack, segfault, ruvseq
>
> You may reply via email or visit
> C: LAPACK/BLAS error with edgeR GLM
>
Here's the output of sessionInfo(), which is mostly the same as your setup:
R version 3.2.1 (2015-06-18)
Platform: x86_64-unknown-linux-gnu (64-bit)
Running under: CentOS release 6.4 (Final)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.10.2 limma_3.24.14
My best guess is that it's a problem with your edgeR or R installations. See if you get better results after re-installing either or both of them. For R itself, I usually compile it fresh from the CRAN source, but depending on how you've set it up, you might need to apt-get remove it or something equivalent before attempting re-installation.
Attachments don't seem to come through when you send it via the support site. It would probably be better to send it to me directly off-list, or to put up a link to DropBox (or whatever site you like to use for file sharing).
My guess is that RUV is actually capturing the biology (i.e. Your genes are not truly negative controls) and this causes a (almost) colinearity in the design. Can you post your design matrix here?
Scanning the system, it's not hard to believe that LAPACK installations are causing a problem. Although I never explicitly installed it myself (until this recent problem arose), I see liblapack.a and liblapack.so files in several areas (see below). It appears the libraries are used by CGAL, scipy, numpy and MATLAB, possibly others.
I will try to track down the source of the problem (prefer not to uninstall these other packages, especially MATLAB), but if anyone sees an obvious workaround, I'd love to hear about it.
Aaron, it looks as if your trick to download the R tarball and compile it myself worked. It's unclear to me why the synaptic package manager yields a different executable; but for now I'm just happy to have it working. Many thanks for your help!
As Davide says, can you give us the design matrix that causes the problem, and a set of counts for one gene that can reproduce the error? DSYRK is called by DPOTRF for the Cholesky decomposition, so I guess that's where the problem is occurring within
mglmLevenberg.Many thanks for your quick reply!
I've roughly narrowed it to the following code that I copied/pasted from the RUVSeq user's guide. I should emphasize that this is not a problem with RUVSeq, as I have another example with a similar issue that doesn't use it at all. However, since you're the one who responded...
library(RUVSeq)
library(zebrafishRNASeq)
library(RColorBrewer)
data(zfGenes)
filter <- apply(zfGenes,1,function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]
x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered), phenoData=data.frame(x, row.names=colnames(filtered)))
uqset <- betweenLaneNormalization(set, which="upper")
# Run RUV adjustment on upper-quantile
uqRUVSet <- RUVg(uqset, spikes, k=1)
pData(uqRUVSet)
colors <- brewer.pal(3, "Set2")
plotRLE(uqRUVSet, outline=FALSE, ylim=c(-4, 4), col=colors[x], main='Old version')
# Run RUV adjustment on spike-ins
spikeRUVSet <- RUVg(set, spikes, k=1)
pData(spikeRUVSet)
plotRLE(spikeRUVSet, outline=FALSE, ylim=c(-4, 4), col=colors[x], main='New version')
# Now use the revised matrix for DE analysis (either one causes error below)
#design <- model.matrix(~x + W_1, data=pData(uqRUVSet))
design <- model.matrix(~x + W_1, data=pData(spikeRUVSet))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
# ERROR:
y = estimateGLMCommonDisp(y, design)
I've copied/pasted the tail of my output below. The W_1 values are slightly different, but I don't imagine it is the source of the problem. I didn't do anything with LAPACK until I got this error the first time; then I thought it must be a library version issue, so I reinstalled the latest LAPACK before I reinstalled R and then the R packages. However, I was unaware of these other LAPACK nuances (didn't even know about the --with-lapack flag) so I will read that documentation -- many thanks for the link.
estimateGLMCommonDisp
Error in mglmLevenberg(y, design = design, dispersion = dispersion, offset = offset, :
BLAS/LAPACK routine 'DSYRK ' gave error code -10
> design
(Intercept) xTrt W_1
Ctl1 1 0 -0.0390975
Ctl3 1 0 0.4088907
Ctl5 1 0 0.4414666
Trt9 1 1 -0.2138095
Trt11 1 1 -0.7527070
Trt13 1 1 0.1552567
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"