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.

> head(de_obj_TimeSeries)

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
  varMetadata: labelDescription
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 
....... 

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.

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.

Thanks for weighing in, John.