Multi-layered genetic approaches to identify approved drug targets

Summary Drugs targeting genes linked to disease via evidence from human genetics have increased odds of approval. Approaches to prioritize such genes include genome-wide association studies (GWASs), rare variant burden tests in exome sequencing studies (Exome), or integration of a GWAS with expression/protein quantitative trait loci (eQTL/pQTL-GWAS). Here, we compare gene-prioritization approaches on 30 clinically relevant traits and benchmark their ability to recover drug targets. Across traits, prioritized genes were enriched for drug targets with odds ratios (ORs) of 2.17, 2.04, 1.81, and 1.31 for the GWAS, eQTL-GWAS, Exome, and pQTL-GWAS methods, respectively. Adjusting for differences in testable genes and sample sizes, GWAS outperforms e/pQTL-GWAS, but not the Exome approach. Furthermore, performance increased through gene network diffusion, although the node degree, being the best predictor (OR = 8.7), revealed strong bias in literature-curated networks. In conclusion, we systematically assessed strategies to prioritize drug target genes, highlighting the promises and pitfalls of current approaches.


In brief
Sadler et al. compared and benchmarked genetically informed approaches combined with network diffusion to prioritize drug target genes. Gene prioritization methods were based on large-scale genetic studies such as genome-wide association and wholeexome studies, as well as tissue-wide and whole-blood expression and protein quantitative trait loci.

INTRODUCTION
Drugs whose targets have genetic support were found to be more likely to succeed in clinical trials. 1,2 Although multiple methods have been proposed to establish such genetic support, leveraging genetic data to find disease genes, and ultimately drug targets, has proven to be challenging. [3][4][5][6] The most straightforward approach maps genome-wide association study (GWAS) signals to the closest genes, with more sophisticated methods incorporating linkage disequilibrium (LD) structure and gene annotation information to compute gene scores. [7][8][9] Over the past decade, large-scale molecular quantitative trait loci (mQTL) datasets facilitated the discovery of disease mechanisms and the identification of potential new drug targets. [10][11][12][13][14][15] Several methods, including Mendelian randomization studies, transcriptome-wide association studies, and colocalization methods have integrated expression and protein QTL data with GWASs to pinpoint likely causal genes for complex traits and diseases. [16][17][18][19][20][21][22] More recently, the availability of high-throughput sequencing data enabled the discovery and analysis of rare variants and their aggregated effects to reveal gene-disease associations. 23,24 Whole-exome sequencing (WES) in the UK Biobank (UKBB) showed that genes prioritized this way are 3.6 times more likely to be targets of drugs approved by the US Food and Drug Administration (US FDA). 25 A B C Figure 1. Overview of the analysis workflow (A) Three different gene prioritization methods were tested in this study. The first one uses GWAS summary statistics as input (GWAS). The second combines molecular QTL and GWAS summary statistics (QTL-GWAS): either expression QTL (eQTL) or protein QTL (pQTL) data. The third leverages individual-level wholeexome sequencing (WES) data (Exome). In the GWAS method, gene p values are based on the sum of squared SNP Z scores (T sum ) that follows a weighted c 2 1 distribution. The QTL-GWAS method integrates QTL and GWAS summary statistics through Mendelian randomization (MR). MR causal effect sizes (b MR ) are calculated from GWAS and mQTL effect sizes (GWAS b and mQTL b, respectively) and gene scores are the corresponding p values. The Exome method aggregates rare variants from WES data. Putative loss-of-function and missense variants with minor allele frequencies (MAF) below 1% are collapsed in burden tests, which results in gene p values. The different approaches were benchmarked for their ability to prioritize drug target genes. the four seed-generating methods and the three networks were applied to 30 traits ( Figure 1) using five different reference sets of target genes (DrugBank, 35 Ruiz et al., 36 ChEMBL, 37 DGIdb, 38 and STITCH 39 ). Overall, we provide an in-depth comparison of all combinations of these approaches, identifying their respective strengths and caveats.

