Search
Question: EDGE: Persistent lrt() and odp() command pi0est() error due to low p values
0
15 months ago by
joe.d.moore.ii0 wrote:

Hi,

I am trying to analyze for differential gene expression within a time series of RT-qPCR measured genes across three treatments. Unfortunately, when I try to implement either LRT or ODP on my "deSet" object to determine significantly differentially expressed genes, I encounter an error, seemingly due to the fact that all 14 genes are differentially expressed, which isn't too surprising -- we chose those genes for a reason.

More unfortunately, I haven't figured out how to either add an argument into odp() or lrt() to alter lambda and I'm not sure how I would go about starting my analysis at the pi0est() command where I could insert an argument for lambda. In short, my analysis has been stymied. I would really appreciate any help, including any discussion about whether EDGE is actually appropriate for RT-qPCR data and what other programs might be more appropriate.

I have pasted the last few lines of relevant code below. Not sure how to add the relevant files, which I know would be useful. I've provided the output from head() (see below). Apologies for other norms I might be breaking. Please call me out -- I'll do better next time.

Thanks so much for any help.

Joe

de_obj_TimeSeries <- build_study(data = TimeSeries_nrmlzdtorrsA_expr, grp = grp, tme = tme, ind = ind, sampling = "timecourse")
ODP_TimeSeries <- odp(de_obj_TimeSeries)
Null iteration:  100
Error in pi0est(p, ...) :
ERROR: The estimated pi0 <= 0. Check that you have valid p-values or use a different range of lambda.

ExpressionSet Summary

ExpressionSet (storageMode: lockedEnvironment)
assayData: 1 features, 47 samples
element names: exprs
protocolData: none
phenoData
sampleNames: 1 2 ... 47 (47 total)
varLabels: tme grp
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

de Analysis Summary

Total number of arrays: 47
Total number of probes: 1

Biological variables:
Null Model: ~grp + ns(tme, df = 2, intercept = FALSE)
<environment: 0x00000000039043d0>
Full Model: ~grp + ns(tme, df = 2, intercept = FALSE) + (grp):ns(tme, df = 2,
intercept = FALSE)
<environment: 0x00000000039043d0>

Individuals:
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
[1,]    1    4    9   16   25   36   49   64   90   110   132    12    26    42    60    80   102   126   152
[,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
[1,]   180   210   242   276    24    50    78   108   140   174   210   248   288   330   374   420    36    74
[,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47]
[1,]   114   156   200   246   294   344   396   450   506   564

Expression data:
1        2        3        4        5        6        7        8        9       10       11       12
0.021175 0.019531 0.016925 0.053386 0.061850 0.074847 0.067439 0.057970 0.004961 0.003640 0.033114 0.009714
13       14       15       16       17       18       19       20       21       22       23       24
0.025729 0.007059 0.017994 0.018889 0.033362 0.024705 0.019246 0.023468 0.003617 0.002721 0.001987 0.026010
25       26       27       28       29       30       31       32       33       34       35       36
0.029406 0.014558 0.021823 0.018298 0.032143 0.061413 0.029892 0.025503 0.005466 0.003626 0.002578 0.006927
37       38       39       40       41       42       43       44       45       46       47
0.005139 0.002057 0.003750 0.002426 0.002264 0.007822 0.003929 0.004465 0.001253 0.000911 0.001164
....... 
modified 4 months ago • written 15 months ago by joe.d.moore.ii0

Since you have so few tests, I recommend using the lrt() function instead of odp() which will simply performs a t-test or F-test.  All you need to do is add lambda=0 as an argument and it will be passed to qvalue():

lrt(de_obj_TimeSeries, lambda=0)

Setting lambda=0 yields the Benjamini-Hochberg estimator.

1
15 months ago by
United States
James W. MacDonald46k wrote:

Using any software package intended for microarray analysis on a qPCR experiment with just 14 genes is probably not ideal. The basic idea for most of these packages is that you don't have much data for each gene, but you have lots of genes, so you could borrow information from all the other genes to help with the comparison of each individual gene.

In your case you seem to have a fair number of observations for each gene, but not that many genes, so there isn't much information to borrow, and each gene might not need to borrow any information anyway. I suppose you could use the limma package, and use trend = TRUE and robust = TRUE for the eBayes step, but if this were my experiment to analyze I would probably just go all old school and just fit 14 individual models using lm.

Thanks, James. Good to hear from someone who has done this before. Old school here I come.

I do use limma myself for qPCR experiments, and it works well even with ~10 genes. It is a worthwhile improvement over lm().

pi0est(), on the other hand, is never going to be at all reliable with just 10 genes. That's why limma uses the Benjamini-Hochberg method and doesn't estimate pi0.

Good to know. Thanks, Gordon.

0
4 months ago by
United States
Storey, John D.60 wrote:

Since you have so few tests, I recommend using the lrt() function instead of odp() which will simply performs a t-test or F-test.  All you need to do is add lambda=0 as an argument and it will be passed to qvalue():

lrt(de_obj_TimeSeries, lambda=0)


Setting lambda=0 yields the Benjamini-Hochberg estimator.

0
4 months ago by
joe.d.moore.ii0 wrote:

Thanks for weighing in, John.