Dear Bioconductors,

I'm trying to use edgeR to compute the differential expression from some
RNA-Seq data, but it seems to be generating some odd results.

I'm running edgeR on data on about 6000 genes, across 5 experimental
conditions. When I compute the differential expression for any two of these
conditions edgeR it is returning a 'nan' fold-change for about 1500 out of
the 6000 genes being tested. Amazingly, it is also returning p-values for
these fold-changes, and the p-values cover a range of values from total
insignificant (>0.5) to very significant (<1E-4)! At first I thought if
might be because there is sometimes no signal for a gene in a given
condition, but 1) other cases of this nature produce '-inf' or 'inf'
fold-changes, not 'nan' fold-hanges, and 2) in some cases edgeR is
calculating a 'nan' fold-change for something that has signal in every
replicate of both conditions! I've checked everything I can think of... I
don't get any errors or warnings, the Norm Factors are sensible, the raw
signal is sensible, all the objects look well-formed (i.e. like they contain
all the bits they should contain).... its very confusing and frustrating.

So, am I doing something wrong? Or is there a deeper problem with this
package?

I'm using R v 2.13.1, edgeR 2.2.5 and the commands I'm using are:

> # create the design matrix
> model.groups<-groups
> model.factors<-as.factor(model.groups)
> model<-model.matrix(~model.factors)
>
> # build DGElist and calculate normalization factors
> x=as.data.frame(data)
> rownames(x)=genenames
> y=DGEList(x,group=groups)
> y=calcNormFactors(y)
> str(y)

Formal class 'DGEList' [package "edgeR"] with 1 slots
  ..@ .Data:List of 3
  .. ..$ :'data.frame': 17 obs. of  3 variables:
  .. .. ..$ group       : Factor w/ 5 levels "cond00","cond06",..: 3 3 3 4 4
4 3 3 1 1 ...
  .. .. ..$ lib.size    : num [1:17] 1866963 2364994 1712838 2782920 2780054
...
  .. .. ..$ norm.factors: num [1:17] 0.992 0.978 1.002 1.036 0.972 ...
  .. ..$ : num [1:6351, 1:17] 7 187 44 0 98 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:17] "cond07.rep0020" "cond07.rep0021"
"cond07.rep0022" "cond08.rep0023" ...
  .. ..$ : Named logi [1:6351] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. ..- attr(*, "names")= chr [1:6351] "38536" "38537" "38538" "38539"
...

>
> # estimate dispersion and fit models
> z=estimateGLMCommonDisp(y, design))
> str(z)

Formal class 'DGEList' [package "edgeR"] with 1 slots
  ..@ .Data:List of 4
  .. ..$ :'data.frame': 17 obs. of  3 variables:
  .. .. ..$ group       : Factor w/ 5 levels "cond00","cond06",..: 3 3 3 4 4
4 3 3 1 1 ...
  .. .. ..$ lib.size    : num [1:17] 1866963 2364994 1712838 2782920 2780054
...
  .. .. ..$ norm.factors: num [1:17] 0.992 0.978 1.002 1.036 0.972 ...
  .. ..$ : num [1:6351, 1:17] 7 187 44 0 98 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:17] "cond07.rep0020" "cond07.rep0021"
"cond07.rep0022" "cond08.rep0023" ...
  .. ..$ : Named logi [1:6351] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. ..- attr(*, "names")= chr [1:6351] "38536" "38537" "38538" "38539"
...
  .. ..$ : num 0.101

>
> fit<-glmFit(z,model,dispersion=z$common.dispersion)
> str(fit)

Formal class 'DGEGLM' [package "edgeR"] with 1 slots
  ..@ .Data:List of 12
  .. ..$ : num [1:6351, 1:5] -12.89 -9.22 -10.72 -11.54 -9.56 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:5] "(Intercept)" "model.factorscond06"
"model.factorscond07" "model.factorscond08" ...
  .. ..$ : int [1:6351] 12 12 12 12 12 12 12 12 12 12 ...
  .. ..$ : Named num [1:6351] 23.3 3.45 30.57 25.89 3.07 ...
  .. .. ..- attr(*, "names")= chr [1:6351] "38536" "38537" "38538" "38539"
...
  .. ..$ : num [1:17, 1:5] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:17] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr [1:5] "(Intercept)" "model.factorscond06"
"model.factorscond07" "model.factorscond08" ...
  .. .. ..- attr(*, "assign")= int [1:5] 0 1 1 1 1
  .. .. ..- attr(*, "contrasts")=List of 1
  .. .. .. ..$ model.factors: chr "contr.treatment"
  .. ..$ : num [1:6351, 1:17] 14.4 14.4 14.4 14.4 14.4 ...
  .. ..$ :'data.frame': 17 obs. of  3 variables:
  .. .. ..$ group       : Factor w/ 5 levels "cond00","cond06",..: 3 3 3 4 4
4 3 3 1 1 ...
  .. .. ..$ lib.size    : num [1:17] 1866963 2364994 1712838 2782920 2780054
...
  .. .. ..$ norm.factors: num [1:17] 0.992 0.978 1.002 1.036 0.972 ...
  .. ..$ : NULL
  .. ..$ : num 0.101
  .. ..$ : num [1:17] 1851365 2313681 1716797 2882487 2703330 ...
  .. ..$ : NULL
  .. ..$ : num [1:6351, 1:17] 6.32 160.53 59.13 12.7 105.68 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:17] "cond07.rep0020" "cond07.rep0021"
