ConReg-R: Extrapolative recalibration of the empirical distribution of p-values to improve false discovery rate estimates

Background False discovery rate (FDR) control is commonly accepted as the most appropriate error control in multiple hypothesis testing problems. The accuracy of FDR estimation depends on the accuracy of the estimation of p-values from each test and validity of the underlying assumptions of the distribution. However, in many practical testing problems such as in genomics, the p-values could be under-estimated or over-estimated for many known or unknown reasons. Consequently, FDR estimation would then be influenced and lose its veracity. Results We propose a new extrapolative method called Constrained Regression Recalibration (ConReg-R) to recalibrate the empirical p-values by modeling their distribution to improve the FDR estimates. Our ConReg-R method is based on the observation that accurately estimated p-values from true null hypotheses follow uniform distribution and the observed distribution of p-values is indeed a mixture of distributions of p-values from true null hypotheses and true alternative hypotheses. Hence, ConReg-R recalibrates the observed p-values so that they exhibit the properties of an ideal empirical p-value distribution. The proportion of true null hypotheses (π0) and FDR are estimated after the recalibration. Conclusions ConReg-R provides an efficient way to improve the FDR estimates. It only requires the p-values from the tests and avoids permutation of the original test data. We demonstrate that the proposed method significantly improves FDR estimation on several gene expression datasets obtained from microarray and RNA-seq experiments. Reviewers The manuscript was reviewed by Prof. Vladimir Kuznetsov, Prof. Philippe Broet, and Prof. Hongfang Liu (nominated by Prof. Yuriy Gusev).


Background
In high-throughput biological data analysis, multiple hypothesis testing is employed to address certain biological problems. Appropriate tests are chosen for the data, and the p-values are then computed under some distributional assumptions. Due to the large number of tests performed, error rate controls (which focus on the occurrence of false positives) are commonly used to measure the statistical significance. False discovery rate (FDR) control is accepted as the most appropriate error control. Other useful error rate controls include conditional FDR (cFDR) [1], positive FDR (pFDR) [2] and local FDR (lFDR) [3] which have similar interpretations as that of FDR. However, appropriate FDR estimation depends on the precise p-values from each test and the validity of the underlying assumptions of the distribution.
The p-values from multiple hypothesis testing, for n hypotheses, can be described by a mixture model g(p) (1) with two components: one component g 0 (p) originates from true null hypotheses and follows uniform distribution U(0, 1) [4], and the other component g 1 (p) results from true alternative hypotheses and follows a distribution confined to the p-values close to 0 [5]. The mixing parameter, π 0 , is the proportion of true null hypotheses in the data. More precisely, g(p) = π 0 g 0 (p) + (1 − π 0 )g 1 (p) where g 0 (p) = 1 denotes the probability density function of a uniform distribution over (0, 1) and g 1 (p) ≈ 0 for p close to 1. Therefore, g(p) will be close to a constant (i.e. π 0 ) for p close to 1. FDR in multiple hypothesis testing for a given p-value threshold a is estimated as π 0 can be estimated from this mixture model in equation (1) as [2] where b is typically chosen to be 0.25, 0.5 or 0.75. These estimates are reasonable under the uniform distribution assumption of a component in this mixture model [6].
However, in many applied testing problems, the pvalues could be under-estimated or over-estimated for many known or unknown reasons. The violation of pvalue distribution assumptions may lead to inaccurate FDR estimation. There are many factors influencing FDR estimation in the analysis of high-throughput biological data such as microarray and sequencing studies. Dependence among the test statistics is one of the major factors [7,8]. Usually in microarray data, there are many groups of genes having similar expression patterns and the test statistics (for example, t-statistic) are not independent within one group. The global effects in the array may also influence the dependence in the data. For example, batch and cluster effects [9,10] always occur in the experiments and sometimes they may be the major cause of incorrectly estimated FDR.
Further, due to the "large p, small n" problem [11] for the gene expression data, some parameters such as mean and variance for each gene cannot be well estimated, or the test assumptions are not satisfied or the distribution of the statistic under null hypotheses may not be accurate. Therefore, many applied testing methods modified the standard testing methods (for example, modifying t-statistic to moderated t-statistic [12]) to increase their usability. As the modified test statistics only approximately follow some known distribution, the approximate p-value estimation may influence the FDR estimation. Resampling strategies may better estimate the underlying distributions of the test statistics. However, due to small sample size and data correlation, the limited number of permutations and resampling bias [13] also influence the FDR estimation.
To address the above problems, we propose a novel extrapolative recalibration procedure called Constrained Regression Recalibration (ConReg-R) which models the empirical distribution of p-values in multiple hypothesis testing and recalibrates the imprecise pvalue calculation to better approximated p-values to improve the FDR estimation. Our approach focuses on p-values as the p-values from true null hypotheses are expected to follow the uniform distribution and the interference from the distribution of p-values from alternative hypotheses is expected to be minimal towards p = 1. In contrast, the estimation of the empirical null distributions of test statistics may not be accurate as their parametric form may not be known beforehand and their accuracy may depend on the data and the resampling strategy used. ConReg-R first maps the observed p-values to predefined uniformly distributed p-values preserving their rank order and estimates the recalibration mapping function by performing constrained polynomial regression to the k highest p-values. The constrained polynomial regression is implemented by quadratic programming solvers. Finally, the p-values will be recalibrated using the normalized recalibration function. FDR is estimated using the recalibrated p-values and theπ 0 can be determined during ConReg-R procedure. We demonstrate that our ConReg-R procedure can significantly improve the estimation of FDR on simulated data, and also the environmental stress response time course microarray datasets in yeast and a human RNA-seq dataset.

