[BioC] DEXseq: fitDispersionFunction error
Elena Sorokin
sorokin at wisc.edu
Fri Feb 17 16:11:51 CET 2012
Dear Alejandro,
I can't seem to find the development version that you were talking
about. I am running the version available for BioC 2.9:
http://watson.nci.nih.gov/bioc_mirror/packages/2.9/bioc/html/DEXSeq.html
...Could you tell me where to download v 1.1.6 of your software?
Thanks,
Elena
On 2/17/2012 6:11 AM, Alejandro Reyes wrote:
> Dear Elena,
>
> I am glad you like the package!
> We recently fixed this bug, it should be fine if you use the most
> recent development
> version of DEXSeq (1.1.6)
>
> This should be solved by then, if not, could you send me your ecs
> object to give it a
> closer look?
>
> Thanks!
>
> Alejandro
>
>
> Greetings all,
>
> I ran into another problem with the DEXseq program, and I was hoping
> that somebody could help me once more. I successfully ran through a
> full analysis of one time point of my dataset, and was hoping ideally to
> extend the analysis to my full time course, using the GLM test. However,
> I ran into this error when running fitDispersionFunction:
>
> > ecs<- fitDispersionFunction(ecs)
> Warning messages:
> 1: In glmgam.fit(mm, disps[good], start = coefs) :
> Too much damping - convergence tolerance not achievable
> 2: In glmgam.fit(mm, disps[good], start = coefs) :
> Too much damping - convergence tolerance not achievable
>
> I couldn't find the answer in the list archive. Do the DEXseq authors or
> anyone else know what this error means?
>
> I'd be truly grateful for any advice about this - I really like the
> DEXseq package, and I'm hoping to use/reference it in my manuscript if I
> can get it to work! My full script is below.
>
> Sincerely,
> Elena
> -------------------------------------
> > library(DEXSeq)
> > annotationfile = file.path("../DEXSeq_annotations.gff")
> > annotationfile
> [1] "../DEXSeq_annotations.gff"
> > samples<- read.table("PLdesign.txt",header=T,row.names=1)
> > samples
> treatment time
> PL_D1hr_1 vehicle one
> PL_D1hr_2 vehicle one
> PL_U1hr_1 drug one
> PL_U1hr_2 drug one
> PL_D2hr_1 vehicle two
> PL_D2hr_2 vehicle two
> PL_U2hr_1 drug two
> PL_U2hr_2 drug two
> PL_D4hr_1 vehicle four
> PL_D4hr_2 vehicle four
> PL_U4hr_1 drug four
> PL_U4hr_2 drug four
> PL_D6hr_1 vehicle six
> PL_D6hr_2 vehicle six
> PL_U6hr_1 drug six
> PL_U6hr_2 drug six
> > fullFilenames<- list.files(full.names=TRUE,pattern="exon.txt")
> > fullFilenames
> [1] "./PL_1hrD1_exon.txt" "./PL_1hrD2_exon.txt" "./PL_1hrU1_exon.txt"
> [4] "./PL_1hrU2_exon.txt" "./PL_2hrD1_exon.txt" "./PL_2hrD2_exon.txt"
> [7] "./PL_2hrU1_exon.txt" "./PL_2hrU2_exon.txt" "./PL_4hrD1_exon.txt"
> [10] "./PL_4hrD2_exon.txt" "./PL_4hrU1_exon.txt" "./PL_4hrU2_exon.txt"
> [13] "./PL_6hrD1_exon.txt" "./PL_6hrD2_exon.txt" "./PL_6hrU1_exon.txt"
> [16] "./PL_6hrU2_exon.txt"
> > ecs<- read.HTSeqCounts(countfiles = fullFilenames,design =
> samples,flattenedfile = annotationfile)
> # Check out ecs object
> > head(counts(ecs))
> ./PL_1hrD1_exon.txt ./PL_1hrD2_exon.txt ./PL_1hrU1_exon.txt
> 2L52.1:001 12 16 7
> 2L52.1:002 32 22 19
> 2L52.1:003 18 16 14
> 2L52.1:004 19 16 20
> 2L52.1:005 28 22 9
> 2L52.1:006 19 15 22
> ./PL_1hrU2_exon.txt ./PL_2hrD1_exon.txt ./PL_2hrD2_exon.txt
> 2L52.1:001 4 7 4
> 2L52.1:002 19 23 16
> 2L52.1:003 20 15 17
> 2L52.1:004 11 8 17
> 2L52.1:005 17 23 27
> 2L52.1:006 20 19 23
> ./PL_2hrU1_exon.txt ./PL_2hrU2_exon.txt ./PL_4hrD1_exon.txt
> 2L52.1:001 2 4 0
> 2L52.1:002 18 13 13
> 2L52.1:003 15 5 15
> 2L52.1:004 12 4 8
> 2L52.1:005 18 14 13
> 2L52.1:006 15 10 13
> ./PL_4hrD2_exon.txt ./PL_4hrU1_exon.txt ./PL_4hrU2_exon.txt
> 2L52.1:001 2 1 1
> 2L52.1:002 9 6 8
> 2L52.1:003 5 3 17
> 2L52.1:004 2 3 7
> 2L52.1:005 17 4 5
> 2L52.1:006 16 4 6
> ./PL_6hrD1_exon.txt ./PL_6hrD2_exon.txt ./PL_6hrU1_exon.txt
> 2L52.1:001 4 2 0
> 2L52.1:002 20 13 4
> 2L52.1:003 8 7 7
> 2L52.1:004 7 4 6
> 2L52.1:005 14 16 9
> 2L52.1:006 12 13 6
> ./PL_6hrU2_exon.txt
> 2L52.1:001 2
> 2L52.1:002 12
> 2L52.1:003 8
> 2L52.1:004 4
> 2L52.1:005 6
> 2L52.1:006 4
> > design(ecs)
> treatment time
> ./PL_1hrD1_exon.txt vehicle one
> ./PL_1hrD2_exon.txt vehicle one
> ./PL_1hrU1_exon.txt drug one
> ./PL_1hrU2_exon.txt drug one
> ./PL_2hrD1_exon.txt vehicle two
> ./PL_2hrD2_exon.txt vehicle two
> ./PL_2hrU1_exon.txt drug two
> ./PL_2hrU2_exon.txt drug two
> ./PL_4hrD1_exon.txt vehicle four
> ./PL_4hrD2_exon.txt vehicle four
> ./PL_4hrU1_exon.txt drug four
> ./PL_4hrU2_exon.txt drug four
> ./PL_6hrD1_exon.txt vehicle six
> ./PL_6hrD2_exon.txt vehicle six
> ./PL_6hrU1_exon.txt drug six
> ./PL_6hrU2_exon.txt drug six
> > head(fData(ecs))
> geneID exonID testable dispBeforeSharing dispFitted
> dispersion
> 2L52.1:001 2L52.1 E001 NA NA NA NA
> 2L52.1:002 2L52.1 E002 NA NA NA NA
> 2L52.1:003 2L52.1 E003 NA NA NA NA
> 2L52.1:004 2L52.1 E004 NA NA NA NA
> 2L52.1:005 2L52.1 E005 NA NA NA NA
> 2L52.1:006 2L52.1 E006 NA NA NA NA
> pvalue padjust chr start end strand transcripts
> 2L52.1:001 NA NA chrII 1867 1911 + 2L52.1
> 2L52.1:002 NA NA chrII 2506 2694 + 2L52.1
> 2L52.1:003 NA NA chrII 2738 2888 + 2L52.1
> 2L52.1:004 NA NA chrII 2931 3036 + 2L52.1
> 2L52.1:005 NA NA chrII 3406 3552 + 2L52.1
> 2L52.1:006 NA NA chrII 3802 3984 + 2L52.1
> # All in all, Exon Count Set looks good
> > ecs<- estimateSizeFactors(ecs)
> > sizeFactors(ecs)
> ./PL_1hrD1_exon.txt ./PL_1hrD2_exon.txt ./PL_1hrU1_exon.txt
> ./PL_1hrU2_exon.txt
> 1.7543472 1.6354950 1.6969669
> 1.5496025
> ./PL_2hrD1_exon.txt ./PL_2hrD2_exon.txt ./PL_2hrU1_exon.txt
> ./PL_2hrU2_exon.txt
> 1.4932349 1.3392481 1.4051111
> 1.1780430
> ./PL_4hrD1_exon.txt ./PL_4hrD2_exon.txt ./PL_4hrU1_exon.txt
> ./PL_4hrU2_exon.txt
> 0.9026051 0.5879963 0.5206059
> 0.5726255
> ./PL_6hrD1_exon.txt ./PL_6hrD2_exon.txt ./PL_6hrU1_exon.txt
> ./PL_6hrU2_exon.txt
> 0.7345296 0.7536120 0.7763199
> 0.7529934
> > formuladispersion<- count ~ sample + ( exon + treatment ) * time
> > ecs<- estimateDispersions( ecs, formula = formuladispersion,
> nCores=11)
> Estimating Cox-Reid exon dispersion estimates using 11 cores. (Progress
> report: one dot per 100 genes)
> ..........................................................................................................................................................>
>
>
> > /usr/local/lib64/R/library
> > ecs<- fitDispersionFunction(ecs)
> Warning messages:
> 1: In glmgam.fit(mm, disps[good], start = coefs) :
> Too much damping - convergence tolerance not achievable
> 2: In glmgam.fit(mm, disps[good], start = coefs) :
> Too much damping - convergence tolerance not achievable
>
> > sessionInfo()
> R version 2.14.0 (2011-10-31)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=C LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] multicore_0.1-7 DEXSeq_1.0.2 Biobase_2.14.0
>
> loaded via a namespace (and not attached):
> [1] hwriter_1.3 plyr_1.7.1 statmod_1.4.14 stringr_0.6
> tools_2.14.0
> >
>
More information about the Bioconductor
mailing list