TreeQTL: hierarchical error control for eQTL findings

TL;DR

TreeQTL improves error control in eQTL studies by using a hierarchical testing approach.

Contribution

TreeQTL introduces a novel hierarchical error control method for eQTL analysis.

Findings

TreeQTL provides better error rate control compared to traditional methods.

The hierarchical approach improves the reliability of eQTL findings.

Abstract

Commonly used multiplicity adjustments fail to control the error rate for reported findings in many expression quantitative trait loci (eQTL) studies. TreeQTL implements a stage-wise multiple testing procedure which allows control of appropriate error rates defined relative to a hierarchical grouping of the eQTL hypotheses.

The R package TreeQTL is available for download at http://bioinformatics.org/treeqtl.

sabatti@stanford.edu

Figures

Organization of eQTL hypotheses in TreeQTL. Hypotheses testing local regulation are on the left and hypotheses testing distal regulation are on the right. Local regulation hypotheses have been grouped by probe and distal regulation hypotheses are grouped by SNP, so that Level 1 rejections will result in discoveries of LeProbes and DeSNPs. Tested hypotheses are colored in red, and rejected hypotheses indicated with a star.

Full text

Introduction

The overarching goal of eQTL analysis is to gain insight into the genetic regulation of gene expression. Typically, this is carried out by testing a vast collection of hypotheses Hvp probing association (or linkage) between the genotype at variant (locus) v and the measured expression for probe p, where v = 1, …, M, p = 1, …, G, M is on the order of hundreds of thousands, and G of tens of thousands. Given the large number of hypotheses tested, the need to adjust for multiplicity is universally recognized and the false discovery rate (FDR) [5] is typically adopted as the target global error rate.

In an effort to improve interpretability, reporting of results is typically organized along more general findings such as the discovery of genes subject to local regulation [7] or regulatory SNPs (eSNPs) [9]. The adopted strategy for multiplicity adjustment needs to offer guarantees on the specific reported discoveries. For example, imagine testing the Hvp hypotheses using the Benjamini Hochberg (BH) [5] rule and defining as an eSNP those variants v for which Hvp is rejected for at least some probe p. While this would control the FDR among the Hvp rejections, it would offer no control of the rate of false eSNP discoveries, as shown by the simulations in [10].

Researchers in the eQTL field have recognized this challenge and have noted that since local regulation is more common than distal [1], hypotheses probing these two mechanisms should be tested separately. However, there is no single standard in the literature for error rates targeted or error-controlling strategies: for example, one finds the notion of per-gene error rates [9, 13] or the application of Bonferroni across genes [14] in local regulation, while for distal effects significance cut-offs vary from 5 × 10^−8^ [8] to 5.78 × 10^−12^ [14]. This makes comparison across studies and replicability challenging.

Approach

1

To overcome the confusion generated by the plurality of approaches and to provide guarantees relative to the discoveries reported, it is useful to recognize the structure among the hypotheses tested in an eQTL study and to introduce some terminology. We distinguish between hypotheses testing local (when the distance between variant v and probe p is less then a threshold) and distal regulation, indicating them with Lvp and Dvp, respectively.

Further, we recognize that we might be interested mainly in identifying which genes appear to have local (distal) regulation prior to obtaining a detailed list of the variants involved in this regulation; or we might want to pinpoint SNPs that appear to have local (distal) effects on multiple genes, even without committing to a comprehensive list of these genes. Recalling that our measurements of expression are mediated by probes and that typically typed genetic variants are SNPs, we use LeProbe to signify a probe whose expression is influenced by some local DNA variants and LeSNP to designate a specific SNP associated to variability in expression for some local probes. With LeAssociation we signify the local association between one specific SNP and one specific probe. For distal regulation, DeProbe, DeSNP, and DeAssociation are similarly defined. Mapping these back to the original collection of hypotheses, we note, for example, that to discover a DeSNP is equivalent to rejecting the intersection hypothesis Dv•=∩p D_vp_.

To control global error rates defined in terms of the reported discoveries (LeSNPs, LeAssociations, etc.), we have implemented in TreeQTL a multi-resolution approach based on results in [3] and whose practical effectiveness in GWAS has been described in [10]. Furthermore, by accounting for the described structure, TreeQTL has the potential to increase power: the effective number of tests is reduced and one can capitalize on the adaptivity of FDR.