Methods
Under the null hypotheses, the p-values are uniformly distributed. Hence, ConReg-R first generates the uniformly distributed p-values within [0, 1] range.

Uniformly distributed p-value generation
Let p i denotes the p-value of the i th test (i = 1, ..., n), without loss of generality, we assume p 1 ≥ p 2 ≥ ... ≥ p n . If we choose a suitable k <n such that the i th null hypothesis H (i) 0 (i ≤ k) is most likely true, then p 1 , ..., p k correspond to the order statistics of k independent uniformly distributed random variables provided p i 's i(i = 1, ..., k) are correctly estimated.
Let p i denote the ideal p-values under H (i) 0 (i ≤ k) , and suppose p k is known. p i (i ≤ k) can be defined as Therefore, p i (i ≤ k) are uniformly distributed over (p k , 1). Then Using (3), (2) becomes Since k is usually large, k/(k -1) is almost 1, therefore p i in (4) can be approximated as We can estimate the recalibration function f(·), to be described below, between {p i } k i=1 and {p i } k i=1 and apply it to all input p-values to output the recalibrated p-values, p cal i (i = 1, . . . , n) i.e.
By Stone-Weierstrass theorem [14], polynomial functions can well approximate any continuous function in the interval [0,1]. Therefore we use polynomial regression to estimate the recalibration function f(·) satisfying appropriate boundary and monotone constraints.

Constrained Regression Recalibration (ConReg-R)
Let y i = p i and x i = p i (i = 1 ... k), and the recalibration polynomial function f(·) is defined as follows, The constraints f (0) = 0, f (1) = 1 and f' (x) >0 should be imposed to ensure the orders of the p-values remain the same after the transformation. Furthermore, the constraint for either f" (x) >0 or f" (x) <0 indicates the function f should also be a monotonic convex or monotonic concave function to deal with the situations with under-estimated or over-estimated p-values separately and helps in good extrapolation. The constraints f (0) = 0 and f (1) = 1 can be easily met by scaling and shifting the regression function. Therefore, the regression function only depends on the other two constraints which can be combined into one constraint during the regression procedure.
Quadratic programming (QP) [15] is employed to estimate the regression function as follows: Let y = (y 1 , ..., y k ) T , b = (b 0 , ..., b t ) T and Equation (7) can be rewritten more succinctly as and the constrains for the first and second order derivatives of f (X) will be Ab ≥ b where b = (0, ..., 0) T is a 2l × 1 vector and is a 2l × (t + 1) matrix, where a 1 , ..., a l are l randomly generated numbers following U(0, 1) to guarantee this constraint is valid in (0, 1), and c is chosen to be 0 (or 1) if f is desired to be convex (or concave respectively).
The least squares procedure for (8) will minimize Minimizing (9) under Ab ≥ b is equivalent to minimizing under Ab ≥ b, where Q = X T X and q = -X T y. Therefore, the constrained polynomial regression problem can be reformulated as a quadratic programming problem.

