A simple approach to integrate high-throughput functional datasets (e.g., from studies of the transcriptome, proteome, or microRNAome) with genome-wide genotype data obtained from microarray- or sequencing-based studies is to select SNPs that meet certain functional criteria as illustrated in the example in Figure 1. In the first step of this example, SNPs are filtered by requiring that they be associated with genes whose expression levels are associated with the phenotype (Zhong et al., 2010). In the next step, we further reduce the number of SNPs by requiring that they be associated with protein levels that are themselves associated with the phenotype. This process can continue using other omics datasets.

To simplify the description, we focus on the case in which only the gene (mRNA) expression data are integrated, which is depicted with the diagram in Figure 2. This triangle approach and variations thereof were proposed by Huang et al. (2007b) and others (Zhong et al., 2010) and applied to an array of cellular phenotypes. The first step of this method aims to identify a set of gene expression traits associated with the given phenotype at an arbitrarily set p-value threshold, p < p_{gene-phenotype}. It is important to emphasize that this threshold, as in the subsequent thresholds to be defined below, is generally set arbitrarily. In practice, these thresholds are used to prioritize genes or SNPs for downstream analyses. Indeed, one aim of our study is to quantify the significance of an association from a triangle method regardless of the choice of thresholds used during the integrative process. The second step of the method is to identify SNPs that are associated with the selected gene expression traits again at an arbitrarily set threshold, p < p_SNP-gene. At a stringent threshold, this step maps the gene expression traits to genomic loci; this step thus identifies the eQTLs for the corresponding genes. Finally, in the last step of the triangle, the resulting SNPs are interrogated for association with the phenotype. Our primary aim is to describe a method to quantify the significance of the SNPs resulting from this multi-step “triangle” approach.

Since the triangle method is a multi-step approach that derives a final SNP set from a series of (potentially) increasingly stringent thresholds, it is reasonable to expect that such an approach should yield a final set with substantially reduced false discovery rates (FDRs) for association with the phenotype. A simple approach to assess the significance of the findings for this subset of SNPs would be to compute the FDR for them (Storey and Tibshirani, 2003). We illustrate the problem of this approach in Figure 3 in which we show the QQ plot of the associations after applying the triangle method to a simulated phenotype, which has no association with genotype. In this particular example, the first threshold p_{gene-phenotype} was set at 0.05 while p_SNP-gene was set at 5 × 10⁻⁶. Circles above the red line represent SNPs with FDR < 0.05. (Strictly speaking, circles with p-values less than the one with the largest p-value that goes above the red line has FDR < 0.05.) As the figure indicates, the triangle method may yield several spurious associations, if we rely on a “naïve” FDR approach. This example shows the need to develop a more sophisticated approach to estimate the significance of results in this integrative context.

We describe here our approach to generating an empirical null distribution of p-values (Figure 4). First, let Y₁, Y₂, Y₃, …, Y_n be simulated phenotypes obtained from permuting the phenotype data. (Typically, n = 1000.) In case covariates are used, they should be relabeled in sync with the phenotype. For each simulated phenotype, we apply the same triangle method. For each Y_i, we derive the set of gene expression traits g_ij that meet the threshold, p-value < p_{gene-phenotype}, where the associations between the phenotype Y_i and gene expression traits are calculated while preserving the correlation structure of all gene expression phenotypes. For each g_ij, we retrieve the set of eQTLs, S_ijk, associated with the gene at the pre-defined threshold, p-value < p_SNP-gene. The subset of these eQTL SNPs that satisfy p-value < p_{SNP-phenotype} provides a set of p-values {Pijk′}, for each simulated phenotype Y_i. Note that each such set {Pijk′} of p-values may differ in count between simulated phenotypes. Note that i indexes simulations, j indexes genes, and k indexes eQTLs.

We utilize these sets of p-values derived from simulated phenotypes to estimate the null distribution of p-values. Having shown the limitation of the use of the traditional FDR for the integrative triangle method, we derive a simple formula to estimate the FDR using this empirical null distribution.

We closely follow Storey’s approach (Storey and Tibshirani, 2003) to estimate the FDR. The difference in our approach is that we do not assume that the null distribution of p-values is uniform. Instead, we use the empirical distribution generated by simulating the phenotype and performing the integrative analysis. We define the significance level t and reject the null hypothesis of no association for all p-values smaller than t. We use the actual values in the observed vector of p-values as cutoff. Thus, for each p-value, t, in the observed vector of p-values, we compute the FDR of the strategy of rejecting all p-values less than or equal to t. Let the number of falsely significant SNPs be denoted as F(t) = #{nullp_i ≤ t, i = 1, …, m} and the number of significant SNPs be denoted as S(t) = #{p_i ≤ t, i = 1, …, m} with m the total number of SNPs after applying the integrative approach. We estimate the FDR as follows:

where E[.] is the expectation operator. The approximate equality in Eq. 1 is proven by Storey (2003).

The denominator is estimated using the observed number of significant SNPs p ≤ t,

The numerator can be factored as P(p ≤ t and null) = P(p ≤ t | null)·P(null). The first factor P(p ≤ t | null) is estimated using the empirical distribution: #{p_sim,i ≤ t, i = 1, …, M₀}/M₀ where the p_sim’s are the p-values generated with the simulated phenotypes and M₀ is the sum (across all simulations, M₀ = Σm_o,s, where m_o,s corresponds to the number of eQTLs selected after applying the triangle method to the simulated phenotype Y_s) of the total number of SNPs selected using the simulated phenotypes. Note that for uniformly distributed p-values, we would have P(p ≤ t | null) = t. We know, however, that when the set of SNPs are derived from the integrative approach, the null p-values may not be distributed uniformly, as illustrated in Figure 5. The second factor P(null) is the proportion of SNPs that are unrelated to the phenotype and may be estimated as the ratio

or may be set to 1 to yield a more conservative estimate of FDR.

In summary, we estimate the FDR based on the empirical distribution as follows:

with λ = 0.5. We can also use the more conservative estimate

In order to ensure increasing FDR for increasing p-values, we define q-values as

A triangle method, similar to the one described here, had been originally applied to cellular sensitivity data for etoposide (Huang et al., 2007b), one of the most widely used anti-cancer agents. Using our empirical FDR approach, we re-analyzed the same phenotype data from the original experiments, which had sought to quantify the cytotoxic effect of the drug on the cell lines using a colorimetric-based assay, as previously described (Huang et al., 2007b). We conducted our study on the 90 HapMap cell lines of European descent (CEU). The quantitative trait used here was IC₅₀, defined as the concentration required for 50% cell growth inhibition.

We provide an R function for estimating the empirical FDR that can be used once the observed and the simulated p-values are generated (http://www.scandb.org/newinterface/empiricalFDR.R). The way these p-values are generated will depend on the specific integration method used, the eQTL mapping database, and the number of components in the “genetic machinery.”

The GWAS of etoposide IC₅₀ did not yield any significant signals, as perhaps expected from the small sample size. Figure 5 shows a QQ plot of the distribution of p-values (as circles). However, we found a highly significant enrichment for gene regulatory signals among the etoposide-associated variants relative to frequency-matched SNPs (Gamazon et al., 2010a). This finding raises the possibility of the use of eQTL annotation to increase the power to detect true associations. We therefore proceeded to incorporate eQTL functional annotation through the integrative triangle method.

Expression levels had been generated by our group on these cell lines for more than 10,000 genes, allowing us to perform associations between etoposide IC₅₀ and gene expression traits; those genes meeting p < 0.05 (see Table S1 in Supplementary Material) were carried forward in the triangular analysis. We then utilized SCAN (Gamazon et al., 2010b), a public repository for the results of our eQTL studies on the HapMap cell lines, to annotate the selected genes showing association with etoposide IC₅₀ with expression-associated SNPs (p < 10⁻⁴; see Table S2 in Supplementary Material). Finally, the selected eQTLs were tested for association with etoposide IC₅₀. Figure 5 shows the QQ plot of association p-values for all SNPs, the QQ plot for the final SNP set derived from the triangle method, and the QQ plot for the triangle method-prioritized SNPs from each of 1000 simulated phenotypes. The figure illustrates that certain eQTL SNPs from this triangle method-derived SNP set attained a (traditional) FDR < 0.05, but also that the triangle method may yield spurious associations using the traditional FDR.

We applied our proposed empirical FDR method to the observed set of p-values from the triangle analysis-derived set of SNPs. To generate an empirical null distribution of p-values, we conducted simulations (see Materials and Methods). Table S3 in Supplementary Material lists the most significant etoposide-associated SNPs based on our empirical FDR method. Note the comparison between traditional FDR and eFDR for the most highly ranked SNPs prioritized by the triangle method (based on unadjusted p-value), showing that traditional FDR inflates the significance of selected SNPs.