[BioC] edgeR 3.0.6: Possible (small) bug in glmLRT when passing contrasts?

[h.soueidan at nki.nl]

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