Hi All
I am using R 2.14.2 version. I have tried running the goseq tutorial code and it fails to perform the probability weighting function calculation (pwf).
I am not sure what is going on, here is the part I am stuck at:
> library("goseq")
Loading required package: BiasedUrn
Loading required package: geneLenDataBase
> library("edgeR")
Loading required package: limma
> table.summary=read.table(system.file("extdata","Li_sum.txt",package='goseq'),
+ sep='\t',header=TRUE,stringsAsFactors=FALSE)
> counts=table.summary[,-1]
> rownames(counts)=table.summary[,1]
> grp=factor(rep(c("Control","Treated"),times=c(4,3)))
> summarized=DGEList(counts,lib.size=colSums(counts),group=grp)
> disp=estimateCommonDisp(summarized)
> disp$common.dispersion
[1] 0.05688364
> tested=exactTest(disp)
> topTags(tested)
Comparison of groups: Treated-Control
logFC logCPM PValue FDR
ENSG00000127954 11.557868 6.680782 2.574972e-80 1.274766e-75
ENSG00000151503 5.398963 8.499425 1.781732e-65 4.410322e-61
ENSG00000096060 4.897600 9.446635 7.983756e-60 1.317479e-55
ENSG00000091879 5.737627 6.281841 1.207655e-54 1.494654e-50
ENSG00000132437 -5.880436 7.951978 2.950042e-52 2.920896e-48
ENSG00000166451 4.564246 8.458263 7.126763e-52 5.880292e-48
ENSG00000131016 5.254737 6.607187 1.066807e-51 7.544766e-48
ENSG00000163492 7.085400 5.126627 2.716461e-45 1.681014e-41
ENSG00000113594 4.051053 8.602951 9.272066e-44 5.100255e-40
ENSG00000116285 4.108522 7.864828 6.422468e-43 3.179507e-39
> genes=as.integer(p.adjust(tested$table$PValue[tested$table$logFC!=0],
+ method="BH")<.05)
> names(genes)=row.names(tested$table[tested$table$logFC!=0,])
> table(genes)
genes
0 1
19535 3208
> pwf=nullp(genes,"hg19","ensGene")
Loading hg19 length data...
Error in pcls(G) : initial parameters not feasible
I will appreciate any help on this matter.
Rocky
I am using R 2.14.2 version. I have tried running the goseq tutorial code and it fails to perform the probability weighting function calculation (pwf).
I am not sure what is going on, here is the part I am stuck at:
> library("goseq")
Loading required package: BiasedUrn
Loading required package: geneLenDataBase
> library("edgeR")
Loading required package: limma
> table.summary=read.table(system.file("extdata","Li_sum.txt",package='goseq'),
+ sep='\t',header=TRUE,stringsAsFactors=FALSE)
> counts=table.summary[,-1]
> rownames(counts)=table.summary[,1]
> grp=factor(rep(c("Control","Treated"),times=c(4,3)))
> summarized=DGEList(counts,lib.size=colSums(counts),group=grp)
> disp=estimateCommonDisp(summarized)
> disp$common.dispersion
[1] 0.05688364
> tested=exactTest(disp)
> topTags(tested)
Comparison of groups: Treated-Control
logFC logCPM PValue FDR
ENSG00000127954 11.557868 6.680782 2.574972e-80 1.274766e-75
ENSG00000151503 5.398963 8.499425 1.781732e-65 4.410322e-61
ENSG00000096060 4.897600 9.446635 7.983756e-60 1.317479e-55
ENSG00000091879 5.737627 6.281841 1.207655e-54 1.494654e-50
ENSG00000132437 -5.880436 7.951978 2.950042e-52 2.920896e-48
ENSG00000166451 4.564246 8.458263 7.126763e-52 5.880292e-48
ENSG00000131016 5.254737 6.607187 1.066807e-51 7.544766e-48
ENSG00000163492 7.085400 5.126627 2.716461e-45 1.681014e-41
ENSG00000113594 4.051053 8.602951 9.272066e-44 5.100255e-40
ENSG00000116285 4.108522 7.864828 6.422468e-43 3.179507e-39
> genes=as.integer(p.adjust(tested$table$PValue[tested$table$logFC!=0],
+ method="BH")<.05)
> names(genes)=row.names(tested$table[tested$table$logFC!=0,])
> table(genes)
genes
0 1
19535 3208
> pwf=nullp(genes,"hg19","ensGene")
Loading hg19 length data...
Error in pcls(G) : initial parameters not feasible
I will appreciate any help on this matter.
Rocky
Comment