Abstract
Drug repurposing may provide a solution to the substantial challenges facing de novo drug development. Given that 66% of FDA-approved drugs in 2021 were supported by human genetic evidence, drug repurposing methods based on genome wide association studies (GWAS), such as drug gene-set analysis, may prove an efficient way to identify new treatments. However, to our knowledge, drug gene-set analysis has not been tested in non-psychiatric phenotypes, and previous implementations may have contained statistical biases when testing groups of drugs. Here, 1201 drugs were tested for association with hypercholesterolemia, type 2 diabetes, coronary artery disease, asthma, schizophrenia, bipolar disorder, Alzheimer’s disease, and Parkinson’s disease. We show that drug gene-set analysis can identify clinically relevant drugs (e.g., simvastatin for hypercholesterolemia [p = 2.82E-06]; mitiglinide for type 2 diabetes [p = 2.66E-07]) and drug groups (e.g., C10A for coronary artery disease [p = 2.31E-05]; insulin secretagogues for type 2 diabetes [p = 1.09E-11]) for non-psychiatric phenotypes. Additionally, we demonstrate that when the overlap of genes between drug-gene sets is considered we find no groups containing approved drugs for the psychiatric phenotypes tested. However, several drug groups were identified for psychiatric phenotypes that may contain possible repurposing candidates, such as ATC codes J02A (p = 2.99E-09) and N07B (p = 0.0001) for schizophrenia. Our results demonstrate that clinically relevant drugs and groups of drugs can be identified using drug gene-set analysis for a number of phenotypes. These findings have implications for quickly identifying novel treatments based on the genetic mechanisms underlying diseases.
Over the last 15 years, the use and effectiveness of genome-wide association studies (GWAS) have developed rapidly, resulting in the discovery of over 270,000 phenotypically associated genetic variants1,2. While translating GWAS findings into use for pharmacological and clinical purposes has been challenging, seminal work by Nelson et al. and King et al. demonstrated that drugs with genetic support for efficacy were 2-5 times more likely to gain approval status compared to compounds with no such evidence3,4. In fact, 66% of FDA-approved drugs in 2021 are supported by human genetic evidence that link the compounds to their gene product targets5. GWAS results have been used as early as 2009 to prioritize novel targets for pharmaceuticals6. For example, Wang et al identified the IL12/IL23 pathway as a novel target for Crohn’s disease, which supported the subsequent development of Ustekinumab, an antibody that binds with IL12 and IL236,7.
However, the development of novel pharmaceuticals is a time and financially intensive process, averaging 12.8 years and $2.6 billion dollars for a single new medication8,9. Despite wide advances across science and technology in the past 60 years, the cost for approved drugs per billion spent has seen an 80-fold increase after adjusting for inflation10. Additionally, the development of novel drugs for psychiatric and neurological conditions has been particularly affected, with several large pharmaceutical companies shutting down research and development for psychiatric disorders entirely11. Given the size of the challenges facing de novo drug development, interest has turned towards drug repurposing, the process of identifying new clinical indications (the phenotype a drug is approved to treat) for currently approved medications12. Moreover, the vast quantities of GWAS results provide excellent data regarding possible gene targets to be leveraged for drug repurposing. To date, many methods have been created to repurpose drugs using GWAS summary statistics, including candidate gene approaches, Mendelian randomization, pharmagenic enrichment scores, transcriptomic signature matching, and drug gene-set analysis13–17.
Drug gene-set analysis, an approach first published by De Jong and Breen18, can be used to identify candidate drugs based on GWAS results. Gene-set, or pathway, analysis involves mapping genetic variants (e.g., single nucleotide polymorphisms [SNPs], copy number variants [CNVs]) to genes, and grouping genes into sets based on, for example, a priori biological functions or co-expression networks, and testing each gene set for association with the phenotype19,20. For drug gene-set analysis, gene sets are created for drugs based on known biological and/or chemical properties. For example, a drug-gene set can be made of the genes that code the proteins targeted by a specific drug, and this gene set can be tested for association with a phenotype. While De Jong & Breen found no gene-sets to be significantly associated with schizophrenia, their nominally significant results included the known schizophrenia drug trifluoperazine, and metoclopramide and neratinib, which target dopamine receptors and tyrosine kinases (both implicated in schizophrenia)18. This method has been continued by Gaspar and colleagues, who developed the drug repurposing tool Drug Targetor, which is largely based on drug gene-set analysis17,21.
However, while drug gene-set analysis has been tested in schizophrenia and bipolar disorder, it has never been verified in non-psychiatric phenotypes and, to our knowledge, no drugs have ever been found to be significantly associated with schizophrenia. Additionally, previous research testing groups of drugs for associations have used methods that violated statistical assumptions of independence17. Here, we apply drug gene-set analysis to eight phenotypes – hypercholesterolemia (i.e., high cholesterol), type 2 diabetes, coronary artery disease, asthma, schizophrenia, bipolar disorder, Alzheimer’s disease, and Parkinson’s disease – and use a novel statistical approach to test if groups of drugs are enriched for genetic signal associated with a phenotype while accounting for the covariance between drugs. Our objectives were first, to assess the ability of drug gene-set analysis at prioritizing clinically relevant drugs for each phenotype, and second, to identify candidate drugs for repurposing for the phenotypes investigated.
Results
We extracted data for drug-gene targets and interactions from the Clue Repurposing Hub and the Drug Gene Interaction Database and created a gene set for each drug (n = 1201) composed of the genes whose protein products are known to be targeted by or interact with the drug in question (Figure 1A)22,23. To test drug-gene set associations, competitive gene-set analysis was performed using the MAGMA software tool for each of the eight phenotypes19 while conditioning on a gene set of all drug target genes in our data (n = 2281; Methods). This was done to ensure that significant drug-gene set associations were due to effects unique to drug pathways, and not due to common properties shared by all drug target genes (e.g., many gene products targeted by drugs are receptors). For the individual drug-gene sets see Supplementary Table S1.
(A) Gene sets were created for every drug (n = 1201) consisting of the genes whose protein products are targeted by or interact with that drug. Since similar drugs share gene targets it is possible for overlap of genes to occur between drugs. MAGMA was used for gene-set analysis to test for significant relationships between drug-gene sets and each phenotype. Drug-gene sets were tested using a conditional competitive gene-set analysis (Methods). (B) Drug groups were tested for enrichment of genetic signal using a modified multiple linear regression model. This was done for ATC III groups (n = 85), mechanism of action drug groups (n = 79), and clinical indication drug groups (n = 118) with n > 5 (Methods). (C) Lastly, gene sets were created for each subcategory of each grouping method (i.e., ATC III codes, clinical indications, and mechanisms of action) and tested for association using the same MAGMA protocol as drugs in panel A (Methods).
Individual drug-gene set results
We identified 22 significant drug-trait associations across the eight phenotypes we investigated using drug gene-set analysis (Table 1 & 2; Fig 1A). Drugs were significantly associated with each of the phenotypes except for asthma and Parkinson’s disease. For the non-psychiatric/neurological phenotypes, many of our results included drugs approved to treat the diseases in question. For instance, all 4 of the drugs associated with type 2 diabetes – glipizide, glyburide, mitiglinide, and repaglinide – are approved for type 2 diabetes24–26 (Tables 1 & 2). Similarly, statin drugs (e.g., pravastatin and atorvastatin) were associated with both coronary artery disease and hypercholesterolemia, which are used to lower LDL cholesterol and prevent myocardial infarctions27. Likewise, lomitapide and fenofibrate, which were also associated with coronary artery disease, are drugs used to treat hypercholesterolemia, a known risk factor for coronary artery disease28–30. We found no significant associations between approved psychiatric drugs and the psychiatric phenotypes that were tested. However, the antiepileptic drug gabapentin was significantly associated with schizophrenia, which is known to have some efficacy in treating psychiatric disorders31. For all individual drug gene-set analysis results see Supplementary Tables S2 – S9.
Drug group enrichment
Drug gene sets were grouped in three ways: by Anatomical Therapeutic Classification (ATC) III code, by clinical indication, and by the mechanism of action. Previous research that tested for enrichment of genetic signal within groups of drug gene sets used drug enrichment curves derived from a Wilcoxon Mann Whitney U test17. However, due to the overlap of genes between drug gene sets, the assumption of data independence needed for a Wilcoxon Mann Whitney U test is not valid. Thus, here we instead employed a multiple linear regression approach, modeling the drug group membership as a predictor of the MAGMA drug gene-set t-statistic, to test whether drugs in that drug group exhibited (on average) more strongly associated drug gene sets than other drugs. Overlap between gene sets is accounted for through the residual covariance matrix of the regression model.
Each drug group with 5 or more drugs within the ATC (n = 85), mechanism of action (n = 79), and clinical indication (n = 118) categories was tested for enrichment of genetic signal using the method above. Bonferroni correction was applied by correcting for the number of drug groups tested within each of the three grouping methods, and then using a family-wise correction for the number of grouping methods used. Therefore, the Bonferroni-corrected p-value threshold was (0.05/85) / 3 = 0.00019 for ATC drug groups, (0.05/79) / 3 = 0.00021 for mechanism of action drug groups, and (0.05/118) / 3 = 0.00014 for clinical indication drug groups. For all drug groups see Supplementary Tables S10 – S12, and for the drug group enrichment results see Supplemental Table S13 – S21.
As seen in Figure 2 (top), we identified 14 drug groups that were significantly enriched for genetic signal across the 8 phenotypes we examined. For the non-psychiatric/neurological phenotypes, we found at least one drug group containing one or more drugs approved for the phenotype tested. For hypercholesterolemia, we identified two drug groups where more than 50% of the drugs were approved for hypercholesterolemia (C10A and HMGCR inhibitors), and we identified the mechanism of action group insulin secretagogues for type 2 diabetes, which contains only drugs approved for type 2 diabetes. None of the drug groups that were associated with our psychiatric and neurological phenotypes contained approved drugs. However, the ATC code N07B – drugs used for addictive disorders – that is associated with schizophrenia contains several drugs (e.g., naltrexone, buprenorphine) that have shown significant improvements in treating both positive and negative symptoms in the disorder32. As seen in the bottom panels of Figure 2, none of the clinical indication drug groups for the phenotypes we examined were significantly associated with our phenotype after correction for multiple tests. The clinical indication groups for Parkinson’s disease, type 2 diabetes, and hypercholesterolemia were nominally enriched for genetic signal for each of their respective phenotypes.
Drug groups were tested for enrichment of genetic signal with each phenotype. The left panels show the −log10 p-values for each drug group, and the right panels show the percentage of approved drugs in the group. The top panels show all drug groups that are significant after multiple testing correction. The bottom panels portray the drug groups containing all approved drugs for each of the 8 phenotypes that were tested. TCL = hypercholesterolemia; T2D = type 2 diabetes; CAD = coronary artery disease; AST = asthma, SCZ = schizophrenia; BIP = bipolar disorder; PKD = Parkinson’s disease; ALZ = Alzheimer’s disease. See Supplementary Tables S13 – S14 for full results.
To examine the genes driving the significant results seen in the top panel of Figure 2, we generated heatmaps of the drug-gene sets for the most significant drug groups for each phenotype; insulin secretagogues for type 2 diabetes (Fig 3A), the clinical indication drug group of urticaria for asthma (Fig 3B), ATC code J02A for schizophrenia (Fig 3C), and ATC code A04A for Alzheimer’s disease (Fig 3D). For the C10A heatmaps associated with hypercholesterolemia and coronary artery disease see Supplemental Figures S1 and S2.
The heatmaps depict the drug-gene sets from drug groups that were significantly associated with a phenotype, with color signifying the gene’s −log10 p-value from MAGMA. Panel (A) depicts drug-gene sets from the insulin secretagogues group associated with type 2 diabetes, panel (B) illustrates the drug-gene sets from the clinical indication group urticaria associated with asthma, panel (C) shows the drug-gene sets from ATC code J02A associated with schizophrenia, and panel (D) portrays the drug-gene sets from ATC code A04A associated with Alzheimer’s disease. Drug gene set names are listed on the y-axis, and gene names are listed on the x-axis. Color = −log10 p-value from gene association tests in MAGMA.
As seen in Figure 3A, the genes with the highest MAGMA −log10 p-values targeted by insulin secretagogues are KCNQ1, TCF7L2, PPARγ, IGF2BP2, and KCNJ11. All these genes are known to be associated with type 2 diabetes and affect biological processes like insulin secretion and fat storage in adipose tissue33,34. For the clinical indication drug group urticaria (Fig 3B), there are several genes with evidence linking them to asthma. Namely, research has indicated that the PPARδ gene targeted by desloratadine may be a promising therapeutic target for asthma, with desloratadine itself – an allergy medication – having shown efficacy in treating seasonal allergic asthma35,36. Moreover, another study found that EHMT2 was differentially expressed in patients with asthma compared to healthy controls37. The targets of drugs in the ATC J02A group also contain genes related to schizophrenia. Antigen levels for DPYD were found to be elevated in individuals with schizophrenia compared with healthy controls, and it is known to be associated with schizophrenia via previous genome-wide association studies38,39. Moreover, there is evidence indicating that splicing variation of CYP2D6 may play a role in schizophrenia as well40. Lastly, the ATC group A04A associated with Alzheimer’s disease targets a number of serotonergic receptors, which have shown promise as therapeutic targets in animal models by reducing levels of amyloid-beta41. Interestingly, all of the drug groups depicted appear to be driven by a variety of genes, as opposed to a small number of genes targeted by all drugs in a group (with the exceptions of PPARγ, KCNJ11, and ABCC8 for the insulin secretagogues).
Next, we compared the gene targets of approved drugs for each phenotype to further investigate why drug-gene set analysis was better at identifying approved drugs for the non-psychiatric/neurological phenotypes that were tested. As seen in Figure 4, there are 38 gene targets of approved drugs (8.44% of drug target genes) that are significant after Bonferroni correction for the non-psychiatric/neurological phenotypes investigated, compared to 17 significant gene targets (3.21% of drug target genes) for drugs approved for the psychiatric and neurological phenotypes. Moreover, each of the four non-psychiatric/neurological phenotypes has at least 4 significant genes targeted by approved drugs, while bipolar disorder and Parkinson’s disease have only 1 significant gene targeted by an approved drug, and Alzheimer’s disease has none. This may explain why none of the significant drug groups associated with the psychiatric and neurological phenotypes contained any drugs approved specifically for those disease and disorders.
The top panel depicts the gene-level Manhattan plot for the targets of approved drugs for the four non-psychiatric/neurological phenotypes (n = 450), and the bottom panel depicts the same for the psychiatric/neurological phenotypes (n = 530). Genes that are labeled are significant after Bonferroni correction for their respective phenotype (since the number of genes tested varied per phenotype, different Bonferroni thresholds were applied to each phenotype; see Methods). Color = phenotype; shape = even and odd chromosomes; TCL = hypercholesterolemia; T2D = type 2 diabetes; CAD = coronary artery disease; AST = asthma, SCZ = schizophrenia; BIP = bipolar disorder; PKD = Parkinson’s disease; ALZ = Alzheimer’s disease.
Merged drug gene-set analysis
Some drug-gene sets have a very small number of gene targets, resulting in lower power to detect associations with phenotypes20. Additionally, many drug categories in our data do not contain enough drugs to test for enrichment. Therefore, to be able to include drug targets with small gene set sizes and to test association for drug categories with low numbers of drugs, we also created gene sets by grouping drugs again by ATC III code, clinical indication, and mechanism of action. For each category within the grouping methods, a single gene set was created made up of all the gene targets for drugs in the group (Figure 1C), resulting in 164 ATC III code gene sets, 515 clinical indication gene sets, and 354 mechanism of action gene sets (Tables S22 – S24). Competitive gene-set analysis was performed with the same protocol used when testing individual drug-gene sets for each of the three types of grouped drug-gene sets (Methods). We identified 32 significant categorical gene set-trait relationships, with every phenotype but asthma and Parkinson’s disease having at least one associated gene set. For all significant categorical drug gene sets see Supplementary Table S25, and for all categorical drug-gene set results see Supplementary Tables S26 – S49.
Drug similarity networks
For each phenotype, we created network graphs to investigate how the top drugs per phenotype cluster based on their Jaccard similarity coefficient42 (Figure 5), which is calculated using the proportion of gene targets shared by both drugs (see Methods). Central nodes in each network graph were defined using all drugs with a MAGMA competitive gene-set analysis p-value < 0.01 per phenotype. All drugs with a Jaccard similarity of 0.25 or greater with a central node were included as peripheral nodes. Figure 5A displays the similarity network graph for drugs associated with hypercholesterolemia and Figure 5B illustrates the similarity of drugs associated with schizophrenia.
Drug similarity was calculated using the Jaccard index to measure the proportion of shared gene targets between every pair of drugs. Drugs nominally associated with our phenotypes (p < 0.01) were used as central nodes, with peripheral nodes being any drug with a Jaccard similarity > 0.25 with any central node. Panel (A) displays the results for hypercholesterolemia, and (B) shows the network for drugs associated with schizophrenia. Triangle = approved drug for the phenotype, square = Bonferroni-significant drug that is not approved for the phenotype, edge width = Jaccard similarity coefficient.
As may be expected, statin drugs approved for hypercholesterolemia tend to cluster together, while the other two significant drugs, lomitapride and probucol, have no similar compounds (Figure 5A). There is also a cluster of PPAR agonists around the significant drug saroglitazar and nominally significant drug bezafibrate, which both target the PPAR α/γ subunits. Figure 5B demonstrates that almost none of the Bonferroni-significant drugs identified using drug gene-set analysis for schizophrenia cluster with approved schizophrenia drugs, with the exception being one approved drug sharing an edge with benidipine. Most nominally significant drugs for schizophrenia (Figure 5B) appear in the same large cluster, which ranges from antiepileptics in one subcluster (e.g., ethosuximide), calcium channel blockers in the middle (e.g., felodipine), and ACE inhibitors on the other side (e.g., quinapril). The exception is gabapentin, which appears by itself connected to only one other drug, spironolactone. Interestingly, spironolactone has been found to improve schizophrenia symptoms in mice by acting as an antagonist to the NRG1-ERBB4 signaling pathway, and is currently being tested in humans as an add-on treatment for SCZ43. The rest of the approved schizophrenia drugs (triangles) appear in the cluster centered around amisulpride, azaperone, domperidone, and perospirone. For network graphs of all phenotypes see the Supplemental materials (Figures S3 – S8).
Discussion
In the current study, we applied drug gene-set analysis to four non-psychiatric/neurological phenotypes and four psychiatric/neurological phenotypes, to investigate the performance of drug gene-set analysis at identifying clinically relevant drugs (i.e., drugs approved to treat the phenotypes we examined), and to identify potential repurposing candidates. We found the first drugs whose targets were significantly enriched for genes associated with hypercholesterolemia, type 2 diabetes, coronary artery disease, schizophrenia, and Alzheimer’s disease using this method, and identified drugs with the same mechanism of action (calcium channel blockers) as those published in the latest genome-wide association study of bipolar disorder44. In all, 22 drug-trait associations were identified across the eight phenotypes we examined. Notably, 4 of the 4 drugs identified for type 2 diabetes are approved to treat type 2 diabetes, and 3 statin drugs were significantly associated with both coronary artery disease and hypercholesterolemia, respectively24– 26,45,46. We also demonstrated that drug gene-set analysis does not have to be used to isolate single drugs for repurposing and can also pinpoint categories of drugs by using two different methods: by testing for the enrichment of individual drugs in a category and by combining gene targets of all drugs in a category into a single gene set, thereby also pointing into the direction of putative novel drug development that fall within an identified category. Specifically, we identified drug groups enriched for genetic signal for hypercholesterolemia, type 2 diabetes, coronary artery disease, and asthma that contained approved drugs for those diseases. Additionally, we found that the clinical indication drug groups for hypercholesterolemia, type 2 diabetes, and Parkinson’s disease were nominally enriched.
Phenotypes that were not enriched for drugs that are already approved for treatment may not reflect a shortcoming of drug gene-set analysis, but more so the gap between the gene products targeted by current treatments and the genes that are most associated with it as measured by genome-wide association studies. Thus, while no approved drugs for the psychiatric/neurological phenotypes examined here were identified, this may be because current treatments for these diseases only target disease symptoms and not the underlying disease mechanisms47–50. The results presented here establish that drug gene-set analysis can identify clinically relevant drugs – with accuracy varying across the phenotypes we examined – and is able to offer new drug candidates for repurposing.
Across our drug gene-set analysis results there is evidence for repurposable drugs for several phenotypes. For instance, we find repeated evidence for the potential treatment of schizophrenia with calcium channel blockers. The calcium channel blockers felodipine, amlodipine, and benidipine were significantly associated with schizophrenia after Bonferroni correction. The use of calcium channel blockers in psychiatry has shown mixed efficacy in the literature. Previous research has tested the use of calcium channel blockers for psychiatric disorders and found them ineffective, although many of these studies used first-generation drugs (e.g., verapamil)51. Recent research has found that calcium channel blocker use is associated with reduced illness severity for schizophrenia and lower rates of psychiatric admission and self-harm during psychosis52,53. Additionally, Lintunen et al. found that calcium channel blocker use was associated with lower risk for psychiatric hospitalization in individuals with schizophrenia, an effect which was specific to the calcium channel blocker subclass of dihydropyridines54. This is a significant distinction, as drugs in this subclass are also known to have higher selectivity and permeability across the blood brain barrier, perhaps making them better candidates for psychiatric use than first generation calcium channel blockers51. Notably, all calcium channel blockers that were identified for both schizophrenia and bipolar disorder, felodipine, amlodipine, benidipine, manidipine, and nilvadipine, are all dihydropyridines55. Likewise, one of the two drugs identified using drug gene-set analysis in the 2021 genome-wide association study on bipolar disorder was nisoldipine, which is also a dihydropyridine calcium channel blocker44. Thus, further research into the use of dihydropyridines in psychiatry may be warranted.
Despite the promising results found here, drug gene-set analysis still leaves several questions unanswered regarding repurposable drugs prioritized with this method. First, drug gene-set analysis does not specifically test for direct associations between a drug and a phenotype. Rather, it tests whether the gene targets of drugs are enriched for genetic signal that is associated with a phenotype. Second, and perhaps most importantly, is that drug gene-set analysis does not indicate the direction of effect, meaning that drugs identified could either alleviate or exacerbate disease progression and symptoms. Additionally, drug gene-set analysis provides no information on the way a drug interacts with its targets (e.g., agonism versus antagonism). These two limitations mean that drug gene-set analysis may best be used in the first steps of a drug repurposing pipeline, or as one of several methods used to triangulate novel drugs for a phenotype56. Drugs prioritized with drug gene-set analysis could be followed up using genomics approaches that consider gene expression (e.g., transcriptomic signature matching) or more statistical approaches like Mendelian randomization15 or local genetic correlation analysis57, using protein levels of the drug-gene targets for each gene in a drug gene set. This may prove tedious for drugs with large gene sets and may fail to capture interaction effects between different proteins. However, transcriptomic signature matching (e.g., correlating gene expression induced by a drug with the gene expression associated with a phenotype) could be an efficient follow-up for large drug-gene sets. Whilst the limitations above are important to note, they do not outweigh the benefits. Given the excessive financial and time challenges facing psychiatric drug development, a cheap in-silico, and quick method for genetically informed drug repurposing, such as drug gene-set analysis, may prove invaluable. The drug gene-set analysis pipeline used here is publicly available (Methods) and can be used with any user provided GWAS summary statistics.
Materials & Methods
Data
Clue Repurposing Hub
The Clue Repurposing Hub is a publicly accessible, annotated data set of FDA-approved, clinical trial, and pre-clinical drugs. We extracted gene targets for every FDA-approved drug from the Clue Repurposing Hub that had a known clinical indication, disease category, and mechanism of action (n = 1201).
Drug Gene Interaction Database
To augment the sets of gene targets downloaded from the Clue Repurposing Hub, we queried the Drug Gene Interaction Database API for additional genes targeted by or interacting with each compound in our data set. This database is publicly accessible.
GWAS
Summary statistics for hypercholesterolemia, coronary artery disease, and Parkinson’s disease were downloaded from GWAS Atlas58–60 (https://atlas.ctglab.nl/). Summary statistics for type 2 diabetes were downloaded from the DIAGRAM Consortium website (http://diagram-consortium.org/downloads.html)33, summary statistics for asthma were retrieved from GWAS Catalog (https://www.ebi.ac.uk/gwas/)61, and those for bipolar disorder and schizophrenia were accessed from the Psychiatric Genomics Consortium website44,61,62 (https://www.med.unc.edu/pgc/download-results/). Lastly, GWAS summary statistics for Alzheimer’s disease (excluding 23andme data) were downloaded from the CTG Lab website63 (https://ctg.cncr.nl/software/summary_statistics).
Analyses
Drug gene-set analysis
Drug gene-sets were created using the gene targets and interactions retrieved from both the Clue Repurposing Hub and Drug Gene Interaction Database for every compound, resulting in 1201 drug gene sets with n >= 2. Additionally, we created 3 types of drug gene-set by grouping drugs by clinical indication (n = 515), by mechanism of action (n = 354), and by ATC III code (n = 164). Competitive gene-set analysis was performed separately with MAGMA v1.09b for each type of drug-gene set, using the 1K Genomes European reference panel19,64. Competitive gene set analysis was computed while conditioning on a gene set of every gene target in our data set (n = 2281) to assess whether drug gene-set associations were specific to the drug-gene set or driven by common properties of drug target genes.
Drug group associations
Groups of drugs were examined for enrichment of genetic signal using a multiple linear regression model to estimate the effect of drug group membership on MAGMA t-statistics. Given that the MAGMA t-statistic is a measure of the degree to which a drug gene set is enriched for genetic signal for a phenotype, this allows us to test whether drugs in that drug group exhibited (on average) greater enrichment than other drugs.
Denoting the vector of drug gene-set t-statistics as Y and encoding a binary indicator variable Gj for drug group j, scored 1 for drugs in that drug group and 0 otherwise, we fit the model Y = β0 + CβC + GjβG + ε = Xβ + ε. Here, C is a covariate matrix containing the drug gene-set size as well as the log of the drug gene-set size, and ε is the residual term with distribution ε ∼ MVN(0, σ2S), where S is the sampling correlation matrix of the MAGMA t-statistics. Combining the intercept and predictor variables in X, the model can be estimated using Generalized Least Squares, yielding the estimate, , with sampling covariance matrix
.For the residual variance we have the estimate
, with
and with N the number of drugs and K = 4 the number of parameters in the regression equation.
For each drug group j, we test the null hypothesis H0 : βG = 0 against the one-sided alternative HA : βG> 0 using the test statistic , where WG is the diagonal element of W corresponding to
and T has a t-distribution with N − K degrees of freedom. Correction for multiple tests was done using a two-stage Bonferroni correction, correcting for the three grouping methods used and then for the number of groups tested within each grouping method. This resulted in a threshold of 0.05/Ng/3 for each grouping method g with Ng the number of groups for g, where Ng = 85 for the ATC III code drug groups, Ng = 79 for the mechanism of action drug groups, and Ng = 118 for the clinical indication drug groups.
Manhattan Plots of Genes Targeted by Approved Drugs
The Bonferroni-significant thresholds for the gene-level Manhattan plots were derived by 0.05/Nx,y, where Nx is the number of genes tested in MAGMA for a phenotype y. The Bonferroni-significant thresholds are as follows: phypercloterolemia = 0.05/17493 = 2.86E-06, pdiabetes = 0.05/11631 = 4.29E-06, pcad 0.05/18180 = 2.75E-06, pasthma= 0.05/18315 = 2.73E-06, pschizophernia = 0.05/18265 = 2.74E-06, pbipolar= 0.05/18262 = 2.74E-06, pparkinsons = 0.05/18367 = 2.72E-06, palzheimers= 0.05/18409 = 2.72E-06.
Drug Similarity Network Graphs
Network graphs were created to investigate how the top results cluster based on the similarity of their gene targets. The similarity of each pair of drugs was evaluated using the Jaccard similarity coefficient, which is calculated using the formula for drug gene sets A and B, where |A ∩ B| denotes the number of genes shared by those gene sets and |A ∪ B| the number of genes that occurs in at least one of the two gene sets. Central nodes for each phenotype were defined as drug-gene sets with a MAGMA gene-set analysis p-value < 0.01, and all drugs with a Jaccard index >= 0.25 with a central node were included as peripheral nodes. Network graphs were created in Python using the package networkx 2.6.3.
URLs
CTG Lab, https://ctg.cncr.nl/software/summary_statistics; GWAS Atlas, https://atlas.ctglab.nl/; GWAS Catalog, https://www.ebi.ac.uk/gwas/; Drug Gene Interaction Database, https://www.dgidb.org/; DIAGRAM Consortium, http://diagram-consortium.org/downloads.html); Psychiatric Genomics Consortium, https://www.med.unc.edu/pgc/download-results/; Clue Repurposing Hub, https://clue.io/repurposing; MAGMA, https://ctg.cncr.nl/software/magma.
Code availability
The software pipeline used to compute drug gene-set analysis and the drug group analyses is publicly available and can be accessed here: https://github.com/nybell/drugsets.
Data Availability
All data used in this study are publicly available and can be found using the links provided in the manuscript.
https://github.com/nybell/drugsets
https://ctg.cncr.nl/software/summary_statistics
http://diagram-consortium.org/downloads.html
Acknowledgments
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 964874. DP is supported by the Netherlands Organization for Scientific Research - Gravitation project ‘BRAINSCAPES: A Roadmap from Neurogenetics to Neurobiology’ (024.004.012), and DP and EU are funded by the European Research Council advanced grant ‘From GWAS to Function’ (ERC-2018-ADG 834057). NYB, EU, EvW, and DP have no conflicts of interest to declare. CdL was funded by F. Hoffman-La Roche AG.
Footnotes
Revision includes updated Supplemental Tables and the proper text citations to reflect these changes. Also added an explicit Code Availability section.