Two further modifications
We use QuadProg package in R to solve the quadratic programming problem [16]. Due to floating point errors [17], Q = X T X tends to be positive semidefinite instead of being positive definite. To get around this, we add a sufficiently small positive value (l = 10 -10 ) to the diagonal of Q to guarantee Q' = Q + lI t+1 is positive definite and Q' replaces Q in (10).
Furthermore, the polynomial function may not accurately fit the data due to the limitation of the polynomial maximal power (usually set the maximal power t = 10). We can add the fraction of the power (i.e. a non-integer power) to increase the accuracy of the fit. For example,

Computational procedure
For any given k, after applying ConReg-R, the estimation ofπ 0 and its variation (error) are given bŷ and (12) where MAD denotes the median absolute deviation. The final regression function and optimal k(k best ) are determined by examiningπ 0 (k) and eπ 0 (k) over k. Figure  1 illustrates how to choose k best from the functionπ 0 (k). Ideally,π 0 (k) is not expected to change over a range of k (as shown by the blue dashed line in Figure 1) such that p 1 , ..., p k are most likely to be from null hypotheses. If k is too large, p 1 , ..., p k may contain too many p-values from alternate hypotheses andπ 0 (k) may be wrongly estimated to be close to 1, in an extreme case if k is chosen to be n thenπ 0 (k) = 1. However, the extrapolation in recalibration procedure may be unreliable if only a small number of p-values (i.e. small k) are used for the regression andπ 0 (k) may fluctuate near the real π 0 (the red curve in Figure 1). Therefore, we aim to choose optimal k(k best ) as a trade-off to include just enough p-values from null hypotheses for the regression to achieve good extrapolation. The k that gives stable estimate (eπ 0 (k) < δ) and the last minimum ofπ 0 (k) is chosen to be the k best . The regression function, extrapolation andπ 0 (k) corresponding to k = k best are chosen for recalibrating p-values and re-estimating FDR.
The following is the computational procedure for a given {p i } n i=1 in descending order: 3. Use quadratic programming to obtain regression function h k , where c can be predefined or estimated by checking whether more than half of points for (p i , p i ) are above the diagonal (line from origin to (1, 1)) (c = 1) or below the diagonal (c = 0).
to satisfy constraints f k (0) = 0, and f k (1) = 1. 5. Repeat steps 2-4 for all k, and compute theπ 0 (k) and eπ 0 (k) for each k. Let k best be the maximal of k which locally minimizesπ 0 under the constraint of small eπ 0 , where the cutoff of eπ 0 and local minimization criteria should be predefined. 6. Choose the final regression function f (.) under k best and output recalibrated p-values. 7. Re-estimate the FDR using recalibrated p-values andπ 0 =π 0 (k best ).
R-code for ConReg-R is attached as Additional file 1.

