TIGA: target illumination GWAS analytics

2.1 NHGRI-EBI GWAS Catalog preprocessing

The February 12, 2021, release of the Catalog references 11 671 studies and 4865 PubMed IDs. The curated associations include 8235 studies and 2706 EFO (experimental factor ontology) mapped traits. After filtering studies to require (i) EFO-mapped trait; (ii) P-value < 5 × 10⁻⁸; (iii) reported effect size (odds ratio or beta); and (iv) mapped protein-coding gene, we obtained 4118 studies, 1521 traits and 10 264 genes. For consistency between studies, only genes mapped by the Ensembl pipeline (The NHGRI-EBI GWAS Catalog FAQ) for genomics annotations were considered (not author-reported). Figures 1 and 2 illustrate the growth of GWAS research as measured by counts of studies and subjects.

2.2 RCRAS = RCR aggregated score

The purpose of TIGA is to evaluate the evidence for a gene–trait association, by aggregating multiple studies and their corresponding publications. The iCite RCR (Hutchins et al., 2016) is a bibliometric statistic designed to evaluate the impact of an individual publication (in contrast to the journal impact factor). By field- and time-normalizing per-publication citation counts, the RCR measures evolving impact, in effect a proxy for scientific consensus. Hence by aggregating RCRs we seek a corresponding measure of scientific consensus for associations—see (1).

where study = GWAS (study accession), gc = gene count (in study), pub = publication (PubMed ID) and sc = study count (in publication).

The log₂() function is used with the assertion that differences of evidence depend on relative, rather than absolute differences in RCR. Division by sc effects a partial count for publications associated with multiple studies. Since RCR ≥ 0, log₂(RCR + 1) ≥ 0 and intuitively, when RCR = 1 and sc = 1, log₂(RCR + 1) = 1. Similarly, division by gc reflects a partial count since studies may implicate multiple genes. This approach is informed by bibliometric methodology, including fractional publication counts, as described elsewhere (Cannon et al., 2017). For recent publications lacking RCR, we used the global median as an estimated prior. Computed thus, RCRAS extends RCR with similar logic, providing a rational bibliometric measure of evidence for scoring and ranking gene–trait associations.

2.3 Association weighting by SNP–gene distance

Mapping genomic variation of single nucleotides (SNPs) to genes is a challenging area of active research (Lamparter et al., 2016; Liu et al., 2010; Mishra and Macgregor, 2015). This project does not contribute to mapping methodology. Rather, TIGA employs mappings provided by the Catalog between GWAS SNPs and genes, generated by an Ensembl pipeline that ‘adds additional SNP-specific information associated with the rsID extracted … This information is retrieved using the Ensembl API and the source of the data is both Ensembl and NCBI’ (The NHGRI-EBI GWAS Catalog GWAS Catalog Curation). It is important to note that this method is unbiased and derived from experimental data and the current human reference genome. TIGA aggregates SNP-trait associations, assessing evidence for gene–trait associations, based on these understandings:

SNPs within a gene are more strongly associated than SNPs upstream or downstream.

Strength of association decreases with distance, or more rigorously stated, the probability of LD between an SNP and protein-coding gene decreases with genomic physical distance. Accordingly, we employ an inverse exponential scoring function, consistent with LD measure (Δ) and coefficient of decay (β) by Wang et al. (2006).

This function, used to weight N_snp to compute a distance-weighted SNP count N_snpw, is plotted together with the observed frequencies of mapped gene distances in Supplementary Figure S1, to illustrate how the extant evidence is weighted—see (2).

where d = distance in base pairs and k = ‘half-life distance’ (50k).

2.4 Multivariate ranking

Multivariate ranking is a well-studied problem which needs to be addressed for ranking GWAS associations. We evaluated two approaches, namely non-parametric μ scores (Wittkowski and Song, 2010) and meanRank, and chose the latter based on benchmark test performance. meanRank aggregates ranks instead of variables directly, avoiding the need for ad hoc parameters. Variable ties imply rank ties, with missing data ranked last. We normalize scoring to (0,100] defining meanRankScore as follows.

Variables of merit used for scoring and ranking gene–trait associations:

• N_snpw: N_snp weighted by distance inverse exponential described above.

• pVal_mLog: max[-Log(pValue)] supporting gene–trait association.

• RCRAS: RCR aggregated score (iCite RCR-based), described above.

Variables of merit and interest not currently used for ranking:

• OR: median (odds ratio, inverted if <1) supporting gene–trait association.

• N_beta: simple count of beta values with 95% confidence intervals supporting gene–trait association.

• N_snp: SNPs involved with gene–trait association.

• N_study: studies supporting gene–trait association.

• study_N: mean (SAMPLE_SIZE) supporting gene–trait association.

• geneNtrait: total traits associated with the gene.

• traitNgene: total genes associated with the trait.

N_snp, N_study, geneNtrait and traitNgene are counts of the corresponding unique entities. From the variables selected via benchmark testing the meanRankScore is computed using (3):

where rank_i = rank of ith variable and N = number of variables considered.

