[BioC] DEXseq: fitDispersionFunction error
Elena Sorokin
sorokin at wisc.edu
Thu Feb 16 21:45:23 CET 2012
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