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 .......
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():
Setting lambda=0 yields the Benjamini-Hochberg estimator.