Significant independent associations detected for 128 UK Biobank traits

To compare the power of gene-based tests in empirical data, we evaluated the numbers of significant independent gene-based associations detected for each method across 128 approximately independent GWAS traits in the UK Biobank (selection procedures are described in Materials and methods). The number of independent associations is calculated for each trait by selecting the most significant gene-based association p-value, masking all gene-based tests that include variants within 1 Mbp of variants for the selected gene, and repeating until all genes with Bonferroni-adjusted p-value ≤ 5% are either selected or masked. This procedure ensures that all selected genes are separated by at least 1 Mbp, and provides a conservative estimate of the number of significant independent signals. The omnibus test (“GAMBIT”) detected significantly more associations than other gene-based association methods overall (Fig 5A), and consistently detected more associations than other methods across a wide range of traits (Fig 5B). We also compared the numbers of significant associations for each method without filtering or LD pruning (Fig 6). Statistics that incorporate many variants over a broad region for each gene (e.g., dTSS-weighted tests) yield substantially more significant associations, as expected.

Fig 5

UK Biobank analysis: Numbers of significant independent associations detected.

Numbers of independent gene-based associations (at Bonferroni-corrected 5% significance level) detected by each method across 128 UK Biobank traits. Panel A: Total number of significant independent associations across traits (delineated by horizontal black lines) for each gene-based test; Wilcoxon signed-rank p-values (top) for paired comparisons between no. associations detected by the omnibus test (“GAMBIT”; red) versus Pascal/SOCS (blue) and single-annotation gene-based tests (green). The omnibus test detects significantly more associations than any individual constituent gene-based test or by Pascal/SOCS across UK Biobank traits. Panel B: Comparison of total numbers of genes detected across individual traits for the omnibus test (y-axis) versus single-annotation tests (x-axis).

Fig 6

UK Biobank analysis: Overlap between gene-based association methods.

Panel A: Total number of significant genes (p-value < 2.5e-6) for each method across all 128 traits. Unlike Fig 5, gene-based associations in Fig 6 are not filtered or LD pruned, and a single significant GWAS variant can produce multiple significant gene-based associations for a given method. Here, a larger number of significant genes does not necessarily suggest greater statistical power. Panel B: The i, j^th heatmap element can be interpreted as the conditional probability that gene-based test i is significant given that gene-based test j is significant, which is estimated as the total number of overlapping significant genes between tests i and j divided by the total number of significant genes for test j.

Concordance with benchmark genes for 25 UK Biobank traits

We compiled lists of benchmark genes from the ClinVar database [21] and the Human Phenotype Ontology (HPO) [22] for 25 traits in the UK Biobank to compare the gene-based analysis methods identifying causal genes; procedures and selection criteria are detailed in Materials and Methods. Results are shown separately using the union and intersection of ClinVar and HPO benchmark genes; the latter gene set is expected to have higher specificity, albeit fewer genes. Performance identifying benchmark genes was assessed by ranking genes separately within each benchmark locus for each UK Biobank trait, where a benchmark locus is defined as the set of all genes within 1 Mbp of a genome-wide significant single-variant association that also is within 1 Mbp of a benchmark gene. To compare the performance of gene ranking methods, we calculated the fraction of loci at which the top-ranked gene coincides with a benchmark gene (Fig 7) and assessed receiver operating characteristic (ROC) and precision-recall curves for each method (S2 Fig).

Fig 7

UK Biobank analysis: Performance identifying benchmark genes.

Percentage of loci at which the benchmark gene (identified from HPO and/or ClinVar) is top-ranked for each gene-based association or gene ranking method. For each method, bars on the left (outlined in black) are calculated for benchmark loci present in both HPO and ClinVar (54 loci), and bars on the right (faded outline) are calculated using the union of all HPO and ClinVar loci (153 loci). Horizontal red lines indicate the expected percentage of top-ranked benchmark genes under the null hypothesis that gene rank and benchmark labels are independent. Error bars indicate 95% confidence intervals. TSS-to-top SNP refers to ranking genes by the distance between TSS and the most significant single variant at each causal locus; the dTSS-weighted gene-based test (dTSS) uses an exponential weight funcion to assign higher weight to variants nearer the TSS for each gene (Methods).

