[BioC] Error in estimateCommonDisp

[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