μ scores were implemented via the muStat (Wittkowski and Song, 2010). Vectors of ordinal variables represent each case, and non-dominated solutions are cases which are not inferior to any other case at any variable. (For TIGA, cases are genes or traits, corresponding with trait queries or gene queries, respectively, and their variables of merit described above.) The set of all non-dominated solutions defines a Pareto boundary. A μ score is defined simply as the number of lower cases minus the number of higher, but the ranking is the useful result. The ranking rule between case k and case k′ may be formalized as in (4). Simply put, case k′ is higher than case k if it is higher in some variable(s) and lower in none

2.5 Benchmark against gold standard

Lacking a suitable gold standard set of gene–trait associations in general, we instead relied on established gene–disease associations from the genetics home reference (Fomous et al., 2006) and UniProtKB (The UniProt Consortium, 2018) databases. This gold standard set was built following a previously described approach (Pletscher-Frankild et al., 2015). It consists of 5366 manually curated associations (positive examples) between 3495 genes and 709 diseases. All other (2 472 589) possible pairings of these genes and diseases were considered negative examples.

To assess the quality of the TIGA gene–trait associations, we mapped the Ensembl gene IDs to STRING v11 identifiers using the STRING alias file (Szklarczyk et al., 2019) and the EFO terms to disease ontology identifiers (Schriml et al., 2019) based on ontology cross-references and the EMBL-EBI ontology Xref service. We then benchmark any individual variable or multivariate ranking of the associations by constructing the receiver operating characteristic (ROC) curve by counting the agreement with the gold standard.

3.1 The TIGA web application

TIGA facilitates drug target illumination by currently scoring and ranking associations between protein-coding genes and GWAS traits. While not capturing the entire Catalog, the TIGA app can aggregate and filter GWAS findings for actionable intelligence, e.g. to enrich target prioritization via interactive plots and hitlists (Fig. 3), allowing users to identify the strongest associations supported by evidence.

Hits are ranked by meanRankScore described in Section 2. Scatterplot axes are effect (OR or N_beta) versus evidence as measured by meanRankScore. This app accepts ‘trait’ and ‘gene’ query parameters via URL, e.g. ?trait=EFO_0004541, ?gene=ENSG00000075073, ?trait=EFO_0004541&gene=ENSG00000075073. Gene markers are colored by target development level (TDL) (Oprea et al., 2018). TDL is a knowledge-based classification that bins human proteins into four ordinal and non-overlapping categories: Tclin, mechanism-of-action designated targets via which approved drugs act (Avram et al., 2020; Santos et al., 2017; Ursu et al., 2019); Tchem are proteins known to bind small molecules with high potency; Tbio includes proteins that have gene ontology (Ashburner et al., 2000) ‘leaf’ (lowest level) experimental terms; or meet two of these conditions: a fractional publication count (Pafilis et al., 2013) above 5, 3 or more gene ‘Reference Into Function’ annotations (Mitchell et al., 2003), or 50 or more commercial antibodies in Antibodypedia (Björling and Uhlén, 2008); Tdark are manually curated UniProtKB proteins that fail to place in any of the previous categories.

3.2 Benchmark against gold standard disease–gene associations

To benchmark the quality of the GWAS associations in TIGA, we focused on the 383 EFO terms that could be mapped to diseases and their 20 458 associations with genes. We evaluated the performance of each variable of merit individually against the manually curated gold standard gene–disease associations.

The resulting ROC curves showed that the best performing variables are RCRAS, N_study, pVal_mLog, N_snpw and N_snp, which have areas under the curve (AUC) higher than 0.6 (Fig. 4A and Supplementary Fig. S1). The three variables RCRAS, pVal_mLog and N_snpw are furthermore complementary, having a maximal pairwise Spearman correlation of 0.325, whereas N_study and N_snp are strongly correlated with the better performing RCRAS and N_snpw, respectively. We thus used these three variables as the basis for calculating two multivariate rankings, namely μ score and meanRankScore. We benchmarked both rankings the same way as the individual variables and found that μ score performs marginally better than meanRankScore score based on their AUC values (Fig. 4B). However, as the meanRankScore outperforms the μ score in the area of interest [0.0, 0.2] and is more than five orders of magnitude faster to calculate, we selected it as the final ranking in TIGA. Corresponding plots for the lower performing variables are provided in the Supplementary materials for completeness.

3.3 Using TIGA for drug target illumination

The main motivation of developing TIGA is to capture GWAS data when illuminating drug targets. Table 1 shows how many targets from each protein family and Illuminating the Druggable Genome (IDG) TDL are covered with associated traits in TIGA, with families as defined by drug target ontology (Lin et al., 2017) level 2. IDG TDL is a knowledge-based classification: Tclin = mechanism-of-action drug targets (Santos et al.); Tchem = small-molecule modulators are known; Tbio = biological function elucidated; Tdark = understudied protein (Oprea et al., 2018). Coverage for the understudied 2469 Tdark proteins is of particular interest. However, the data for other TDLs can also provide unique and complementary evidence, especially in case of Tbio proteins that are biologically characterized but have not been clinically validated.