GAMBIT omnibus tests had the highest performance identifying benchmark genes among the gene ranking methods considered, particularly for the stricter gene set, although the difference was not statistically significant relative to most other gene ranking methods (Fig 7). Gene-based tests using coding variants alone had the second-highest performance (Fig 7; S2 Fig), which may reflect the enrichment for coding associations within the benchmark gene set (S3 Fig) caused by benchmark gene selection criteria (described in Materials and methods). Due to the over-representation of coding associations, Fig 7 may underestimate the impact of incorporating heterogeneous regulatory annotations for associated loci without an established benchmark gene.

Further inspection revealed a number of loci of biological or clinical interest. In the analysis of skin cancer in the UK Biobank, three melanin or melanogenesis-related genes (TYR, OCA2, and MC1R) and telomerase reverse transcriptase (TERT) were top-ranked by the omnibus test, but not top-ranked based on TSS-to-top-SNP distance, while all other benchmark genes for skin cancer were top-ranked by both methods or by neither. At the TERT locus, the lead GWAS variant was intronic, whereas the lead variants for TYR, OCA2, and MC1R were nonsynonymous. Unsurprisingly, the latter three benchmark genes were also top-ranked based on coding variant gene-based p-values; however, only TERT was top-ranked based on CT-TWAS.

Similarly, APOB, which encodes an apolipoprotein and is associated with autosomal dominant forms of hypercholesterolemia, was top-ranked by the omnibus test but not by TSS-to-top-SNP distance for disorders of lipoid metabolism in the UK Biobank. Despite being >150 Kbp from the intergenic lead GWAS variant, APOB was also top-ranked by all single-annotation gene-based tests individually. Conversely, TSHR, which encodes a thyroid horomone receptor, was top-ranked based on TSS-to-top-SNP distance but not by the omnibus test for thyrotoxicosis. In this case, the lead GWAS variant was intronic, and CT-TWAS was the only single-annotation gene-based test that ranked TSHR as the top gene at its locus; in this example, while the omnibus test for TSHR was significant, it was outranked by CEP128 at the locus. A complete table of results for benchmark genes is provided in Supplementary Materials.

Linear-type gene-based tests (L-type)

The oldest and most widely used gene-based tests are linear combinations of genotypes across variants [25, 26, 50], here referred to as L-type tests. We define the L-type test statistic as T_L = (w^⊤R_Zw)^−1/2 w^⊤ Z, where w is a vector of single-variant weights, Z is a vector of single-variant association statistics (where each Z_j follows the standard normal distribution under the null hypothesis), and R_Z is the correlation matrix of z-scores. Under the null hypothesis of no association, T_L follows the standard normal distribution. The L-type test statistic T_L can be computed from GWAS summary statistics (single-variant z-scores, or effect sizes and standard errors) and covariance estimates, and can be written either as linear combinations of single-variant association statistics or as linear combinations of genotypes [4, 51].

Examples of L-type tests include burden tests, which calculate burden scores as a weighted sum of rare, putatively deleterious mutations [25, 50]; the cohort allelic sums test (CAST) [52]; and TWAS/PrediXcan tests [6–8], which aggregate eQTL variants using predictive weights estimated from external data sets, e.g. from the GTEx project [40]. These can be viewed as tests of association between GWAS trait and an explicit proxy variable constructed as a linear combination of genotypes. Importantly, L-type tests rely on prior knowledge regarding the directions of effect across variants [27, 50]. For example, the signed weights used in burden tests often reflect the hypothesis that rare deleterious alleles increase risk for disease, and the predictive weights used in TWAS/PrediXcan reflect the hypothesis that gene expression mediates the associations between genotypes and complex trait.

Quadratic-type gene-based tests (Q-type)