Overview of the analysis
In this study, we calculated gene prioritization scores and tested their ability to identify drug targets across 30 traits (Figure 1). We focused on three types of method, termed GWAS, QTL-GWAS, and Exome, that allow the computation of gene scores provided genetic association data ( Figure 1A).
The GWAS method takes as input GWAS summary statistics together with a matching LD reference panel. Gene p values are calculated based on the sum of squared test statistics for SNPs falling into the gene region. 9 The QTL-GWAS methods integrate GWAS summary statistics with mQTL data for the gene of interest. We calculated gene scores using (1) eQTL-GWAS data from the largest available whole-blood eQTL study (eQTLGen study, n = 31,684) 13 as well as tissue-wide eQTL data from the GTEx Consortium v.8 (n = 65-573 for 48 tissue types) 40 and (2) pQTL-GWAS data from the largest available plasma pQTL study (deCODE study, n = 35,559). 14 Integration was done by performing MR analyses using either the protein or the transcript as exposure and the GWAS trait as outcome.
If not specified otherwise, the eQTL-GWAS method refers to the tissue-wide analysis in which the eQTLGen and GTEx data are combined by considering the tissue for which the MR effect was the most significant (STAR Methods). While the GWAS and QTL-GWAS methods focus on common genetic variants, the Exome method considers only rare variants from WES data with minor allele frequencies (MAFs) below 1%. Gene scores were based on gene burden tests that aggregate putative lossof-function and missense variants, and we used the resulting p values from the WES analysis in the UKBB. 25 To allow for a fair comparison with the Exome method while also exploiting disease-specific consortium GWAS summary statistics with maximized case counts, we calculated gene prioritization scores for the GWAS and QTL-GWAS methods using both consortium GWAS and UKBB GWAS data that matched Exome sample sizes (Tables S1 and S2; STAR Methods).
Disease genes may not coincide with drug target genes, but they may be in close proximity in terms of molecular interaction ( Figure 1B). Through diffusion based on random walks, we leveraged network connectivity to prioritize neighbors of disease genes, which may be drug targets. We tested this hypothesis on three different network types: the STRING PPI network, which relies on literature interactions, among other data types 32 ; a gene coexpression network based on 31,499 RNA-seq samples (CoXRNAseq) 33 ; and a gene coexpression network based on single-cell RNA-seq and proteomics data (FAVA). 34 Gene prioritization scores were obtained following diffusion at six different restart parameter values (r = 0, 0.2, 0.4, 0.6, 0.8, 1) (STAR Methods).
Finally, prioritized disease genes, defined as the top 1% of genes identified through the 12 combinations of gene prioritization and network diffusion methods (5% for combinations involving the pQTL-GWAS method to account for the smaller set of testable genes), were then tested for enrichment with the five drug target genes using Fisher's exact test ( Figure 1C). Background genes were defined as all genes that could be tested by the respective method, and sensitivity analyses were performed on background genes testable for all methods. Second, we calculated the area under the receiver operating characteristic curve (AUC) values, which has the advantage of not requiring any thresholds. To compute a combined enrichment score per method, we aggregated results across traits and drug databases termed overall odds ratios (ORs) or overall AUC values (STAR Methods).

Concordance of prioritized genes among gene scoring methods
We first analyzed whether genes prioritized by the GWAS, QTL-GWAS, and Exome methods were concordant ( Figure 2). For each of the 30 traits, we calculated gene scores for the testable autosomal protein-coding genes (GWAS, $19,150; eQTL-GWAS, $12,550 (blood) and $16,250 (tissue-wide); pQTL-GWAS, $1,870; Exome, $18,800). In the tissue-wide eQTL-GWAS method, the tissue with the most significant MR p value was selected. In Figure S1, we show the proportion of genes mapped to a particular tissue category. The contributions of glandularendocrine, neural central nervous system (CNS), and whole-blood (eQTLGen) tissue categories were the highest (respective means of 15.3%, 12.8%, and 12.6% across the 30 traits ; Tables S3 and  S4). Although each trait had genes mapped to nearly all tissues, a few distinctive patterns could be observed: cardiac muscle tissues contributed the most to atrial fibrillation (16.4%); vascular tissues the most to coronary artery disease (16.5%), followed by diastolic (11.1%) and systolic (9.9%) blood pressure; and the neural CNS the most to schizophrenia (16.9%) and bipolar disease (16.6%).
The concordance of prioritized genes among pairs of methods is summarized in Figure 2. For each trait, we calculated Fisher's exact tests between the top prioritized genes at thresholds ranging in the top 0.1%-10% (STAR Methods). The overlap was the highest between the GWAS and the eQTL-GWAS methods (Figure 2A). At 1%, the median OR was 212.2, which dropped to 51.0 and 22.1 at 5% and 10%, respectively. The overlap of prioritized genes was the lowest with the Exome method. The top 1% GWAS vs. Exome and eQTL-GWAS vs. Exome overlaps (based only on UKBB GWAS summary statistics), yielded median ORs of 1.7 and 1.9, respectively, which dropped to 1.0 at 10% for both methods (Figures 2B and 2C).