Figures 3 and 5 illustrate a typical use case, the plot and gene list for trait ‘high-density lipoprotein cholesterol measurement’, which monitors blood levels of high-density lipoprotein cholesterol as a risk factor for heart disease. Figure 6 shows the provenance for one of the associated genes, GIMAP6 ‘GTPase IMAP family member 6’ with the scores and studies for this gene–trait association, including links to the Catalog and PubMed. GIMAP6 is an understudied (Tbio) member of the GTPases of immunity-associated protein family (GIMAP). Although literature-based evidence of its role in cholesterol homeostasis is scarce (Hoffmann et al., 2018; Richardson et al., 2020), this finding is substantiated by significantly increased circulating HDL cholesterol levels in GIMAP6-knock-out female mice (https://bit.ly/3uznvCU), suggesting that loss of GIMAP6 function may be linked with hypercholesterolemia-associated disorders.

Figures 7 and 8 illustrate another target illumination example, for trait ‘HbA1c measurement’ (glycated hemoglobin, signifying prolonged hyperglycemia), highly relevant to the management of type 2 diabetes mellitus (Rahbar et al., 1969; Saudek and Brick, 2009). Figure 8 shows the provenance for one of the associated genes, SLC25A44 ‘Solute carrier family 25 member 44’ with the scores and studies for this gene–trait association, including links to the Catalog and PubMed. SLC25A44 is an understudied (Tdark) branched-chain amino acid (BCAA) transporter that acts as metabolic filter in brown adipose tissue, contributing to metabolic health (Yoneshiro et al., 2019) and may be involved in subcutaneous white adipose BCAA catabolism (Lee et al., 2021).

4.1 Target illumination

The explicit goal of the NIH IDG program (Oprea et al., 2018) is to ‘map the knowledge gaps around proteins encoded by the human genome’. TIGA is fully aligned with this goal, as it evaluates the GWAS evidence for disease (trait)–gene associations. TIGA generates GWAS-centric trait–gene association dataset using an automated, sustainable workflow amenable for integration into the Pharos portal (Nguyen et al., 2017; Sheils et al., 2021). The Open Targets platform (Ghoussaini et al., 2021; Ochoa et al., 2021) uses Catalog data and other sources, assisted by supervised machine learning, to identify probable causal genes, and validate therapeutic targets by aggregating and scoring disease–gene associations for ‘practicing biological scientists in the pharmaceutical industry and in academia’. OpenTargets associations are enhanced, yet limited by the training data and knowledge sources reflecting current understandings of genetics. In contrast, TIGA provides aggregated evidence solely from the Catalog, reflecting experimental results with minimal bias, thus interpretable in terms of provenance and methodology, and more suitable for some downstream consumers and use cases.

4.2 From information to useful knowledge

In data-intensive fields such as genomics, specialized tools facilitate knowledge discovery, yet interpretation and integration can be problematic for non-specialists. Accordingly, this unmet need for integration and interpretation requires certain layers of abstraction and aggregation, which depend on specific use cases and objectives. Our target audience is drug discovery scientists for whom the aggregated findings of GWAS, appropriately interpreted, can provide additional value as they seek to prioritize targets. This clear purpose serves to focus and simplify all aspects of its design. Our approach for evidence aggregation is simple, easily comprehensible, and based on what may be regarded as axiomatic in science and rational inductive learning: First and foremost, evidence is measured by counting independent confirmatory results.

Interpretability concerns exist throughout science, but GWAS is understood to present particular challenges (Gallagher and Chen-Plotkin, 2018; Lambert and Black, 2012; Marigorta et al., 2018; Visscher et al., 2017). The main premise of GWAS is that genotype–phenotype correlations reveal underlying molecular mechanisms. While correlation does not imply causation, it contributes to plausibility of causation. Genomic dataset size adds difficulty. The standard GWAS P-value significance threshold is 5 × 10⁻⁸ based on overall P-value 0.05 and Bonferroni multiple testing adjustment for 1–10 million tests/SNPs (Marigorta et al., 2018). The statistical interpretation is that the family-wise error rate, or overall probability of a type 1 error, is 5%, but associations to mapped genes require additional interpretation. Motivated by, and despite these difficulties, it is our belief that GWAS data can be rationally interpreted and used by non-specialists, if suitably aggregated. Accordingly, TIGA is a rational way to suggest and rank research hypotheses, with the caveat that the identified signals may be accompanied by experimental noise and systematic uncertainty.

4.3 Designing for downstream integration

Biomedical knowledge discovery depends on integration of sources and data types which are heterogeneous in the extreme, reflecting the underlying complexity of biomedical science. These challenges are increasingly understood and addressed by improving data science methodology. However, provenance, interpretability and confidence aspects are underappreciated and rarely discussed. As in all signal propagation, errors and uncertainty accrue and confidence decays. Here, we proposed the use of simple, transparent and comprehensible metrics to assess the relative confidence of disease–gene associations, via the unbiased meanRank scores. Figure 9 summarizing TIGA sources and interfaces, illustrates its well-defined role. Continuous confidence scores support algorithmic weighting and filtering. Standard identifiers and semantics support rigorous integration. Limiting provenance to the Catalog and its linked publications, semantic interpretability is enhanced.