Variance component tests and quadratic forms of single-variant association statistics comprise another widely used class of gene-based association methods, here referred to as Q-type (quadratic) tests. Q-type tests include VEGAS (or SOCS), defined as the sum of squared single-variant z-scores [28, 53]; the C-alpha test [54]; and SKAT, a weighted quadratic form of single-variant association statistics [27]. We define the Q-type test statistic as T_Q = Z^⊤ diag(w)Z, where diag(w) is a diagonal weight matrix and Z is a vector of single-variant association z-scores; under the null hypothesis of no association, T_Q follows a mixture chi-squared distribution with mixture proportions equal to the eigenvalues of diag(w)^1/2 R_Z diag(w)^1/2, where R_Z is the correlation matrix of z-scores. In contrast to L-type tests, Q-type tests aggregate single-variant association statistics without prior knowledge or assumptions pertaining to the directions of effects across variants [27, 50]. While less tractable than L-type, analytical p-values for Q-type tests can be calculated using a variety of techniques to approximate the tail probabilities of multivariate normal quadratic forms (e.g., [55, 56]), which are far more efficient than permutation procedures or Monte Carlo methods [28, 57]. Q-type tests are most appropriate when a sizable proportion of variants are hypothesized to have non-zero effects of unknown and inconsistent direction [50].

Maximum chi-squared statistic as a gene-based test (M-type)

Perhaps the simplest gene-based test is the maximum chi-squared statistic across variants (or equivalently, the minimum p-value), here referred to as M-type tests. Analytical p-values for M-type tests can be calculated by directly integrating the multivariate normal density of z-scores within the hypercube given by $$ \mathit{x}\in {\mathbb{R}}^m:{max}_k|x_k|\le {max}_j\left|Z_j\right|$$ where m is the number of variants, or approximated by adjusting the minimum p-value across variants by the effective number of tests [28, 29]. M-type tests are most appropriate when only one or a small fraction of variants are hypothesized to have non-trivial effects. We note that the M-type test accounts for the correlation structure across variants, whereas Tippett’s method [58] assumes independent p-values; thus, the M-type test reduces to Tippet’s method when z-scores are uncorrelated.

Aggregated cauchy association test (ACAT)

The aggregated Cauchy association test (ACAT), a recently proposed method to combine multiple dependent p-values, can be used to construct gene-based tests by transforming single-variant association p-values using the Cauchy quantile and cumulative distribution functions, and computing a p-value

Harmonic mean p-value (HMP)

Another recently proposed method to combine multiple dependent p-values, the Harmonic Mean P-value (HMP; [31]), can similarly be used to construct gene-based tests by weighting p-values from single-variant association tests. The unadjusted HMP p-value is defined

While this statistic can be anti-conservative when directly interpreted as a p-value, Wilson (2019) showed that 1/p_HMP follows a Landau distribution (with scale and location parameters given in Table 1), which can be used to compute an asymptotically exact HMP p-value. The Landau density function is

Generalizations and extensions

The simple forms of gene-based tests described above can be related and combined through a variety generalizations and extensions. Q-type and M-type can both be viewed as special cases of a statistic (∑_j w_j|Z_j|^p)^1/p, which is equivalent to Q-type when p = 2 and to M-type when p → ∞; this generalization has been used, for example, in the aSPU gene-based test [61]. Similarly, Q-type and L-type can both be viewed as special cases of a statistic $$ {\mathit{Z}}^{\top }(\pi \text{diag}\left({\mathit{w}}_1\right)+(1-\pi ){\mathit{w}}_2{\mathit{w}}_2^{\top })\mathit{Z}$$, which is equivalent to Q-type when π = 1 and L-type when π = 0; this generalization has been used, for example, in the SKAT-O gene-based test [50]. Finally, ACAT and HMP can be used to combine p-values across multiple gene-based test forms [30, 31].

Coding variants

Gene-based tests for coding variants are calculated by aggregating variants separately within each coding subclass (e.g., missense, nonsense, and synonymous) using Q-type (by default) test statistics. These stratified tests are then combined across subclasses using a p-value combination procedure (HMP or ACAT) to calculate a coding omnibus test for each gene.

