[BioC] edgeR glm fit error

[Gordon K Smyth]

Hi Colin,

Thanks for the nice reproducible data example.

The problem that you see is caused by lack of convergence of the 
interative glm algorithm for some of your transcripts.  Here are two 
possible solutions.

Firstly, your design has the form of a one-way layout, so you could 
analyse your data using classic exact edgeR instead of the glm code.  In 
this style of analysis, you would make pairwise comparisons between the 
time points.  I have tried this one your data, and it works fine:

   > times <- factor(times,levels=c("zero","one","three","six"))
   > d$samples$group <- times
   > d2 <- estimateCommonDisp(d)
   > d2$common.disp
   [1] 0.2359880
   > d2 <- estimateTagwiseDisp(d2)
   Using grid search to estimate tagwise dispersion.
   > topTags(exactTest(d2,c("zero","one"),common=FALSE))
   Comparison of groups: one-zero
           logConc      logFC       PValue          FDR
   19847 -29.88350  40.265106 7.849699e-50 1.458945e-45
   6370  -31.78645  36.459199 8.076247e-37 7.505257e-33
   20014 -31.57373  36.884641 1.569212e-36 7.655628e-33
   13189 -16.39203 -17.244752 1.647612e-36 7.655628e-33
   3387  -31.89198  36.248152 6.213699e-36 2.309756e-32
   6146  -19.46744  11.313569 4.111695e-31 1.273666e-27
   17831 -30.38242 -39.267266 9.207825e-31 2.444809e-27
   791   -12.50480   7.443195 8.829879e-25 2.051402e-21
   11698 -16.87556   7.510317 1.779005e-24 3.673842e-21
   17562 -31.98794 -36.056237 2.586138e-23 4.806596e-20

If you really do need the glm functionality, because you want to do 
something other than pairwise comparisons, then we will have to make the 
glm code work for you.  We will probably never succeed in writing glm 
fitting code that converges for every possible data set, yet for your data 
we can make the functions work if (i) we recognise that the design is a 
oneway layout or (ii) use better starting values.  If you need this, let 
us know and we will send you some newer code.

Best wishes
Gordon


> Date: Sat, 26 Feb 2011 15:08:28 -0500
> From: Colin Maxwell <csm29 at duke.edu>
> To: bioconductor at r-project.org
> Subject: [BioC] edgeR glm fit error
>
> Hello,
> I'm trying to use edgeR to analyze some RNA-seq time series data. I have
> four time points. The first and last time points have three replicates,
> while the middle two have two replicates. The following gives the error I
> get:
>
> require(edgeR)
> counts <- read.csv("http://www.duke.edu/~csm29/counts.csv",row.names=1)
> d <- DGEList(counts)
> d <- calcNormFactors(d)
> d <- d[rowSums(d$counts)>9,]
> times <- rep(c("zero","one","three","six"),c(3,2,2,3))
> design <- model.matrix(~factor(times))
> d <- estimateCRDisp(d, design)
>
> Error in beta[k, ] <- betaj[decr, ] :
> NAs are not allowed in subscripted assignments
>
> traceback()
> 3: mglmLS(y, design, dispersion, offset = offset)
> 2: adjustedProfileLik(spline.disp[i], y.filt, design = design, offset =
> offset.mat.filt +
>      lib.size.mat.filt)
> 1: estimateCRDisp(d, design)
>
> sessionInfo()
> R version 2.12.0 (2010-10-15)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/C/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_2.0.3
>
> loaded via a namespace (and not attached):
> [1] limma_3.6.9
>
> The function seems to be getting hung up on one or two genes when I recover
> the error. However, when I remove those genes from the data, the problem is
> still there. Any help would be much appreciated!
>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}