Dependence simulation
Data dependence is one of the major causes for over-or under-estimated p-values. We simulated an expression data, with dependence, Z = (z ij )(i = 1, ..., n, j = 1, ..., r) with n (n = 10000) genes and r(r = 10) replicates using the formula as follows, To compare the result, we also simulated a data set with no dependence using the same procedure but with the dependence effect d ij = 0. One sample t-test was performed to generate p-values. Figure 2 shows the p-value density histograms for π 0 = 0.7 and π 0 = 0.9. As can be seen in the plots B and D in Figure 2, the p-value histograms from independent data have constant frequency for p ≥ 0.5 and the density near 1 indicates theπ 0 . However, the p-value histograms from dependent data (the plots A and C in Figure 2) do not have such constant frequency and p-value density increases as p-value increases in the neighborhood of 1. The density near 1 exceeds the respective π 0 .
ConReg-R used the above p-values as input and output the recalibrated p-values. The results are summarized in Figure 3. For the independent data sets, the algorithm chose k = 0.71n for π 0 = 0.7 and k = 0.64n for π 0 = 0.9 since it locally minimizedπ 0 under  π 0 (k). The p-values do not significantly change after regression. As such, the regression curves almost overlap with the diagonals, and the input p-value histogram and the output p-value histogram are very similar to each other. The FDR estimation errors (the absolute difference between FDR estimated by p-values and real FDR) also do not significantly change after applying ConReg-R and the estimation of FDR is very close to the real FDR. However, for the dependent data sets, the algorithm chose k = 0.62n for π 0 = 0.7 and k = 0.88n for π 0 = 0.9. The regression curves are all below the diagonals and the output p-value histograms after applying ConReg-R appears more like the ones obtained for the independent data. The accuracy of estimated FDR after applying ConReg-R is substantially improved.
To study more complicated dependency situations, we generated dependent datasets with random dependence effect [8] as follows, where r is the correlation constant (here we set r = 0.5) which determines the correlation coefficient between genes. Here b i denotes the biological effect, and d j denotes the random dependence effect. Set b i = 1, if i ≤ n(1 -π 0 ) and b i = 0 if i > n(1 -π 0 ), and d j~N (0, 1). Let ε ij N(0, 1) be the background noise. The result for π 0 = 0.7 and π 0 = 0.9 are shown in Additional File 2, Figure S1. Similar to the simulations of fixed dependence effect, the estimated FDR after applying ConReg-R is closer to real FDR. The results of our procedure for 100 repeated simulations are summarized in the boxand-whisker plots in Figure 4. As shown in this figure, for the independent data sets, the FDR estimation errors (the mean absolute difference between real FDR and the FDR estimated by p-values using Benjamini-Hochberg method) after applying ConReg-R is slightly higher. However, it is still acceptable since most simulation resulted in errors below 0.05. For the dependent data sets with fixed and random dependence effects, the FDR estimation errors after applying ConReg-R are significantly less than those without applying ConReg-R. The FDR estimation for π 0 = 0.9 is even closer to real FDR after applying ConReg-R compared with the result for π 0 = 0.7 because of more p-values used for regression and less number of p-values for extrapolating in datasets of π 0 = 0.9.

Combined p-values simulation
In many analyses, more than one dataset are involved and a meta-analysis by combining p-values from different studies or datasets is needed to estimate the overall significance for each gene. For example, (i) to find genes which are significant in at least one experiment, minimal p-values will be of interest; (ii) to identify genes which are significant across all the experiments, the maximal p-values will be of interest; and (iii) in order to detect genes which are significant on average, the product of p-values will be appropriate. The distribution of combined p-values will not be uniform even under true null hypotheses [18]. For currently used meta-analysis methods, such as "minimal", "maximal" or "product", we can obtain the transformation functions to recalibrate the combined p-values to satisfy the condition of p-values are uniform distributed under true null hypotheses. However, for other more complicated meta-analysis methods, the transformation function cannot be determined accurately leading to under-or over-estimation of significance, and ConReg-R can provide the polynomial function approximation for the unknown transformation.
The two p-value histograms for each π 0 = 0.7 and π 0 = 0.9, and for each of five different methods are plotted in Figure 5. It can be seen from Figure 5 that the p-value histograms after theoretical transformation have constant frequency after 0.5 and the p-value density near 1 indicates theπ 0 . However, the p-value histograms from "Min", "Square", "Prod" shifted towards 0 and the p-value histograms from "Max", "Sqroot" shifted towards 1.
ConReg-R used the above combined p-values as input and the results are shown in Additional File 2, Figure S2 (π 0 = 0.7) and Additional File 1, Figure S3 (π 0 = 0.9). From Figure S2 and Figure S3, the regression curves are monotonic concave functions for "Min", "Square", "Prod" and monotonic convex functions for "Max", "Sqroot". The histograms after applying ConReg-R are also very similar to the theoretical transformed p-value histograms. The FDR estimation improved significantly after applying ConReg-R. It shows that the estimated FDR after applying ConReg-R is more likely to be the real FDR. The results of using our procedure for 100 repeated simulations are summarized in Figure 6. The FDR estimation errors after applying ConReg-R are significantly less than those obtained without applying ConReg-R.

