[BioC] Error in estimateCommonDisp

[Gordon K Smyth]

Dear Xiaowan Lu,

Thank you for providing me with your data offline.  I runs perfectly for 
me with no problems.  My session is given below.

Looking more closely at your code, I can see a number of misconceptions 
about how to use R and edgeR.  For one thing, you are setting the library 
sizes to be the number of rows in your files.  This is incorrect, and 
there is no need to set this parameter yourself anyway.  Inputing grossly 
incorrect information in this way can lead to the error that you saw.

Could you please try using edgeR in the way recommended in the User's 
Guide.

Best wishes
Gordon

> library(edgeR)
> NH00 <- scan("RefSeqNH00",what=1)
Read 23301 items
> NH10 <- scan("RefSeqNH10",what=1)
Read 23301 items
> LC00 <- scan("RefSeqLC00",what=1)
Read 23301 items
> LC10 <- scan("RefSeqLC10",what=1)
Read 23301 items
> counts <- cbind(NH00,LC00,NH10,LC10)
> group <- factor(c("VEH", "VEH", "E2", "E2"))
> dge <- DGEList(counts=counts,group=group)
> dge <- calcNormFactors(dge)
> dge <- estimateCommonDisp(dge,verbose=TRUE)
Disp = 0.02692 , BCV = 0.1641

> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: i386-w64-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252
[3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C
[5] LC_TIME=English_Australia.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] edgeR_3.4.0  limma_3.18.0



On Fri, 18 Oct 2013, Gordon K Smyth wrote:

> Dear Xiaowan Lu,
>
> The estimateCommonDisp() function in edgeR is used a great deal by us and 
> others without giving an error, so it is likely that there is some special 
> feature of your data that is causing a problem.  Although you have given some 
> code in your email, it doesn't provide us any way to reproduce the error you 
> have or to track down the cause.
>
> First, the error message that you give could not have arisen from the code 
> given in your email.  The error message has arisen from the 
> edgeRFindChangedGenes() function in the GROseq package, which is commented 
> out in the code you give.
>
> It appears that you are not using edgeR directly but are instead using the 
> GROseq package.  GROseq is not available from either Bioconductor or CRAN. If 
> you get errors from GROseq functions, you should write directly to the 
> authors of that package.
>
> If you want to pursue this further you could email your data and edgeR code 
> offline to me or to one of the other edgeR authors.  However you need to 
> provide us with a data example that we can run ourselves using edgeR 
> commands.  We cannot debug functions in the GROseq package.
>
> Best wishes
> Gordon
>
> ---------------------------------------------
> Professor Gordon K Smyth,
> Bioinformatics Division,
> Walter and Eliza Hall Institute of Medical Research,
> 1G Royal Parade, Parkville, Vic 3052, Australia.
> http://www.statsci.org/smyth
>
>
>> Date: Fri, 11 Oct 2013 11:28:52 -0400
>> From: xiaowan.lu at uhnresearch.ca
>> To: <bioconductor at r-project.org>
>> Subject: [BioC] Error in estimateCommonDisp
>> Message-ID: <FNVshYDb.1381505332.8512780.xiaowan.lu at uhnresearch.ca>
>> 
>> To whom it may concern,
>> 
>> I am trying to use the function estimateCommonDisp(object) in edgeR package 
>> to run some analysis, but I keep getting the error:
>> 
>> Error in if (any(input.mean < 0)) stop("input.mean must be non-negative") :
>>  missing value where TRUE/FALSE needed
>> Calls: edgeRFindChangedGenes ... estimateCommonDisp -> equalizeLibSizes -> 
>> q2qnbinom
>> 
>> I double checked the object I used for the function fun and it looked fine 
>> for me. I looked at input.mean, it returns whole bunch of NAs, I am not 
>> sure what went wrong, please help.
>> 
>> Here I attached the code I used:
>> 
>> require(GROseq)
>> require(edgeR)
>> source("FitModel.R")
>> Gbed <- read.table("FinalTranscripts.bed")
>> tID <- paste(Gbed[[1]],Gbed[[2]],Gbed[[6]], sep="")
>> Gbed <- data.frame(Chr=as.character(Gbed[[1]]), 
>> Start=as.integer(Gbed[[2]]), End=as.integer(Gbed[[3]]), 
>> Str=as.character(Gbed[[6]]), RefSeqID=as.character(tID), 
>> MGI=as.character(tID))
>> 
>> G <- LimitToXkb(Gbed, SIZE=13000)
>> 
>> 
>> load("SOAP.RData")
>> nrNH00 <- NROW(NH00)
>> nrNH10 <- NROW(NH10)
>> nrNH40 <- NROW(NH40)
>> nrNH160<- NROW(NH160)
>> nrLC00 <- NROW(LC00)
>> nrLC10 <- NROW(LC10)
>> nrLC40 <- NROW(LC40)
>> nrLC160<- NROW(LC160)
>> 
>> RefSeqNH00 <- CountReadsInInterval(f= G, p= NH00[,c(1:3,6)])
>> RefSeqLC00 <- CountReadsInInterval(f= G, p= LC00[,c(1:3,6)])
>> RefSeqNH10 <- CountReadsInInterval(f= G, p= NH10[,c(1:3,6)])
>> RefSeqLC10 <- CountReadsInInterval(f= G, p= LC10[,c(1:3,6)])
>> 
>> #ChangedGenes <- edgeRFindChangedGenes(RefSeqNH00, nrNH00, 
>> RefSeqLC00,nrLC00, RefSeqNH10,nrNH10, RefSeqLC10,nrLC10, 
>> FILENAME="T/T.ModelFit-10m.png", G=Gbed)
>> 
>> nZI <- !(RefSeqNH00 == 0 & RefSeqLC00 == 0 & RefSeqNH10 == 0 & RefSeqLC10 
>> == 0)
>> cat("nZI=", nZI, "/n")
>> 
>> d <- list()
>> d$counts <- as.matrix(data.frame(NH00= RefSeqNH00[nZI], LC00= 
>> RefSeqLC00[nZI], NHE2= RefSeqNH10[nZI], LCE2= RefSeqLC10[nZI]))
>> dim(d$counts) <- c(NROW(d$counts), NCOL(d$counts))
>> 
>> colnames(d$counts) <- c("NH00", "LC00", "NHE2", "LCE2")
>> 
>> d$samples$files <- c("NH00", "LC00", "NHE2", "LCE2")
>> d$samples$group <- as.factor(c("VEH", "VEH", "E2", "E2"))
>> d$samples$description <- c("VEH", "VEH", "E2", "E2")
>> d$samples$lib.size <- c(nrNH00, nrLC00, nrNH10, nrLC10)
>> d <- new("DGEList",d)
>> d <- estimateCommonDisp(d)
>> 
>> 
>> Xiaowan Lu
>> Bioinformatics Research Assistant,
>> CO-OP Student
>> University Health Network
>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}