Enrichment of prioritized genes for drug targets
Next, we assessed the extent to which prioritized genes overlapped with drug target genes. For each trait, we conducted enrichment analyses for the GWAS, eQTL-GWAS, pQTL-GWAS, and Exome methods using our five definitions of drug target genes.
In Figure 3A, we show the resulting ORs for the DrugBank/ DGIdb database combination. Across methods, genetic support for drug targets was the highest for low-density lipoprotein (LDL)  Figure 2) and using the same background genes and drug target definitions as in (A). (D) ORs calculated for the five drug target definitions and for all four methods (legend in B). The OR was set to 1 for traits with no identified drug target genes. In (C) and (D), the boxplots bound the 25th, 50th (median, center), and 75th quantiles. Whiskers range from minima (Q1 À (1.5 3 IQR)) to maxima (Q3 + (1.5 3 IQR)) with points outside representing potential outliers.  Figure 3B. Restricting disease genes to the top 0.1% for all methods increased the average ORs without changing the method ranking, with average ORs of 3.68, 4.02, 2.40, and 1.44 for the GWAS, eQTL-GWAS, Exome, and pQTL-GWAS methods, respectively. We further analyzed whether identified drug targets were the same across methods and found that prioritized drug target genes were similar between GWAS and eQTL-GWAS methods (average Jaccard index of 0.39), were less so between eQTL-GWAS and pQTL-GWAS methods (blood tissues; average Jaccard index of 0.15), and were very different from Exome identified targets (average Jaccard index of 0.06 between GWAS and Exome and between eQTL-GWAS and Exome methods). Average AUC values across traits were 53.4%, 51.9%, 50.5%, and 49.9% for the GWAS, eQTL-GWAS, Exome, and pQTL-GWAS methods ( Figure 3C). While the number of drugs reported per indication was similar across databases (average of 43.9, 41.8, and 40.4 for Ruiz et al., ChEMBL, and DrugBank, respectively), the average number of reported drug targets was much higher for Ruiz Figures 3D and S2), likely due to the low average number of reported drug targets, which leads to very high ORs when drug targets figured among the prioritized genes (e.g., for LDL and total cholesterol), but for many traits drug target genes were not among the prioritized genes (e.g., for type 1 diabetes, atopic dermatitis, and inflammatory bowel disease).

Examples of drug target prioritization ranks
In Figure 4, we highlight drug targets and their gene prioritization ranks for a few examples (complete list in Table S9). Major antihypercholesterolemic drug targets PCSK9 (evolocumab, alirocumab), HMGCR (statins), and NPC1L1 (ezetimibe) were top ranked by all methods (except for no pQTLs being available for HMGCR and NPC1L1; Figure 4A). HCN4, the target of the antiarrhythmic drug dronedarone, was prioritized as a disease gene for atrial fibrillation only through the GWAS method. Although highly expressed in the atrial appendage and left ventricle of the heart, no eQTL was reported for this gene ( Figure 4B). Several antiepileptic drugs target SCN1A, which was highly prioritized by the GWAS and eQTL-GWAS methods, with the strongest MR effect found in the nucleus accumbens (basal ganglia) of the brain (Figure 4C). The antiplatelet drug dipyrimadole used in the prevention and treatment of vascular diseases such as stroke and coronary artery disease is listed to target 23 genes of the PDE superfamily in ChEMBL. Of these, four (PDE4D, PDE3A, PDE3B, PDE6B) were ranked in the top 1% by the exome method for stroke ( Figure 4D). None of the other methods prioritized any of these 23 genes. For coronary artery disease, another superfamily member (PDE5A) had a low ranking (<2%) by the GWAS and QTL-GWAS methods, supported by solid GWAS and e/pQTL colocalization ( Figure 4E).

