Dear edgeR authors and users,
I am using the currently-released version of edgeR, edgeR_2.6.0. When I tried to estimate common dispersions and trended dispersions for my RNA-seq data, I found the errors. I am listing them as below.  Any suggestions are more than welcome!
### cnts contains the count table, info contains the sample information, and libsize is the sizes of libraries
> de <- DGEList(cnts,group=info$type,lib.size=libsize[,4])
> dim(de)
[1] 22205    46

### Normalization
> de <- calcNormFactors(de)

### Filtering out the genes of which the expression levels are low for at least 45 out ### of 46 samples
> sig.thr <- 1
> idx <- !rowSums(cpm(de)<1)>= nrow(info)-sig.thr
> de <- de[idx,]
> dim(de)
[1] 16710    46

### Set up a multi-factor model and design matrix is a full rank matrix
> type <- info$type
> batch <- info$batch
> design <- model.matrix(~batch+type)
> is.fullrank(design)
[1] TRUE
> design
   (Intercept) batchT typeN typeP typeX
1            1      0     0     0     1
2            1      0     0     0     1
3            1      0     0     0     1
4            1      0     0     0     1
5            1      0     0     0     1
6            1      0     0     0     1
7            1      0     0     0     1
8            1      0     0     1     0
9            1      0     0     0     0
10           1      0     0     1     0
11           1      0     0     1     0
12           1      0     1     0     0
13           1      0     1     0     0
14           1      0     1     0     0
15           1      1     0     0     1
16           1      1     0     0     1
17           1      1     0     0     1
18           1      1     0     0     1
19           1      1     0     0     1
20           1      1     0     0     1
21           1      1     0     0     0
22           1      1     0     0     0
23           1      1     0     0     0
24           1      1     0     0     1
25           1      1     0     0     0
26           1      1     0     0     0
27           1      1     0     0     0
28           1      1     0     0     0
29           1      1     0     0     0
30           1      1     0     0     0
31           1      1     0     1     0
32           1      1     0     1     0
33           1      1     0     1     0
34           1      1     0     1     0
35           1      1     0     1     0
36           1      1     0     1     0
37           1      1     0     1     0
38           1      1     0     1     0
39           1      1     0     1     0
40           1      1     0     1     0
41           1      1     0     1     0
42           1      1     0     1     0
43           1      1     0     1     0
44           1      1     0     1     0
45           1      1     0     1     0
46           1      1     0     1     0
attr(,"assign")
[1] 0 1 2 2 2
attr(,"contrasts")
attr(,"contrasts")$batch
[1] "contr.treatment"

attr(,"contrasts")$type
[1] "contr.treatment"

### the ERROR for estimateGLMCommonDisp
>  estimateGLMCommonDisp(de,design)
Error in beta[k, ] <- betaj[decr, ] :
  NAs are not allowed in subscripted assignments

> traceback()
12: mglmLS(y, design = design, dispersion = dispersion, start = start,
        offset = offset, ...)
11: glmFit.default(y, design = design, dispersion = dispersion, offset = offset)
10: glmFit(y, design = design, dispersion = dispersion, offset = offset)
9: adjustedProfileLik(par^4, y, design, offset)
8: f(arg, ...)
7: function (arg)
   -f(arg, ...)(0.540181513475453)
6: optimize(f = fun, interval = interval^0.25, y = y, design = design,
       offset = offset, maximum = TRUE, tol = tol)
5: dispCoxReid(y, design = design, offset = offset, ...)
4: estimateGLMCommonDisp.default(y = y$counts, design = design,
       offset = offset, method = method, ...)
3: estimateGLMCommonDisp(y = y$counts, design = design, offset = offset,
       method = method, ...)
2: estimateGLMCommonDisp.DGEList(de, design)
1: estimateGLMCommonDisp(de, design)

### the ERROR for estimateGLMTrendedDisp
> estimateGLMTrendedDisp(de,design)
Error in if (any(decr)) { : missing value where TRUE/FALSE needed
 > traceback()
16: mglmLS(y, design = design, dispersion = dispersion, start = start,
        offset = offset, ...)
15: glmFit.default(y, design = design, dispersion = dispersion, offset = offset)
14: glmFit(y, design = design, dispersion = dispersion, offset = offset)
13: adjustedProfileLik(par^4, y, design, offset)
12: f(arg, ...)
11: function (arg)
    -f(arg, ...)(1.08036302695091)
10: optimize(f = fun, interval = interval^0.25, y = y, design = design,
        offset = offset, maximum = TRUE, tol = tol)
9: dispCoxReid(y, design = design, offset = offset, ...)
8: estimateGLMCommonDisp.default(y[bin, ], design, method = method,
       offset[bin, ], min.row.sum = 0, ...)
7: estimateGLMCommonDisp(y[bin, ], design, method = method, offset[bin,
       ], min.row.sum = 0, ...)
6: binGLMDispersion(y, design, min.n = min.n, offset = offset, method = method.bin,
       ...)
5: dispBinTrend(y, design, offset = offset, method.trend = "spline",
       ...)
4: estimateGLMTrendedDisp.default(y = y$counts, design = design,
       offset = offset, method = method, ...)
3: estimateGLMTrendedDisp(y = y$counts, design = design, offset = offset,
       method = method, ...)
2: estimateGLMTrendedDisp.DGEList(de, design)
1: estimateGLMTrendedDisp(de, design)

### sessionInfo
> sessionInfo()
R version 2.14.0 (2011-10-31)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US       LC_NUMERIC=C         LC_TIME=en_US
 [4] LC_COLLATE=en_US     LC_MONETARY=en_US    LC_MESSAGES=en_US
 [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C
[10] LC_TELEPHONE=C       LC_MEASUREMENT=en_US LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] edgeR_2.6.0  limma_3.10.3

Thanks in advance,

Hui

Hui Yao, Ph.D.
Principal Statistical Analyst
MD Anderson Cancer Center








	[[alternative HTML version deleted]]

