Description Usage Arguments Details Value Author(s) References See Also Examples
Calculates the classical local false discovery rate for multiple parallel t-statistics.
1 2 |
xdat |
the matrix of expression values, with genes as rows and samples as columns |
grp |
a grouping variable giving the class membership of each sample, i.e. each column in |
test |
a function that takes |
p0 |
if supplied, an estimate for the proportion of non-differentially expressed genes; if not supplied, the routine will estimate it, see Details. |
nperm |
number of permutations for establishing the null distribution of the t-statistic |
nr |
the number of equidistant breaks into which the range of test statistics is broken for calculating the fdr. |
seed |
if specified, the random seed from which the permuations are started |
null |
optional argument for passing in a pre-calculated null distribution, see Details. |
zlim |
if no |
sv2 |
positive number controlling the initial degree of smoothing for the densities involved, with smaller values indicating more smoothing; see Details. |
err |
positive number controlling the convergence of the smoothing procedure, with smaller values implying more iterations; see Details. |
verb |
logical value indicating whether provide extra information. |
... |
extra arguments to function |
This function calculates the local false discovery rate (fdr, as opposed to global FDR) introduced by Efron et al., 2001. The underlying model assumes that for a given grouping of samples, we study a mixture of differentially expressed (DE) and non-DE genes, and that consequently, the observed distribution of test statistics is a mixture of test statistics under the alternative and the null statistic, respectively. The densities involved are estimated nonparametrically and smoothed, using a permutation argument for the null distribution.
By default, the null distribution is generated and stored only within the function, and is not available outside. If someone wants to study the null distribution, or wants to re-use the same null distribution efficiently while e.g. varying the smoothing parameter, the argument null
allows the use of an externally generated null distribution, created e.g. using the PermNull
function.
Theoretically, the function should support any kind of function along the lines of tstatistics
, however, this has not been tested, and the current setup is very much geared towards t-tests.
We use non-Gaussian mixed model smoothing for Possion counts for smoothing the density estimates, see Pawitan, 2001, and smooth1d
.
Basically, a data frame with one row per gene and two columns: tstat
, the test statistic, and fdr.local
, the local false discovery rate. This data frame has the additional class attributes fdr1d.result
and fdr.result
, see Examples. This is the bad old S3 class mechanism employed to provide plot and summary functions.
Additional information is provided by a param
attribute, which is a list with the following entries:
p0 |
the proportion of non-differentially expressed genes used when calculating the fdr. |
p0.est |
a logical value indicating whether |
fdr |
the smoothed fdr values calculated for the original intervals. |
xbreaks |
vector of breaks for the test statistic defining the interval for calculation. |
A. Ploner
Efron B, Tibshirani R, Storey JD, Tusher V (2001) Empirical Bayes Analysis of a Microarray Experiment. JASA, 96(456), p. 1151-60.
Pawitan Y.(2001) In All Likelihood, Oxford University Press, ch. 18.11
plot.fdr1d.result
, summary.fdr.result
, OCshow
, tstatistics
, smooth1d
, fdr2d
, PermNull
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | # We simulate a small example with 5 percent regulated genes and
# a rather large effect size
set.seed(2000)
xdat = matrix(rnorm(50000), nrow=1000)
xdat[1:25, 1:25] = xdat[1:25, 1:25] - 1
xdat[26:50, 1:25] = xdat[26:50, 1:25] + 1
grp = rep(c("Sample A","Sample B"), c(25,25))
# A default run
res1d = fdr1d(xdat, grp)
res1d[1:20,]
# Looking at the results
summary(res1d)
plot(res1d)
res1d[res1d$fdr<0.05, ]
# Averaging estimates the global FDR for a set of genes
ndx = abs(res1d$tstat) > 3
mean(res1d$fdr[ndx])
# Extra information
class(res1d)
attr(res1d,"param")
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.