Heritability of drug target transcripts and proteins
Previous drug target enrichment analyses have shown that drug target genes are more likely to have lower residual variance intolerance scores (RVISs), i.e., are less tolerant to change. 1 Furthermore, limited overlap between eQTL and GWAS hits has been found, and it has been suggested that GWAS and eQTL genes are under different selective constraints. 41 Hence, under the assumption that drug target genes are more likely to be key (core) GWAS genes, we expected that drug target genes are less likely to harbor QTLs. To test this hypothesis, we assessed whether drug target transcript or protein levels are less amenable to regulation by common genomic variations, which could explain the lower than expected performance of QTL-GWAS approaches.
To this end, we compared the cis heritability of drug target genes vs. non-drug target genes that were measured in the respective studies (i.e., also those with no reported e/pQTLs; STAR Methods), where lower heritability would point toward a negative selection. 42 We conducted the analysis per trait and for each of the five drug target gene definitions; however, we could not observe a clear difference between cis heritabilities of drug target and non-drug target genes ( Figure S5). While this means that we cannot explain why the QTL-GWAS approach does not perform better, it may also imply that drug target genes are not necessarily typical GWAS genes or so-called core genes.
Network diffusion to prioritize drug target genes Finally, we assessed whether network diffusion can identify drug target genes for which there is no direct genetic evidence. Gene scores from prioritization methods defined the initial distribution Genes that were not testable by a given method are reported as NA (no e/pQTL means that the gene was measured, but had no QTL), and a range of ranks (i.e., 1-52) indicates tied p values.
(C) Top plot shows the prioritization ranks of SCN1A (sodium voltage-gated channel alpha subunit 1), a drug target gene of several antiepileptic drugs. Bottom plot shows Mendelian randomization (MR) effects (red dots) with 95% CI (black bars) across tissues in which there was a significant eQTL. (D) Antiplatelet drug dipyrimadole and gene prioritization ranks of its multiple drug targets (a non-exhaustive selection) of the phosphodiesterase (PDE) superfamily.
(E) Top plot shows the gene prioritization ranks of PDE5A, another reported target for dipyrimadole. Bottom plot shows the regional SNP associations (Àlog 10 (p)) with coronary artery disease (CAD; GWAS, green), PDE5A protein (pQTL, red), and PDE5A transcript (eQTL, yellow) (red dashed lines indicate the significance thresholds of the respective SNP association, and gray shading marks the position of PDE5A). Bottom row illustrates the position and strand direction of the genes in the locus.
Cell Genomics 3, 100341, July 12, 2023 7 Article ll OPEN ACCESS p 0 of the diffusion process. This process is regulated by a restart parameter r, whereby lower values result in a stronger diffusion (i.e., genes can be prioritized even when distant from initial disease genes; STAR Methods). The stationary distribution was calculated for six different restart parameters, ranging from no diffusion (r = 1) to complete diffusion (r = 0), and for each of the three networks: the STRING PPI network, 32 an RNA-seq coexpression network (CoXRNAseq), 33 and a coexpression and proteomics network (FAVA). 34 Since the set of testable proteins ($1,870) is enriched for drug target genes (two-sided binomial test: p = 1.3eÀ47 for DrugBank/DGIdb; complete results in Table S15; STAR Methods), the AUC values were artificially inflated upon projecting the gene scores onto the network, and pQTL-GWAS results are hence not discussed.
We further assessed which method's AUC values benefited the most from network diffusion. To allow fair comparison with the Exome methods, we used UKBB GWAS data for the GWAS and eQTL-GWAS methods. Across all diffusion parameters r, overall AUC values were significantly higher for GWAS compared with eQTL-GWAS in the STRING and FAVA network (p diff < 4.45eÀ4), but not any different in the RNA-seq coexpression (CoXRNAseq) network (p diff > 0.05). A nominally significant difference in favor of GWAS compared with Exome was observed only in the STRING network at r values of 0.4, 0.6, and 0.8 (p diff of 0.0262, 7.36eÀ3, and 0.0146, respectively). No statistical differences were observed between the eQTL-GWAS and the Exome method except for a nominally significant difference in favor of eQTL-GWAS at r = 0.2 in the CoXRNAseq network (p diff = 0.0113).
When investigating the network connectivity, we observed that drug target genes were significantly more likely to be hub genes, i.e., to have more connections in the network in comparison with other genes (Figures 5C and S8). This observation was particularly strong in the STRING network (mean log-degree = 13.0 vs. 12.3, p diff = 6.6eÀ284 for DrugBank/DGIdb), but also present in the coexpression networks (D log-degree = 0.064, p diff = 0.011 for CoXRNAseq; D log-degree = 0.3, p diff = 6.6eÀ11 for FAVA). As a consequence, the network's node degree (a gene's number of connections to other genes adjusted by the edge weight) was found to be a good predictor of drug targets, and the best performance was found for the network degree in STRING (overall AUC = 77.6%, overall OR = 8.71). Given this bias, we generated random initial disease gene scores and determined to what extent genetically informed p 0 distributions performed better compared with random p 0 distributions. Although the GWAS, eQTL-GWAS, and Exome methods had significantly higher AUC values compared with random score distributions for any given r value in the STRING network (p diff < 1.62eÀ7; Table S12), the performance of a mildly diffused (r = 0:8) random score (which is unaware of the target disease) performed significantly better than any disease gene prioritization method without diffusion (p diff of 4.18eÀ6, 3.58eÀ10, and 2.10eÀ12 compared with GWAS, eQTL-GWAS, and Exome, respectively). In line with this observation, the network degree was still significantly better than gene prioritization methods at a stronger diffusion of r = 0.2 (p diff of 8.98eÀ6, 9.87eÀ13, and 1.89eÀ11 compared with GWAS, eQTL-GWAS, and Exome, respectively).

