Assume that there are a total m genes under study with P-values of p
i
, i = 1,…,m. From these, we calculate the local FDRs [4–7] and the q-values [3]: fdr
i
and q
i
, for i = 1,…,m, respectively, using false discovery rate analysis package in R, such as fdrtool (specifying statistic = “pvalue”, plot = FALSE). Assume that among them there are a total of r (r > 0) genes with q-values at most as large as 0.05. We declare those genes significant with FDR controlled at 5 % level, and put them in an S set: S = {i : q
i
≤ 0.05}.
As the unit of analysis for an FDR control is a P-value rather than a study subject, we propose a P-value-based bootstrap method to characterize the variability of FDR control. Whereas the usual bootstrap method samples with replacement of the study subjects, our P-value-based bootstrap method samples with replacement directly of the P-values. This is computationally much more efficient, because the P-values in our method do not need to be re-computed from scratch for each bootstrapped sample as in the usual study-subject-based bootstrapping.
To be precise, the j th gene of a bootstrapped sample is G
j
= [m × U + 1], where U is the uniform(0,1) distribution and [x] returns the largest integer not exceeding x. It has a P-value of \({p}_j^{*}={p}_{G_j}.\) From this new set of P-values: p
*
j
for j = 1,…,m., we calculate a new set of local FDRs: fdr
*
j
for j = 1,…,m. Note a star is superscripted to avoid confusion.
There is no guarantee that each and every gene in the original data will be represented in the bootstrapped sample. Put those ‘missing’ genes in a set: M = {i : i ≠ G
j
for j = 1, …, m}. For an i ∉ M, we simply let its bootstrapped local FDR (superscripted B) be fdr
B
i ∉ M
= fdr
*
j
, where j is any value satisfying G
j
=i. For an i ∈ M, we use linear interpolation to estimate its bootstrapped local FDR. First, we find its left and right ‘flanking’ genes. The left flanking genes are those that have the largest P-value (but no larger than p
i
) in the bootstrapped sample, that is, the set: \(\mathrm{L}=\left\{j:{p}_j^{*}=\underset{p_k^{*}\le {p}_i}{ \max}\left({p}_k^{*}\right)\right\}\). The right flanking genes are those that have the smallest P-value (but no smaller than p
i
) in the bootstrapped sample, that is, the set: \(\mathrm{R}=\left\{j:{p}_j^{*}=\underset{p_k^{*}\ge {p}_i}{ \min}\left({p}_k^{*}\right)\right\}\). If L is non-empty, we randomly pick one member in it, say u, and let p
L
= p
*
u
and fdr
L
= fdr
*
u
. If L is empty, we let p
L
= fdr
L
= 0. If R is non-empty, we randomly pick one member in it, say v, and let p
R
= p
*
v
and fdr
R
= fdr
*
v
. If R is empty, we let p
L
= fdr
L
= 1. Now we can use the linear interpolation. If p
L
≠ p
R
, the bootstrapped local FDR for this i ∈ M is \({\mathrm{fd}\mathrm{r}}_{i\in \mathrm{M}}^{\mathrm{B}}=\frac{\mathrm{fd}{\mathrm{r}}_R\times \left({p}_k-{p}_L\right)+\mathrm{f}\mathrm{d}{\mathrm{r}}_L\times \left({p}_R-{p}_k\right)}{p_R-{p}_L}\). If p
L
= p
R
, we let fdr
B
i ∈ M
= fdr
R
(fdr
L
= fdr
R
in this situation anyway).
In a bootstrapped sample, we calculate the bootstrapped q-value by simply averaging the bootstrapped local FDRs pertaining to the r significant genes, that is, \({q}^{\mathrm{B}}=\frac{1}{r}\times {\displaystyle \sum_{i\in \mathrm{S}}{\mathrm{fdr}}_i^{\mathrm{B}}}\). Next, we simulate a binary ‘false discovery indicator’ (1: false positive; 0: true positive) for each and every significant gene. The simulation is done according to an independent Bernoulli distribution with the corresponding bootstrapped local FDR as the parameter. The bootstrapped total number of false positives is then simply the summation of these false discovery indicators, and the bootstrapped false discovery proportion (FDP), that number divided by r, that is, \({\mathrm{FDP}}^{\mathrm{B}}=\frac{1}{r}\times {\displaystyle \sum_{i\in \mathrm{S}} Bernoulli\left({\mathrm{fdr}}_i^{\mathrm{B}}\right)}\). Note that of the r significant genes, the q
B is the average bootstrapped false discovery probability, and the FDPB, the bootstrapped proportion of false positives.
A total of 10,000 bootstrapped samples were generated to estimate the bootstrapped standard errors for the local FDRs, q-value and FDP, respectively. For independent genes, the 95 % bootstrapped percentile confidence intervals for local FDR and q-value at various P-value cutoffs can maintain the coverage probabilities close to the nominal value of 0.95, but for correlated genes, the coverage is below 0.95 (Additional file 1). In practice, it is difficult to tell whether the genes under study are independent of one another or are correlated. Therefore, the bootstrapped standard errors presented in this paper should better be regarded as lower bounds of the variability of the FDR control.