UTR variants

Gene-based tests for UTR variants are calculated by aggregating variants separately within the 3’ and 5’ UTR regions using Q-type (by default) test statistics, and applying a p-value combination procedure (HMP or ACAT) to calculate a UTR omnibus test for each gene.

dTSS weights

One of the most common heuristics to infer likely causal genes at non-coding GWAS loci is to rank genes by distance between their transcription start site (TSS) and the most significant single GWAS variant. This strategy is appealing given the strong enrichment of regulatory variants near TSS.

To incorporate distance-to-TSS (dTSS) and capture association signals at regulatory variants that are not well-annotated in gene-based analysis, we define the dTSS weights for gene k as $$ w_{jk}\left(\alpha \right)=e^{-\alpha \left|d_{jk}\right|}$$, where d_jk is the genomic distance (number of base pairs) between variant j and the TSS for the gene of interest. Larger values of the parameter α confer more weight to variants nearer the TSS. In practice, we only include variants within a specified window (e.g., 500kbp) of the TSS of the corresponding gene. While dTSS weights can be used in any weighted gene-based test (e.g., Q-type tests), ACAT and HMP are particularly well-suited due to their linear computational complexity, as dTSS-weighted tests often involve thousands of variants per gene.

The optimal α value is expected to vary across loci, and likely depends on local gene density and other factors. However, ACAT and HMP can be applied again to calculate omnibus p-values by combining dTSS-weighted gene-based test p-values p_k(α_i) across multiple values α₁, α₂, … [30, 59]. By default, GAMBIT calculates overall dTSS-weighted test statistics by aggregating across α values 10⁻⁴, 5 × 10⁻⁵, 10⁻⁵, 5 × 10⁻⁶.

Enhancer-target gene weights

To capture association signals across regulatory elements that have been assigned to one or more target gene, we weight variants in regulatory elements by element-to-target-gene confidence scores, and aggregate variants for each gene using either ACAT, HMP, or Q-type gene-based test statistics. For example, we define the regulatory-element weighted Q-type test statistic as $$ T_k^{\mathrm{R}}={\sum }_i{\sum }_{j=1}^{m_i}{w}_{ik}{Z}_{ij}^2$$ where m_i is the number of variants in the i^th regulatory element, w_ik is the confidence weight between element i and gene k, and Z_ij is the j^th variant in the i^th regulatory element.

eQTL weights

Given a vector of weights b_kt to predict expression levels for gene k in a given tissue or cell type t as a linear combination of normalized genotypes, we define the z-score TWAS test of association between predicted expression level and GWAS trait as $$ S_{kt}={\mathit{b}}_{kt}^{\top }\mathit{Z}/\sqrt{{\mathit{b}}_{kt}^{\top }{\mathbf{R}}_Z{\mathit{b}}_{kt}}$$ where Z is the vector of single-variant GWAS z-scores and R_Z is the correlation matrix of z-scores.

To aggregate test statistics across multiple tissues or cell-types, which we refer to as Cross-Tissue TWAS (CT-TWAS), we considered three approaches:

Q-type Cross-tissue Test (CT-Q): Calculating the sum of squared tissue-specific test statistics, $$ {\sum }_t{S}_{kt}^2$$, which has a mixture chi-squared distribution under the null hypothesis of no association,

M-type Cross-tissue Test (CT-M): Calculating an analytic p-value for the maximum absolute test statistic max_t|S_kt| using the multivariate normal joint density of tissue- or cell-type-specific test statistics S_k1, S_k2, … under the null hypothesis of no association, and

ACAT or HMP Cross-tissue Test (CT-A or CT-H): Combining tissue- or cell-type-specific p-values p_kt = 2Φ(−|S_kt|) using the ACAT method or HMP respectively.