Examples of prioritized genes through network diffusion
In the following, we describe several examples for which drug targets figured among the top 1% genes only after network diffusion (complete list in Table S13). Amyloid-beta precursor protein (APP) targeted by the monoclonal antibody aducanumab in the treatment of Alzheimer's disease (AD) was ranked 506 (top 2.7%) prior to and 152 (top 0.8%) after diffusion on the STRING network (r = 0.6; Figure 6A) based on the eQTL-GWAS method. Prioritization was largely influenced by its interacting neighbor apolipoprotein E (APOE), which was the top 5 ranked gene for AD by the eQTL-GWAS method and among the top 6 genes (tied p values) by the GWAS method. Although rare mutations in APP are a known cause of AD, 43 the Exome method did not highly prioritize this gene (>top 10%), likely because of low statistical power due to the younger and healthier nature of the UKBB cohort. Indeed, APP was among the top 1% for the GWAS method, leveraging the AD consortium data, but did not reach the top 10% when restricting the analysis to the UKBB. Tumor necrosis factor (TNF), a drug target in the treatment of inflammatory diseases such as psoriasis, was ranked 1,558th (top 8%; Exome-psoriasis) prior to and 182nd (top 0.98%; r = 0.6) post-propagation in the STRING network ( Figure 6B). While initially the drug target F2 (coagulation factor II, thrombin) for venous thromboembolism (VTE) ranked only in the top 2%, it moved up to the top 1% regardless of the network used for diffusion at r = 0.6 (top 0.9%, 0.6%, and 0.7% for STRING, CoXRNAseq, and FAVA, respectively). In the STRING and CoXRNAseq networks, this boost could largely be attributed to the interacting fibrinogen genes (FGA, FGB, and FGG) that ranked in the top 0.06% ( Figure 6C).  Figure 3B, and at an r value of 0, the gene prioritization rank is based simply on the degree of the network nodes. At r < 1, the background genes are the genes reported in the respective network. The star next to the pQTL-GWAS method signals that the set of testable genes for this method is enriched for drug target genes, and therefore, higher AUC values were obtained when adding background genes with zero-valued initial scores. (B) Odds ratios (ORs) between prioritized genes (top 1%) and drug target genes for each network type and method at different r values across the 30 traits (same drug target and background genes as in A). The OR was set to 1 for traits with no identified drug target genes. (C) Histograms showing the degree distribution of drug target genes and non-drug target genes in each network. The difference in log-degree (D) and the p values from two-sided t tests are shown at the top. In (A) and (B), the boxplots bound the 25th, 50th (median, center), and 75th quantile. Whiskers range from minima (Q1 À (1.5 3 IQR)) to maxima (Q3 + (1.5 3 IQR)) with points above or below representing potential outliers.

