Statistical significance for genomewide studies

Abstract

Some of the earliest genomewide studies involved testing for linkage at loci spanning a large portion of the genome. Because a separate statistical test is performed at each locus, traditional p-value cutoffs of 0.01 or 0.05 had to be made stricter to avoid an abundance of false positive results. The threshold for significance in linkage analysis is usually chosen so that the probability of any single false positive among all loci tested is ≤0.05. This strict criterion is used mainly because one or very few loci are expected to show linkage in any given study (1, 2). Because of the recent surge in high-throughput technologies and genome projects, many more types of genomewide studies are now underway. The analyses of these data also involve performing statistical tests on thousands of features in a genome. As opposed to the linkage case, it is expected that many more than one or two of the tested features are statistically significant. Guarding against any single false positive occurring is often going to be much too strict and will lead to many missed findings. The goal is therefore to identify as many significant features in the genome as possible, while incurring a relatively low proportion of false positives.

We are specifically concerned with situations in which a well defined statistical hypothesis test is performed on each of thousands of features represented in a genome. These “features” can be genes, all nucleotide words of a certain length, single-nucleotide polymorphism markers, etc. Several motivating examples are given below. For each feature, a null hypothesis is tested against an alternative hypothesis. In this work, we say that a feature is truly null if the null hypothesis is true, and a feature is truly alternative if the alternative hypothesis is true. If a feature is called significant, then the null hypothesis is rejected in favor of the alternative hypothesis. The goal is to propose and estimate a measure of significance for each feature that meets the practical goals of the genomewide study and that is easily interpreted in terms of the simultaneous testing of thousands of features.

We propose that the recently introduced q value (3, 4) is a well suited measure of significance for this growing class of genomewide tests of significance. The q value is an extension of a quantity called the “false discovery rate” (FDR) (5), which has received much recent attention in the statistics literature (6–11). A FDR method has been used in detecting differential gene expression in DNA microarray experiments (12), which can be shown to be equivalent to the method in ref. 5 under certain assumptions. Also, ideas similar to FDRs have appeared in the genetics literature (1, 13).

Similarly to the p value, the q value gives each feature its own individual measure of significance. Whereas the p value is a measure of significance in terms of the false positive rate, the q value is a measure in terms of the FDR. The false positive rate and FDR are often mistakenly equated, but their difference is actually very important. Given a rule for calling features significant, the false positive rate is the rate that truly null features are called significant. The FDR is the rate that significant features are truly null. For example, a false positive rate of 5% means that on average 5% of the truly null features in the study will be called significant. A FDR of 5% means that among all features called significant, 5% of these are truly null on average.

The q value provides a measure of each feature's significance, automatically taking into account the fact that thousands are simultaneously being tested. Suppose that features with q values ≤5% are called significant in some genomewide test of significance. This results in a FDR of 5% among the significant features. A p-value threshold of 5% yields a false positive rate of 5% among all null features in the data set. In light of the definition of the false positive rate, a p-value cutoff says little about the content of the features actually called significant. The q values directly provide a meaningful measure among the features called significant. Because significant features will likely undergo some subsequent biological verification, a q-value threshold can be phrased in practical terms as the proportion of significant features that turn out to be false leads.

Here we show that the FDR is a sensible measure of the balance between the number of true positives and false positives in many genomewide studies. We motivate our proposed approach in the context of several recent and prominent papers in which awkwardly chosen p-value cutoffs were used in an attempt to achieve at least qualitatively what the q value directly achieves. We also introduce a fully automated method for estimating q values, with an initial treatment of dependence issues between the features and guidelines as to when the estimates are accurate. The proposed methodology is applied to some gene expression data taken from cancer tumors (14), supporting previously shown results and providing some additional information.

Appendix

Remark A: FDR, Positive FDR (pFDR), and the q Value. In this article, we have used FDR and FDR = E[F/S] somewhat loosely. It will almost always be the case that S = 0 with positive probability, which implies that E[F/S] is undefined. The quantity E[F/S|S > 0]·Pr(S > 0) was proposed as a solution to this problem (5), which is the result of setting F/S = 0 whenever S = 0 in the original E[F/S]. This quantity is technically called the FDR in the statistics literature. In our case we want to place a measure of significance on each feature, which is done under the assumption that the feature is called significant. Thus, the inclusion of Pr(S > 0) is somewhat awkward. An alternative quantity, called the pFDR, was recently proposed (23), which is simply defined as pFDR = E[F/S|S > 0]. The q value is most technically defined as the minimum pFDR at which the feature can be called significant (24). Because m is large in genomewide studies, we have that Pr(S > 0) ≈ 1 and FDR ≈ pFDR ≈ E[F]/E[S], so the distinction is not crucial here. Also, the estimate we use is easily motivated for either quantity (4, 10).

Suppose that each feature's statistic probabilistically follows a random mixture of a null distribution and an alternative distribution. Then under a fixed significance rule, the pFDR can be written as Pr(feature i is truly null|feature i is significant), for any i = 1,..., m (3). Similarly, the false positive rate can be written as Pr(feature i is significant|feature i is truly null), for any i = 1, 2,..., m. Notice the similarity between the pFDR and false positive rate: the arguments have simply been swapped in the conditional probabilities. This connection is the motivation for calling our proposed quantity q value. Indeed, the p value of a feature is technically defined to be the minimum possible false positive rate when calling that feature significant (26). Likewise, the q value is based on the minimum possible pFDR.