Yeast Environmental Response Data
Yeast environmental stress response data generated by [20,21] for nearly 6000 genes of yeast (S. cerevisiae) was aimed at understanding how yeast adopts or reacts to various stresses present in its environment. We selected 10 datasets: (1) Heat shock from 25°C to 37°C response; (2) Hydrogen peroxide treatment; (3) Menadione exposure; (4) DTT exposure response; (5) Diamide treatment response; (6) Hyper-osmotic shock response; (7) Nitrogen source depletion; (8) Diauxic shift study; and, (9-10) two nearly identical experiments on stationary phase. We used Limma (Linear Models for Microarray Data) [12] package in R to compute p-values for responsiveness of genes for each dataset.
The p-value distribution for each dataset is shown in Additional File 2, Figure S4. As can be seen in Figure S4, the majority of the p-value histograms do not have similar frequency after p = 0.5, and the density near p = 1 is less than π 0 = 0.5. This implies that the p-values were under-estimated and the number of significantly responsive genes under these environmental Error π 0 = 0.9(indep.) π 0 = 0.9(dep.fix) π 0 = 0.9(dep.random) Figure 4 Boxplots of FDR estimation errors (the mean difference between real FDR and the FDR estimated by p-values) for 100 simulations of independent and dependent datasets at π 0 = 0.7 and π 0 = 0.9 before (input) and after applying ConReg-R (calibrated). Table 1 Combined p-values methods [18].  The density histograms for "Min", "Max", "Sqroot", "Square" and "Prod" datasets at π 0 = 0.7 and π 0 = 0.9. (Th.) indicates the density histograms for each method after theoretical transformation. The gray horizontal line indicates the π 0 for each plot.
stresses should be less than observed. We applied Con-Reg-R on the p-values of each dataset. Our result shows that the histograms of recalibrated p-values obtained by applying ConReg-R are better than without recalibration, and π 0 estimations are all above 0.5 ( Figure S4). We use a true positive set of 270 genes from [22] to compute true FDR (FDR r ). This is the intersection of core environmental stress response genes obtained by co-regulation study in [21] and the yeast orthologs of S. pombe stress response genes. These 270 genes have been used as the true positive sets in other studies [23,24]. The true FDR is calculated based on this 270 gene list and we calculated the improvement of FDR estimation (FDR im ) for each dataset after applying ConReg-R. The FDR im is defined as followed: where FDR 1 i (respectively, FDR 0 i ) is the estimated FDR by recalibration (respectively, input) p-values for gene i(i = 1 ... n); and FDR r i is the true FDR for gene i. The improvements in FDR estimation for all 10 datasets are shown in Figure 7. After applying ConReg-R, FDR estimation improved by 15% to 25% which means that the FDR estimation will be closer to the real FDR.
We performed the meta-analysis for 10 datasets to detect the core environmental stress response genes using "maximal" method. The combined p-values are computed by the maximal p-values across 10 datasets, and then transferred to meta analysis p-values by transformation function in Table 1. The p-value density histograms for meta-analysis before and after applying ConReg-R are shown in Figure 8. The meta-analysis p-values show better distribution after first applying ConReg-R to each dataset and then perform the meta analysis. And FDR estimation improved by 38.5% after applying ConReg-R.

Significance Analysis of Differential Expression from RNA-seq Data
The next-generation sequencing technologies have been used for gene expression measurement. In [25], the authors compared RNA-seq and Affymetrix microarray experiments and claimed that the sequencing data identified many more differentially expressed genes between human kidney and liver tissue samples than microarray data using the same FDR cutoff. In total, 11,493 significant genes were identified by RNA-seq (3380 more genes than Affymetrix), only 6534 (56.9%) genes were also identified by Affymetrix experiments. By checking the p-value histograms for RNA-seq dataset, we found that majority of p-values are very significant and its frequencies are very non-uniform for p >0.5. However, the p-value histogram for Affymetrix datasets is close to uniform for p >0.5 (Additional File 2, Figure S5).
( Figure 9). The FDR estimation is improved by 20% after applying ConReg-R if we used significant genes identified by affymetrix experiments as the true positive set.

Discussion
ConReg-R focuses on the uniformity of p-values under null hypotheses and uses constrained polynomial regression to recalibrate the empirical p-value distribution to more well-defined p-value distribution. Therefore, the FDR estimation can be improved after the recalibration since the assumption of FDR estimation is that the input p-values should follow such an ideal empirical p-value distribution under null hypothesis. If the input p-values follow the properties of ideal empirical p-values distribution, the regression function tends to be diagonal line (i.e., y = x) and the p-values do not change considerably after recalibration. Though our method is discussed in the context of global FDR control, it is equally applicable to the other FDR like controls such as local FDR. Our method does not provide any new FDR control, but inputs better  Figure 8 The p-value density histograms for meta-analysis ("Max") before and after applying ConReg-R using yeast environmental response datasets. The gray horizontal line indicates π 0 = 0.5 for each plot.

