Dear list, I am analyzing a multi-factor time course RNA-Seq experiment with batch-effects using edgeR. After I perform a glmLRT, I cannot use topTags to display the most significant DE tags if I passed contrasts. Here's what I did: [....] # fit is the value from a previous glmFit, design is the value from a model.matrix call # 1) We first call by using a coef > lrt1<-glmLRT(fit,coef=10) > topTags(lrt1) #Work as expected [....] # 2) We call with a constrast > cr<-makeContrasts(time48="TreatA.Time48-TreatB.Time48",levels=design) > lrt2 <-glmLRT(fit,contrast=cr) > topTags(lrt2) Error in abs(object$table$logFC) : Non-numeric argument to mathematical function I don't know if this is the expected behavior (I would say no). I browsed a bit the edgeR source code. Here's what I believe might the culprit : In glmFit.R : [....]#In glmLRT <- fun... # line 161, when using contrasts coef <- 1:ncontrasts if(ncontrasts < ncol(contrast)) contrast <- contrast[,qrc$pivot[coef]] ** logFC <- (glmfit$coefficients %*% contrast)/log(2) if(ncontrasts>1) { [....] # Line 191 tab <- data.frame( *** logFC=logFC, logCPM=(glmfit$abundance+log(1e6))/log(2), [....] logFC is matrix, inheriting its column name from the matrix product at ** (if only one contrast was passed, will be the name of the contrast). Then, when building the value data frame at ***, the name given in the LHS of the assignment is ignored (because logFC is a matrix). The resulting tab has no logFC column. One possible fix would be to cast logFC as a vector when there's only one contrasts. A simple work-around is to rename the first column of the returned value to "logFC", that makes topTags happy. I don't think this bug is data-dependant. Best wishes, Hayssam. PS: Dear Gordon Smyth, thanks a lot for the invaluable help provided on this list! And for edgeR :) > sessionInfo() R version 2.15.1 (2012-06-22) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] grid stats graphics grDevices utils datasets methods [8] base other attached packages: [1] gplots_2.11.0 MASS_7.3-22 KernSmooth_2.23-8 [4] caTools_1.13 bitops_1.0-4.2 gtools_2.7.0 [7] ggplot2_0.9.3 DESeq_1.10.1 lattice_0.20-10 [10] locfit_1.5-8 Biobase_2.18.0 BiocGenerics_0.4.0 [13] edgeR_3.0.6 limma_3.14.3 gdata_2.12.0 [16] stringr_0.6.2 reshape2_1.2.2 data.table_1.8.6 -- Hayssam Soueidan, Ph.D Bioinformatics and Statistics The Netherlands Cancer Institute Plesmanlaan 121, 1066 CX Amsterdam The Netherlands

Dear Hayssam, Yes, thanks, I managed to introduce this bug in edgeR 3.0.6 a week ago. Ryan Thompson alerted me to the same bug a few days ago. It was fixed in edgeR 3.0.7 last Friday. If you update edgeR, the bug should be gone. Best wishes Gordon > Date: Sun, 16 Dec 2012 04:16:15 +0100 > From: <h.soueidan at="" nki.nl=""> > To: <bioconductor at="" r-project.org=""> > Subject: [BioC] edgeR 3.0.6: Possible (small) bug in glmLRT when > passing contrasts? > Message-ID: <b3c23afb-76d7-4c9e-b98d-22eef45d9ef6 at="" nki.nl=""> > Content-Type: text/plain; charset="us-ascii" > > Dear list, > > I am analyzing a multi-factor time course RNA-Seq experiment with > batch-effects using edgeR. After I perform a glmLRT, I cannot use > topTags to display the most significant DE tags if I passed contrasts. > Here's what I did: > > [....] > # fit is the value from a previous glmFit, design is the value from a model.matrix call > > # 1) We first call by using a coef >> lrt1<-glmLRT(fit,coef=10) >> topTags(lrt1) #Work as expected > [....] > > # 2) We call with a constrast >> cr<-makeContrasts(time48="TreatA.Time48-TreatB.Time48",levels=design) >> lrt2 <-glmLRT(fit,contrast=cr) >> topTags(lrt2) > > Error in abs(object$table$logFC) : > Non-numeric argument to mathematical function > > I don't know if this is the expected behavior (I would say no). > > > I browsed a bit the edgeR source code. Here's what I believe might the culprit : > > In glmFit.R : > > [....]#In glmLRT <- fun... > > # line 161, when using contrasts > > coef <- 1:ncontrasts > if(ncontrasts < ncol(contrast)) contrast <- contrast[,qrc$pivot[coef]] > ** logFC <- (glmfit$coefficients %*% contrast)/log(2) > if(ncontrasts>1) { > [....] > > # Line 191 > tab <- data.frame( > *** logFC=logFC, > logCPM=(glmfit$abundance+log(1e6))/log(2), > [....] > > > logFC is matrix, inheriting its column name from the matrix product at ** (if only one contrast was passed, will be the name of the contrast). Then, when building the value data frame at ***, the name given in the LHS of the assignment is ignored (because logFC is a matrix). The resulting tab has no logFC column. > > One possible fix would be to cast logFC as a vector when there's only one contrasts. > > A simple work-around is to rename the first column of the returned value to "logFC", that makes topTags happy. > > I don't think this bug is data-dependant. > > Best wishes, > Hayssam. > PS: Dear Gordon Smyth, thanks a lot for the invaluable help provided on this list! And for edgeR :) > > > >> sessionInfo() > R version 2.15.1 (2012-06-22) > Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) > > locale: > [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] grid stats graphics grDevices utils datasets methods > [8] base > > other attached packages: > [1] gplots_2.11.0 MASS_7.3-22 KernSmooth_2.23-8 > [4] caTools_1.13 bitops_1.0-4.2 gtools_2.7.0 > [7] ggplot2_0.9.3 DESeq_1.10.1 lattice_0.20-10 > [10] locfit_1.5-8 Biobase_2.18.0 BiocGenerics_0.4.0 > [13] edgeR_3.0.6 limma_3.14.3 gdata_2.12.0 > [16] stringr_0.6.2 reshape2_1.2.2 data.table_1.8.6 > > > > -- > Hayssam Soueidan, Ph.D > Bioinformatics and Statistics > The Netherlands Cancer Institute > Plesmanlaan 121, 1066 CX Amsterdam > The Netherlands > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}