Remark B: General Algorithm for Estimating q Values. There is a tradeoff between bias and variance in choosing the λ to use in . For well formed p values, it should be the case that the bias of decreases with increasing λ, the bias being the smallest when λ → 1 (4). Therefore, the method we use here is to estimate . In doing so, we will borrow strength across the over a range of λ, giving an implicit balance between bias and variance.

Consider Fig. 3, where we have plotted versus λ for λ = 0, 0.01, 0.02,..., 0.95. By fitting a natural cubic spline to these data (solid line), we have estimated the overall trend of as λ increases. We purposely set the degrees of freedom of the natural cubic spline to 3; this means we limit its curvature to be like a quadratic function, which is suitable for our purposes. It can be seen from Fig. 3 that the natural cubic spline fits the points quite well. The natural cubic spline evaluated at λ = 1 is our final estimate of π₀. For a variety of simulations and forms of dependence (data not shown), this method performed well, often eliminating all bias in .

Fig. 3.

The versus λ for the data of Hedenfalk et al. (14). The solid line is a natural cubic spline fit to these points to estimate .

The following is the general algorithm for estimating q values from a list of p values.

1. Let p₍₁₎ ≤ p₍₂₎ ≤... ≤ p_(m) be the ordered p values. This also denotes the ordering of the features in terms of their evidence against the null hypothesis.

2. For a range of λ, say λ = 0, 0.01, 0.02,..., 0.95, calculate

   

3. Let $\widehat{f}$be the natural cubic spline with 3 df of on λ.

4. Set the estimate of π₀ to be

   

5. Calculate

   

6. For i = m – 1, m – 2,..., 1, calculate

   

7. The estimated q value for the ith most significant feature is $\widehat{q}$(p_(i)).

Remark C: Analysis of the Hedenfalk et al. Data. The data from ref. 14 can be obtained at http://research.nhgri.nih.gov/microarray/NEJMSupplement. The data consist of 3,226 genes on n₁ = 7 BRCA1 arrays and n₂ = 8 BRCA2 arrays, along with some arrays from sporadic breast cancer, which we did not use. If any gene had one or more measurement exceeding 20, then this gene was eliminated. A value of 20 is several interquartile ranges away from the interquartile range of all of the data and did not seem trustworthy for this example. This left m = 3,170 genes.

We tested each gene for differential expression between these two tumor types by using a two-sample t statistic. Let the log₂ expression value from the jth array and the ith gene be denoted by x_ij. Then $\bar{x}$_i2 = 1/n₂∑ _j∈BRCA2x_ij and are the sample mean and variance for gene i among the arrays taken from BRCA2 tumors. We can similarly define $\bar{x}$_i1 and to be the sample mean and variance for the ith gene among the BRCA1 tumor arrays. The two-sample t statistic for the ith gene, allowing for the possibility that the tumors have different variances, is then

for i = 1, 2,..., 3,170.

We next calculated null versions of t₁, t₂,..., t₃₁₇₀ when there is no differential gene expression. Because it is not clearly valid to assume that the t_i follow a t distribution, we calculate these by a permutation method. Consider all possible ways to assign n = 15 arrays to n₁ = 7 arrays from BRCA1 and n₂ = 8 arrays from BRCA2. Under the assumption that there is no differential gene expression, the t statistic should have the same distribution regardless of how we make these assignments. Specifically, the labels on the arrays are randomly scrambled, and the t statistics are recomputed. Therefore, for B = 100 permutations of the array labels we get a set of null statistics , b = 1,... B. The p value for gene i, i = 1, 2,..., 3,170 was calculated by

We estimated the q values for differential gene expression between the BRCA1 and BRCA2 tumors by using the algorithm presented above. All results, including the computer code used to analyze the data, can be found at http://genomine.org/qvalue/results.html.

Remark D: Theoretical Properties. Some mathematical results hold under “weak dependence” of the p values (or features in the genome). These mathematical results indicate that our method yields conservative q-value estimates. The conservative property is desirable because one does not want to underestimate the true q values. (For the same reason one would not want to underestimate a p value.)

Suppose that with probability 1, we have S(t)/m → G(t) and F(t)/m₀ → G₀(t) for each t ∈ [0, 1] as m → ∞, where G and G₀ are continuous functions. In words, this says that the empirical distribution functions of the observed p values and null p values converge pointwise to some continuous functions. Weak dependence is defined as dependence that allows this pointwise convergence. (As a rule of thumb, the more local the dependence is, the more likely it is to meet the weak dependence criterion.) Also suppose that G₀(t) ≤ t (i.e., uniform distribution or more conservative), and that m₀/m converges. If we constrain (which should usually be the case), then it can be shown that for any δ > 0,

which means that the estimated q values are simultaneously conservative for the true q values, even when taking the worst-case scenario over [δ, 1] for arbitrarily small δ. Also, we can conclude that

which means that if we call all genes with q values ≤ α, then in the long run the FDR will be ≤ α. The proofs of these claims follow from minor modifications to some of the main results in ref. 10.

Notes

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: FDR, false discovery rate; pFDR, positive FDR.