[BioC] DEXseq: fitDispersionFunction error
Alejandro Reyes
alejandro.reyes at embl.de
Fri Feb 17 13:11:49 CET 2012
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