Recently, I use edgeR in R to find DEG, When I loaded my data and estimated the dispersion with estimateCommonDisp, there are always some warnings for my data: estimateCommonDisp(y,verbose=TRUE) Disp = 99.99477 , BCV = 9.9997 There were 50 or more warnings (use warnings() to see the first 50) warnings() Warning messages: 1: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 2: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 3: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 4: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 5: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 6: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 7: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 8: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' 9: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' ??????.. But when I change to estiamteGLMCommonDisp, everything looks good, and here is the results: estimateGLMCommonDisp(y,verbose=TRUE) Disp = 0.69044 , BCV = 0.8309 So, For my data (just have two groups), can I use the estimateGLMCommonDisplay to estimate the dispersion and then use the exactTest to find the DEGs? Thanks and looking forward to your reply. Tongwu Here is my scripts in R: setwd("~/Work/Axeq_RNA_seq/DEG/edgeR/") x<-read.delim("../All_htseqCounts.txt2",head=T,check.names=F,row.names =1) head(x) target<-read.table("../design_type.list",head=T) library(edgeR) ## For R and S ### head(x[,1:17]) y=DGEList(counts=x[,1:17], group=target$Sensitivity[1:17],) dim(y) ## Filtering keep<-rowSums(cpm(y)>1) >=5 y<-y[keep,] dim(y) ## Re-compute the library sizes y$samples$lib.size <- colSums(y$counts) ## Normalizing y<-calcNormFactors(y) y$samples ## Data exploration plotMDS(y) ## Estimating the dispersion y<-estimateGLMCommonDisp(y,verbose=TRUE) ##Disp = 0.69044 , BCV = 0.8309 y <- estimateGLMTrendedDisp(y) y <- estimateGLMTagwiseDisp(y) plotBCV(y) ## Differential expression et<-exactTest(y) top<- topTags(et,n=100) top cpm(y)[rownames(top),] ## total number of DE genes at 5% FDR is given by ## S/R summary(de<-decideTestsDGE(et,adjust.method="BH",p.value=0.05)) ##[,1] ##-1 310 (up-regulaiton in R) ##0 12881 ##1 457 detags<-rownames(y)[as.logical(de)] plotSmear(et,de.tags=detags) abline(h=c(-1,1),col="blue") -- output of sessionInfo(): > sessionInfo() R version 2.15.2 (2012-10-26) 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] stats graphics grDevices utils datasets methods base other attached packages: [1] edgeR_3.0.6 limma_3.14.3 loaded via a namespace (and not attached): [1] tools_2.15.2 -- Sent via the guest posting facility at bioconductor.org.

Hi Tongwu, Right, so that is an absurdly high Dispersion/BCV from estimateCommonDisp(). > So, For my data (just have two groups), can I use the estimateGLMCommonDisplay to estimate the dispersion and then use the exactTest to find the DEGs? Short answer is: I probably wouldn't do that. At least not yet. Let's back up a minute and figure out what is going on -- I guess what is not clear from your code is what your grouping (you say two groups) actually gives. Specifically, what does d$samples look like (i.e. after specifying group=target$Sensitivity[1:17]) ? I notice that you've called estimateGLMCommonDisp() with no design matrix. This means it would treat all 17 samples as replicates and calculate dispersion as if this were a single group. This is probably not what you want. In general, using estimateGLMCommonDisp() with no design matrix in a 2-group situation would tend to overestimate the dispersion, but here it is the opposite, at least relative to estimateCommonDisp() -- I wouldn't expect qCML and Cox-Reid dispersion estimates to be the same, but they should be similar. So, that's strange. I wonder if there are some super highly variable features that kills the qCML estimation? I have not experienced this before, but have you done other checks like an MDS plot ? Looking for highly variable features ? Best regards, Mark On 14.03.2013, at 16:06, Tongwu [guest] <guest at="" bioconductor.org=""> wrote: > > Recently, I use edgeR in R to find DEG, When I loaded my data and estimated the dispersion with estimateCommonDisp, there are always some warnings for my data: > estimateCommonDisp(y,verbose=TRUE) > Disp = 99.99477 , BCV = 9.9997 > There were 50 or more warnings (use warnings() to see the first 50) > warnings() > Warning messages: > 1: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 2: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 3: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 4: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 5: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 6: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 7: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 8: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > 9: In condLogLikDerSize(y, r, der = 0L) : value out of range in 'lgamma' > ??????.. > But when I change to estiamteGLMCommonDisp, everything looks good, and here is the results: > estimateGLMCommonDisp(y,verbose=TRUE) > Disp = 0.69044 , BCV = 0.8309 > So, For my data (just have two groups), can I use the estimateGLMCommonDisplay to estimate the dispersion and then use the exactTest to find the DEGs? > Thanks and looking forward to your reply. > Tongwu > Here is my scripts in R: > setwd("~/Work/Axeq_RNA_seq/DEG/edgeR/") > x<-read.delim("../All_htseqCounts.txt2",head=T,check.names=F,row.nam es=1) > head(x) > target<-read.table("../design_type.list",head=T) > library(edgeR) > ## For R and S ### > head(x[,1:17]) > y=DGEList(counts=x[,1:17], group=target$Sensitivity[1:17],) > dim(y) > ## Filtering > keep<-rowSums(cpm(y)>1) >=5 > y<-y[keep,] > dim(y) > ## Re-compute the library sizes > y$samples$lib.size <- colSums(y$counts) > ## Normalizing > y<-calcNormFactors(y) > y$samples > ## Data exploration > plotMDS(y) > ## Estimating the dispersion > y<-estimateGLMCommonDisp(y,verbose=TRUE) > ##Disp = 0.69044 , BCV = 0.8309 > y <- estimateGLMTrendedDisp(y) > y <- estimateGLMTagwiseDisp(y) > plotBCV(y) > ## Differential expression > et<-exactTest(y) > top<- topTags(et,n=100) > top > cpm(y)[rownames(top),] > ## total number of DE genes at 5% FDR is given by ## S/R > summary(de<-decideTestsDGE(et,adjust.method="BH",p.value=0.05)) > ##[,1] > ##-1 310 (up-regulaiton in R) > ##0 12881 > ##1 457 > detags<-rownames(y)[as.logical(de)] > plotSmear(et,de.tags=detags) > abline(h=c(-1,1),col="blue") > > -- output of sessionInfo(): > >> sessionInfo() > R version 2.15.2 (2012-10-26) > 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] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] edgeR_3.0.6 limma_3.14.3 > > loaded via a namespace (and not attached): > [1] tools_2.15.2 > > -- > Sent via the guest posting facility at bioconductor.org. > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor ---------- Prof. Dr. Mark Robinson Bioinformatics, Institute of Molecular Life Sciences University of Zurich http://tiny.cc/mrobin