Summary of findings
We conducted a comprehensive benchmarking between different genetically informed approaches (GWAS, QTL-GWAS, and Exome) combined with network diffusion to prioritize drug target genes. The strength of our analysis lies in the side-byside comparison of gene prioritization methods that individually have proven to be successful in identifying drug targets. In line with previous reports, we find a 1.3-to 2.2-fold enrichment for drug targets among (the top 1%) prioritized genes. 1,2 Recently, methods have emerged that combine multiple genetic predictors to derive an aggregate score, often using machine-learning techniques. 27,44,45 These scores have demonstrated high enrichment for drug targets but reveal little about underlying molecular mechanisms. Our aim was to disentangle the importance of Article ll OPEN ACCESS the choice of the ground truth (i.e., drug target genes) and the input data (such as mQTLs, WES) in combination with different molecular networks to highlight added benefits while also exposing weaknesses compared with using GWAS data alone.
Comparison of gene prioritization methods Adjusting for differences in background genes and data origins, GWAS yielded higher AUC than eQTL-and pQTL-GWAS, but no significant difference was found with Exome. Genes prioritized by the Exome method were different from those identified by the GWAS and QTL-GWAS methods, which was also reflected in the identified drug targets. While this could imply that rare and common variant genetic architectures are complementary, differences could also be due to power issues. Possibly, with increased sample size, the implicated genes will converge, but the extent to which they can be perturbed by regulatory vs. rare coding variants might remain different. Considering ORs, we lacked the statistical power to claim significant differences between methods, since the number of drug targets among the top 1% prioritized genes can be very low. Overall enrichment ORs for drug targets were 2.17, 2.04, 1.81, and 1.31 for the GWAS, eQTL-GWAS, Exome, and pQTL-GWAS methods, respectively. Although ORs for the pQTL-GWAS method may seem lower, it should be noted that testable proteins (i.e., proteins with pQTLs) accounted for $10% of GWAS-testable genes. On the same background genes, ORs for the tissuewide and blood-only eQTL-GWAS methods were 1.38 and 1.22, respectively. For the AUC metric, no significant difference between eQTL-GWAS and pQTL-GWAS was found. In the method comparisons, we considered multiple drug target gene definitions. The number of targets per drug drastically differed between ChEMBL and the DGIdb or STITCH database due to differences in their construct. Drug target genes in the ChEMBL database are manually curated and should not contain false positives, but it remains debatable whether one should consider only primary or also secondary target genes. For instance, ChEMBL lists only HMGCR as a drug target for statins, whereas the DGIdb database also includes APOA5, APOB, and APOE, among others. For this reason, we considered different databases and present enrichment results for both broad and narrow drug target definitions, as well as aggregates.

Benefits and pitfalls of network diffusion
Network diffusion was beneficial for prioritizing drug target genes with weaker genetic support. A remarkable increase in drug target identification was achieved when using the STRING PPI network. However, this improvement may be due to a circularity in the data generation process, whereby drug target genes are more researched and hence have more chance to be found to interact with other proteins, i.e., they tend to look more hub-like. Although genetically informed gene sets performed better than random ones, the genes prioritized by their node degree in the STRING network resulted in the highest AUC values overall. Thus, care has to be taken when relying on literature-derived gene-gene interactions, as aggressive diffusion will point to the same drug target genes, irrespective of the disease, due to the intrinsic bias stemming from under-and overstudied proteins. While the STRING network resource remains of immense value to identify interacting proteins, non-random missing of network edges leads to a biased network structure, which makes this resource less suitable as input for discovering new drug targets. The improvements made with coexpression networks, which do not suffer from publication/curation biases, were minor in comparison. Although significant with the AUC metric, ORs were not significantly increased with a diffusion of r = 0.6 compared with no diffusion for any of the methods.
Limitations of the study Several limitations should be considered. First, we do not take into account the directionality of therapeutic and genetic effects, i.e., whether the drug is an agonist or antagonist. Although found to be less performant than GWAS, QTL-GWAS methods have the advantage of specifying directionality, as opposed to gene scores from the GWAS approach, which ignores SNP effect directions. Second, the mQTL datasets used cover only a small fraction of possible intermediate traits through which SNPs exert their disease-inducing effects. 46 Third, we focus only on common genetic variants when associating transcript and protein levels. With the advent of coupled rare variant-protein level data, either from populations enriched for rare variants or sequencing data, 14,47 more powerful QTL-GWAS methods are likely to emerge that combine mechanistic insights gained from QTL approaches while capturing rare variant associations previously missed. Fourth, drug target data are sparse, which limits the statistical power in benchmarking analyses. Given the required resources to test a drug target in clinical settings, focusing on top ranking genes is of most interest. This scenario is best described with a threshold that defines highly prioritized genes for enrichment analyses. However, ROC curves that quantify the performance at all prioritization thresholds (i.e., use all data at hand) were better powered to detect subtle differences between methods. Resulting AUC values are relatively low (51%-54%), which may be because ranks of genes with non-significant p values are likely unreliable, but these dominate most of the ROC curve. Related to this, even for low false positive rates, there is room for improvement of the gene prioritization methods.
Combining prioritization methods could increase AUC values, as suggested by the distinct drug target sets identified by GWAS and Exome methods, as could the integration of additional functional genomic annotations. 27,44 Finally, our analysis compares methods using historical drug discovery data as the ground truth. These data are highly biased, with G-protein-coupled receptors being targets of a third of FDA-approved drugs. 48 Many other genes may be effective targets, but have never been tested in clinical trials. Thus, our results may not reflect how well the tested genetic approaches uncover true disease genes, but rather how well they identify targets that were historically prioritized in drug development processes. Since the emergence of robust GWASs, more and more clinical trials are motivated by genetically informed targets. Thus, drug target databases will tend to overlap better with GWAS-inspired genes, leading to artificially higher overlap.