CT-Q and CT-M require the cross-tissue correlation matrix R_S with elements $$ \left[r_S\right]{}_{t{t}^{\prime }}=\text{corr}({\mathit{b}}_{kt}^{\top }\mathit{Z},{\mathit{b}}_{k{t}^{\prime }}^{\top }\mathit{Z})={\mathit{b}}_{kt}^{\top }{\mathbf{R}}_Z{\mathit{b}}_{k{t}^{\prime }}/\sqrt{\left({\mathit{b}}_{kt}^{\top }{\mathbf{R}}_{\mathit{Z}}{\mathit{b}}_{kt}\right)\left({\mathit{b}}_{k{t}^{\prime }}^{\top }{\mathbf{R}}_Z{\mathit{b}}_{k{t}^{\prime }}\right)}$$, which can be computed in O(m² n + mn²) time where m is the number of tissues or cell-types and n is the number of eQTL variants. By contrast, CT-A and CT-H p-values can be computed in O(m) time, since ACAT and HMP do not require the correlation matrix to be computed. By default, GAMBIT implements CT-A; in our analysis of UK Biobank data, CT-M, CT-H, and CT-A generally perform similarly, while CT-Q tends to detect fewer significant associations.

Combining single-annotation test statistics

In early versions of GAMBIT, we combined gene-based p-values across annotation classes for each gene using standard family-wise error rate (FWER) and false detection rate (FDR) controlling procedures. However, two recently proposed methods, the Aggregated Cauchy Association Test (ACAT; [59]) and Harmonic Mean P-value (HMP; [31]), provide more powerful approaches to combine multiple dependent p-values, and we have therefore implemented both of these methods in GAMBIT. Unlike gene-based tests such as SKAT, which are formulated as parametric tests in a generalized linear mixed model, these p-value combination methods are essentially non-parametric, assuming only that p-values are uniformly distributed under the null hypothesis. Like the exponential combination (EC) procedure [62], ACAT and HMP are powerful under sparse alternatives; however, unlike EC, they enable efficient analytic p-value calculation, with computation time linear in the number of p-values. We calculated ACAT p-values (defined above) using standard formula for the Cauchy quantile function and CDF. To calculate asymptotically exact HMP p-values, we adapted a C++ routine from the ROOT System [63] for computing the Landau distribution CDF following the derivations of Wilson (2019) for the asymptotic distribution of the HMP.

Enhancer-target annotation sources

To identify regulatory genetic elements and their putative target genes, we used pre-computed annotation data sets from three existing methods: Joint Effects of Multiple Enhancers (JEME) [11], GeneHancer [35], and RoadmapLinks [10, 34, 64]. GeneHancer provides a global confidence score between each enhancer element and one or more putative target genes, while JEME and RoadmapLinks provide tissue- or cell-type-specific enhancer-target confidence scores. For the latter two data sets, we calculated overall enhancer-target confidence scores across tissues and cell types as the soft maximum (LogSumExp function) of tissue- or cell-type-specific scores for each enhancer-target pair. Descriptive statistics for each enhancer annotation dataset are provided in S2 Table.

eQTL predictive weight annotation sources

To incorporate eQTL variants in gene-based analysis, we used pre-computed tissue-specific predictive weights for eGene expression estimated using GTEx v7 [47] from TWAS/FUSION (including elastic net and LASSO models) [8] and PredictDB [6, 7]. We generated a GAMBIT eWeight annotation files incorporating all available tissues and cell types for each data resource and predictive model. Descriptive statistics for each eQTL variants weight dataset are provided in S1 Table.

Coding variant and gene annotation sources

We annotated coding variants, TSS locations, and UTR variants using TabAnno 419 [32] and EPACTS [33] based on GENCODE v14 [65].

Simulating GWAS summary statistics