Methods

2

TreeQTL is a hierarchical testing procedure that distinguishes two levels of discoveries within each class of local or distal regulation. In Level 1, the users can specify their primary interest as the identification of either eProbes or eSNPs. Given this choice, all the pair-wise association hypotheses are given a position in a tree similar to that shown in Figure 1: each eProbe or eSNP hypothesis indexes the collection of simple hypotheses by whose intersection it is defined. Level 1 hypotheses are tested controlling for global errors within the two regulation classes. The more granular eAssociation hypotheses in Level 2 are tested only when they correspond to a global hypothesis rejected in Level 1.

Organization of eQTL hypotheses in TreeQTL. Hypotheses testing local regulation are on the left and hypotheses testing distal regulation are on the right. Local regulation hypotheses have been grouped by probe and distal regulation hypotheses are grouped by SNP, so that Level 1 rejections will result in discoveries of LeProbes and DeSNPs. Tested hypotheses are colored in red, and rejected hypotheses indicated with a star.

TreeQTL takes as input the p-values for each hypothesis in Level 2: these may be obtained via Matrix eQTL [11] and their correctness is of crucial importance. The p-values for the Level 1 hypotheses are computed using Simes’ rule [12] on the families they index. This summary of the evidence for the global null hypotheses is relatively robust to dependence [4]. Users can, however, input alternative p-values for Level 1, such as those obtained via permutation.

Testing proceeds from Level 1, where the default option is to control the FDR at a user-specified level q1 using the BH [5]. This guarantees control of the FDR for Level 1 discoveries when all hypotheses are independent or under positive regression dependence. The Benjamini-Yekutieli procedure [6], robust to any dependence between tests, is available as an alternative. A third, more stringent option is to control the family-wise error rate (FWER) using the Bonferroni correction.

In Level 2, significance is established using BH within each set of eAssociation hypotheses corresponding to an eSNP or eProbe identified at Level 1, at the more stringent level needed to account for selection [3]. The expected average proportion of false discoveries across the selected Level 2 families is then controlled to the user-specified target level q2.

Example application

3

To demonstrate feasibility, we applied TreeQTL to whole-blood data from the pilot phase of the GTEx project [2]: genotype data at 6,820,472 SNPs and expression levels for 30,115 genes are available across 156 subjects. Following the steps in [2], p-values were obtained by applying Matrix eQTL [11] to quantile-normalized gene expression, adjusting for both known and unknown technical covariates (note that power for distal analysis is likely to be limited with this sample size). In applying TreeQTL, we identified eSNPs as the discovery of interest in Level 1, and set q1 = q2 = 0.05 for local and q1 = q2 = 0.01 for distal regulation. This led to the discovery of 199,032 LeSNPs, involved in 325,529 LeAssociations, and 164,860 DeSNPs involved in 216,933 DeAssociations. By comparison applying BH at levels q = 0.05 and q = 0.01 across all local and distal hypotheses, respectively, led to the identification of 250,945 LeSNPs (with 389,507 LeAssociations) and 179,625 DeSNPs (with 216,683 DeAssociations). TreeQTL, then, results in a similar number of eAssociations, but these are grouped under a small number of eSNPs.

Conclusion

4

By analyzing local and distal regulation separately, TreeQTL has less stringent cut-offs for tests probing local effects, where it is expected that a larger number of hypotheses will be non-null. By grouping hypotheses relative to the same SNP or the same probe, TreeQTL capitalizes on the inherent heterogeneity of the problem. For example, while few SNPs will act as distal regulators of expression and might influence multiple probes, the vast majority of SNPs will not have such an effect: testing all the eAssociation hypotheses relative to one SNP together, separately from those concerning other loci, allows greater power for SNPs with true regulatory roles and reduces false positives. Finally, the hierarchical structure of TreeQTL assures control of the FDR for eSNP and eProbe discoveries. While the current version of TreeQTL implements methodology relative to studies involving only one tissue, future releases will incorporate approaches for the more complex structure of multi-tissue investigations.