Conclusion
To conclude, we systematically evaluated major gene prioritization approaches for their ability to identify approved drug target Cell Genomics 3, 100341, July 12, 2023 11 Article ll OPEN ACCESS genes. Our analyses highlight the power of harnessing multiple data sources by capitalizing on QTLs for mechanistic insights, sequencing data for rare variant associations, GWASs when mQTL signals are missing, and network propagation to leverage gene-gene interactions.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

GWAS data
We used the largest (to-date), publicly available GWAS summary statistics for each analyzed condition (Table S1). GWAS data came mostly from consortia specific to the respective disease, and were often a meta-analysis comprising the UKBB. Twenty-four out of the 30 conditions were case/control studies, the remaining 6 being continuous traits: diastolic and systolic blood pressure (DBP and SBP, respectively 57 ), low-density lipoprotein and total cholesterol (LDL and TC, respectively 58 ), estimated glomerular filtration rate (eGFR 59 ) and heel bone mineral density ( 58 ) proxying chronic kidney disease (CKD) and osteoporosis, respectively. For four traits with low case count in the UK Biobank (< 20,000; chronic obstructive pulmonary disease (COPD), endometriosis, pneumonia and psoriasis) and no large-scale GWAS meta-analysis available, we performed a meta-analysis between the UK Biobank 58 and FinnGen 60 using METAL. 51

GWAS gene scores
We used PascalX 9,49 to compute gene scores based on GWAS summary statistics. The software takes as input GWAS p values, gene annotations and LD structure. SNPs are assigned to genes and their squared z-scores are summed. This sum, under the null, was shown to follow a weighted chi-square distribution with weights being defined by the local LD structure from which gene p values can be derived. 9 We applied PascalX with default parameters (gene ± 50 kB) on protein-coding genes using the Ensembl identifiers and annotations (Ensembl GRCh37.p13 version) and the UK10K reference panel. 61 Across traits, $19,150 protein-coding genes could be tested which were ranked by their PascalX p value.

Molecular QTL-GWAS gene scores
We integrated molecular quantitative trait loci (QTL) and GWAS summary statistics using Mendelian randomization (MR) implemented in the smr-ivw software. 22,62,50 The exposure (transcript or protein levels) and outcome disease were instrumented with independent genetic variants, also called instrumental variables (IVs; r 2 < 0:01) and used to calculate putative causal effect estimates of the exposure on the outcome (b MR ). IVs were required to be strongly associated to the exposure (P QTL < 1e-6) and had to pass the Steiger filter ensuring no significantly stronger effect on the outcome than on the exposure. 63 We used expression QTLs (eQTLs) from the eQTLGen consortium 13 (whole blood; n = 31,684) and tissue-specific QTLs from the GTEx v8 release 40 (European ancestry; n = 65-573 for 48 tissue types; Table S3) to estimate causal transcript-trait effects. In the eQTLGen dataset there were $ 12,550 protein-coding genes with at least 1 IV which increased to $16,250 when integrating the GTEx dataset. MR results from both datasets (whole blood from eQTLGen and 49 tissues from GTEx) were aggregated by considering the MR causal effect with the lowest p value across tissues (Tables S3 and S4). Protein QTLs (pQTLs) from the deCODE study 14 (whole blood; n = 35,559) were used to estimate protein-trait causal effects with $1,870 proteins having at least 1 IV. Prior to the analysis, e/pQTL and GWAS data were harmonized, palindromic SNPs were removed as well as SNPs with an allele frequency difference > 0.05 between datasets. All transcripts and proteins were mapped to Ensembl identifiers as provided by eQTLGen, GTEx and deCODE.

Exome gene scores
We used gene burden test results computed on WES data from the UK Biobank. 25 We extracted gene-trait associations based on putative loss of function (pLOF) and deleterious missense variants with MAF < 1% (M3.1 nomenclature in original publication) with phenotypes matching the investigated conditions as indicated in Table S2. Associations were provided for $18,800 genes which were ranked by the association p value and retrieved by the provided Ensembl identifier.