Figure 9
A Venn diagram summarizing the overlap between significantly differentially expressed genes identified by sequencing (left circle) and microarray (right circle) technologies. The numbers in black are the numbers reported in [25]. The numbers in red are the numbers after applying ConReg-R. In ConReg-R, the cutoff of eπ 0 and local minimization criteria can be changed for choosing the suitable k after checking the plots ofπ 0 and eπ 0 . From the combined p-values simulation, the regression function may not fit the data well for the "Square" case. The fractional power, such as 1/2, can be added in the polynomial function to improve the fit.
The assumption that the regression function is convex or concave can be useful to deal with the imprecise p-values whose distribution is biased towards 1 or 0 respectively. These are the most common cases of the p-values being under-estimated or over-estimated. However, in some cases, the p-values can be under-estimated in one range of p-values and over-estimated in the remaining range of p-values. These p-value distributions may have peak or valley in the middle of the p-value range or even have multiple peaks. The regression function will then no longer be convex nor concave. Regression function to handle this situation is currently under study. Our ConReg-R can be generalized to an iterative weighted least squares method (e.g. decreasing weights from 1 to n). The weight program in the current version of the ConReg-R is assigning a weight of 1 for all p-values from 1 to k and a weight of 0 for the rest. Furthermore, different optimization schemes also need to be experimented. These will be explored in our future work. The distribution of p-values from multiple testing can be modeled by the mixture of uniform distribution and some other well-defined distribution such as Beta distribution [5]. The parametric recalibration method is under development. The discrete p-values from some nonparametric tests cannot be modeled by mixture model and new procedure should be explored to resolve this kind of problem.

Review
Major: Summary. Mathematical part of the manuscript is major, but it is not convicted. Choosing the model (2) is not justified. The used model assumes a uniform distribution of the tail of FDR distribution is not correct due to experimental data; the most empirical FDR distribution exhibits the non-homogeneous fat tail as it was showed in Suppl. 1.
Authors' response: Thanks for your comment. The p-value distribution under true null hypotheses follows the uniform distribution if the data satisfy all the assumptions of the hypothesis testing method [4]. However, the FDR distribution may not have similar property. The model (2) is based on the properties of the p-value distribution, but not on the FDR distribution. For experimental data, many factors influence those assumptions and the p-value distribution may not have uniformly distributed tail which is central to the paper. We clarified this in paragraphs 2, 3 & 4 in page 2.
The final FDR distribution of combined p-value function is not be uniform "even under true null hypotheses", because it is a mixture distribution of the samples from different populations. If the samples are taken from significantly different distributions with different sample size and possible batch effect (not removed), does not allowed you to combine the datasets. If the samples are taken from the same distribution, why should not the p-value distributions be uniform (assuming that they were uniform in the individual datasets before merging)?
Authors' response: In combined p-values simulation, we assume p-values follow the uniform distribution for each dataset under true null hypotheses. For each gene, the overall p-value was obtained by combining its p-values across experiments using different combination methods which are typical non-parametric meta-analysis methods. The combined p-value distributions are not uniform which explains why another distribution is used to derive the overall p-value [18]. We clarified this in paragraphs 1 and 2 in page 9.
Questionable statements p.3 "Dependence among the test statistics is one of the major factors". It would be nice if you give examples of the test statistics you are mentioning.
Authors' response: We added the "t-statistic" as an example in text (paragraph 1 in page 3). "large p, small n problem". Could you describe the meaning this problem?
Authors' response: "large p, small n" problem is the problem that the number of variables (p) is much bigger than sample size (n) which always occurs in gene expression studies. We added one reference [11] for this problem (paragraph 2 in page 3). The n and p used in this context are different from that used in the discussion throughout the paper.
"ConReg-R first maps the observed p-values to uniformly distributed p-values". From the context of this paragraph it is unclear where the set of 'uniformly distributed p-values' is obtained from.
Authors' response: We changed to "predefined uniformly distributed p-values" (paragraph 3 in page 3).
p.4 "Let pi' denote the ideal p-values". How did you decide that the ideal distribution of the p-values is a linear function, as defined in eq2 and is data-and test-independent?