[BioC] edgeR Warings
Tongwu [guest]
guest at bioconductor.org
Thu Mar 14 16:06:21 CET 2013
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.
More information about the Bioconductor
mailing list