Drug target genes
We extracted drug target genes from public resources by combining drug-indication and drug-target links from various databases. A given disease/indication was linked to a drug if the drug was indicated to be prescribed for the selected indication and subsequently, the target genes of these drugs were extracted. For drug-indication pairs we consulted DrugBank, Ruiz 64 using PubChem IDs as intermediates.. d Search Tool for Interacting CHemicals (STITCH) 5.0 39 : Aggregated drug-protein interaction data from high-throughput experiments data, manually curated datasets and prediction methods. Only high confidence drug-protein relationships (confidence score R 700) of the type ''inhibition" and ''activation" were considered. STITCH uses PubChem Chemical Identifiers (CID) for drugs and mapping to DrugBank IDs was done through the chemical sources file provided by STITCH. Protein Ensembl identifiers were mapped to gene Ensembl identifiers using biomaRt (GRCh37, v2.50.3) 54 . d ChEMBL 37 (download: May 2022): ChEMBL provides drug targets which have been manually curated from literature. Drug targets are identified by ChEMBL IDs with mapping to UniProt Accessions provided by ChEMBL. UniProt identifiers were then mapped to gene Ensembl identifiers through the UniProt REST API 65 .
In this analysis we considered drug target genes resulting from the following combinations: DrugBank/DGIdb, DrugBank/STITCH, Ruiz/DGIdb, Ruiz/STITCH, and ChEMBL/ChEMBL. The number of drugs and drug target genes per indication is shown in Table S6.
Transcript and protein level heritabilities Transcript and protein level cis-heritabilities were estimated from QTL effects using a restricted maximum likelihood method (reml) with the LDAK-thin heritability model. The LDAK heritability model assumes that the expected heritability contributed by each SNP depends on its MAF and LD. The analysis was conducted with the LDAK software (v5.2; reml method 55 ) based on all SNPs in proximity of the transcript/protein ( ± 500 kB) and the UK10K reference panel. 61 We set the -power to À0.25 and the -ignore-weights flag to YES to specify the LDAK-thin heritability model. The analysis was restricted to high-quality SNPs which were defined as being nonambiguous, having a sample size > 5,000 and a MAF R 0.01.
Protein heritabilities were based on the deCODE plasma protein dataset 14 and transcript heritabilities for whole blood on the eQTL-Gen dataset. 13 Of the 14,022 protein-coding transcripts in eQTLGen, reml converged for 12,218. Likewise, 3,716 of the 4,502 autosomal proteins in deCODE converged (estimated cis-heritabilities are in Table S10). Genes not converging were omitted in cis-heritability downstream analyses.
To calculate the difference in heritabilities between drug target and non-drug target genes, we considered all transcripts and proteins measured in the respective study which were classified accordingly. Per trait, the difference in heritability was then calculated through a two-sided t-test. Heritability tests were only performed for traits with at least three drug targets within the respective set of measured transcripts/proteins.

Networks
To calculate network diffusion scores, we used the following three networks: d Search Tool for Retrieval of Interacting Genes/Proteins (STRING) v11 32 : The protein-protein (PPI) interaction network results from predictions based on genomic context information, coexpression, text-mining, experimental biochemical/genetic data and curated databases (curated pathways and protein-complex knowledge). Protein Ensembl identifiers were mapped to gene Ensembl identifiers using biomaRt (GRCh37, v2.50.3). 54 We use interaction confidence scores as edge weights. d CoXRNAseq 33 : This network was constructed by first performing a principal component analysis on the gene coexpression correlation matrix of 31,499 RNA-seq samples. Reliable principal components were retained from which the final network was constructed via Pearson correlations. We filtered pairwise interactions to only retain those with z-scores above 4. Genes were reported with Ensembl identifiers and z-scores were used as edge weights. d Functional Associations using Variational Autoencoders (FAVA) 34 : This network is based on single cell RNA-seq read-count data from the Human Protein Atlas and proteomics data from the PRoteomics IDEntifications (PRIDE) database. First, the high-dimensional expression data was reduced into a latent space using variational autoencoders. From this latent space, the network was derived via pairwise Pearson correlations. Each reported interaction has a score which we use as edge weight (final network reports interactions with scores above 0.15). Protein Ensembl identifiers were mapped to gene Ensembl identifiers using biomaRt.