"cond07.rep0022" "cond08.rep0023" ...
  .. ..$ : Named num [1:6351] -12.67 -9.33 -10.36 -11.66 -9.64 ...
  .. .. ..- attr(*, "names")= chr [1:6351] "38536" "38537" "38538" "38539"
...

>
> # liklihood ratio statistics
> results=glmLRT(z, fit, coef = 4)
> str(results)

Formal class 'DGELRT' [package "edgeR"] with 1 slots
  ..@ .Data:List of 10
  .. ..$ :'data.frame': 17 obs. of  3 variables:
  .. .. ..$ group       : Factor w/ 5 levels "cond00","cond06",..: 3 3 3 4 4
4 3 3 1 1 ...
  .. .. ..$ lib.size    : num [1:17] 1866963 2364994 1712838 2782920 2780054
...
  .. .. ..$ norm.factors: num [1:17] 0.992 0.978 1.002 1.036 0.972 ...
  .. ..$ : Named logi [1:6351] FALSE FALSE FALSE FALSE FALSE FALSE ...
  .. .. ..- attr(*, "names")= chr [1:6351] "38536" "38537" "38538" "38539"
...
  .. ..$ : num 0.101
  .. ..$ :'data.frame': 6351 obs. of  4 variables:
  .. .. ..$ logConc     : num [1:6351] -12.67 -9.33 -10.36 -11.66 -9.64 ...
  .. .. ..$ logFC       : num [1:6351] 0.7057 -0.0157 0.5583 -0.3311 -0.1913
...
  .. .. ..$ LR.statistic: num [1:6351] 1.38671 0.00165 1.83369 0.46732
0.24042 ...
  .. .. ..$ p.value     : num [1:6351] 0.239 0.968 0.176 0.494 0.624 ...
  .. ..$ : num [1:6351, 1:5] -12.89 -9.22 -10.72 -11.54 -9.56 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:5] "(Intercept)" "model.factorscond06"
"model.factorscond07" "model.factorscond08" ...
  .. ..$ : num [1:6351, 1:4] -12.63 -9.22 -10.52 -11.63 -9.62 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:6351] "38536" "38537" "38538" "38539" ...
  .. .. .. ..$ : chr [1:4] "(Intercept)" "model.factorscond07"
"model.factorscond08" "model.factorscond10"
  .. ..$ : num [1:17, 1:5] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:17] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr [1:5] "(Intercept)" "model.factorscond06"
"model.factorscond07" "model.factorscond08" ...
  .. .. ..- attr(*, "assign")= int [1:5] 0 1 1 1 1
  .. .. ..- attr(*, "contrasts")=List of 1
  .. .. .. ..$ model.factors: chr "contr.treatment"
  .. ..$ : num [1:17, 1:4] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:17] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr [1:4] "(Intercept)" "model.factorscond07"
"model.factorscond08" "model.factorscond10"
  .. ..$ : num 0.101
  .. ..$ : chr "model.factorscond06"

These are some examples of the results:

--------------------
edgeR out:

LR.statistic:  0.975075487675
logConc:  -13.3931239073
p.value:  0.323417620577
logFC:  NaN
geneID:  38956

Raw data:

Condition 1:
replicate 1: 2,372,532 good reads, female: 8.0
replicate 2: 3,966,968 good reads, female: no signal
replicate 3: 1,389,571 good reads, male: no signal

Condition 2:
replicate 1: 3,102,608 good reads, male: 5.0
replicate 2: 3,451,983 good reads, male: no signal
replicate 3: 2,892,192 good reads, male: no signal
replicate 4: 3,620,124 good reads, male: 5.0
replicate 5: 2,640,968 good reads, female: no signal
--------------------

--------------------
edgeR out:

LR.statistic:  3.57045101322
logConc:  -13.6684814046
p.value:  0.0588163345523
logFC:  NaN
geneID:  38959

Raw data:

Condition 1:
replicate 1: 2,372,532 good reads, female: 5.0
replicate 2: 3,966,968 good reads, female: 5.0
replicate 3: 1,389,571 good reads, male: no signal

Condition 2:

replicate 1: 3,102,608 good reads, male: no signal
replicate 2: 3,451,983 good reads, male: no signal
replicate 3: 2,892,192 good reads, male: 7.0
replicate 4: 3,620,124 good reads, male: no signal
replicate 5: 2,640,968 good reads, female: no signal
--------------------

--------------------
edgeR out:


LR.statistic:  1.81602638091
logConc:  -11.265720531
p.value:  0.177786996458
logFC:  NaN
geneID:  38965

Raw data:

Condition 1:

replicate 1: 2,372,532 good reads, female: 22.0
replicate 2: 3,966,968 good reads, female: 27.0
replicate 3: 1,389,571 good reads, male: 5.0

Condition 2:

replicate 1: 3,102,608 good reads, male: 26.0
replicate 2: 3,451,983 good reads, male: 26.0
replicate 3: 2,892,192 good reads, male: 9.0
replicate 4: 3,620,124 good reads, male: 36.0
replicate 5: 2,640,968 good reads, female: 55.0
--------------------


-- 
Cheers,

Nick Schurch

Data Analysis Group (The Barton Group),
School of Life Sciences,
University of Dundee,
Dow St,
Dundee,
DD1 5EH,
Scotland,
UK

Tel: +44 1382 388707
Fax: +44 1382 345 893

	[[alternative HTML version deleted]]