We simulated GWAS traits under the model $$ \mathit{Y}={\mathbf{1}}_n{\beta }_0+\tilde{\mathbf{G}}\mathit{\beta }+\mathit{\epsilon }$$, where $$ \mathit{Y}\in {\mathbb{R}}^n$$ is a quantitative trait for a GWAS sample of size n, 1_n is the n × 1 vector of 1’s and $$ {\beta }_0\in \mathbb{R}$$ is the trait intercept, $$ \tilde{\mathbf{G}}\in {\mathbb{R}}^{n\times m}$$ is the centered and scaled genotype matrix where each column has mean 0 and variance 1, $$ \mathit{\beta }\in {\mathbb{R}}^{m\times 1}$$ is a vector of causal genetic effects, and $$ \epsilon \in {\mathbb{R}}^n$$ is an i.i.d. trait residual with $$ \mathbb{E}\left({\epsilon }_i\right)=0$$ and $$ \text{Var}\left({\epsilon }_i\right)={\sigma }_{\epsilon }^2$$. We scale the parameters $$ {\sigma }_{\epsilon }^2$$ and β so that Y_i has unit marginal variance.

We define the vector of single-variant association statistics (equivalent to t-test statistics from simple linear regression) for variants k = 1, 2, …, m as

We simulated GWAS association statistics Z by calculating $$ \widehat{\mathbf{R}}$$ from the European subset of the 1000 Genomes Project panel, and replacing $$ \widehat{\mathbf{D}}$$ by its limiting value D with elements $$ D_{kk}=1-{\alpha }_k^2$$.

Simulating genetic effects at causal loci

We used empirical functional annotation data to simulate causal genetic effects β, guided by the intuition that a variant’s functional effects ultimately determine its effects on complex traits. While minor allele frequency (MAF) was not explicitly used to select causal variants in simulations, this procedure induces an implicit relationship between MAF and causal status due to the relationship between MAF and functional annotations (S4 Fig). For each simulated causal locus, we selected a causal gene by sampling a single CCDS protein-coding gene, and defined proximal genes as any gene with TSS within 1 Mbp of the causal gene TSS. We then simulated single-variant GWAS summary statistics for all variants associated with any causal and proximal genes by proximity (≤ 1 Mbp) or functional annotations (e.g., eQTL variants).

We simulated causal genetic effects under 5 scenarios: 0) no association (null model), 1) coding association, 2) enhancer association, 3) eGene association, and 4) UTR association. For coding and UTR associations, we first selected the number of causal variants $$ M^{*}={\sum }_{j}I({\beta }_j^2>0)$$ from a Poisson distribution with rate parameter λ = M/4 truncated to 1 ≤ M* ≤ M, where M is the total number of coding (or UTR) variants for the causal gene, and randomly selected M* causal variants from the total set of M coding (or UTR) variants for the causal gene. This procedure results in ~25% of all coding (or UTR variants) having non-zero causal effects, while ensuring that at least one variant is causal. For enhancer associations, we similarly simulated the number of causal enhancers $$ M_e^{*}$$ from a Poisson distribution with rate parameter λ = M_e/4, where M_e is the number of enhancers mapped to the causal gene, and selected causal enhancers using a categorical distribution with probability weights derived from confidence scores between enhancer elements and the causal gene. For eGene associations, we selected a single causal tissue at random, and simulated causal effect sizes proportional to precomputed eQTL weights for the causal gene and tissue. Because eQTL weights are noisy in practice, we used simulated weights $$ \tilde{\mathit{w}}\sim {\mathcal{N}}_{M^{*}}(\mathit{w},\frac{9}{10N}{\widehat{\mathbf{R}}}^{-1})$$ in place of the original weight vector w in TWAS gene-based tests, where N is the GTEx v7 sample size for the causal tissue.

For non-eQTL effects, we simulated the genetic effect for each causal variant β_j from an iid normal distribution, scaled so that the total genetic variance at each locus is equal to $$ h_L^2$$. Because our model has assumed genotypes are scaled with unit variance (and β is scaled accordingly), this simulation approach implicitly assumes that the heritability-scale effect of each causal variant is independent of MAF. This is essentially equivalent to the widely-used model Var(β_j,unscaled) = τ²[2MAF_j(1 − MAF_j)]^a, where a = −1 [26, 66, 67].