Integrating systems biology sources illuminates drug action.

There are significant gaps in our understanding of the pathways by which drugs act. This incomplete knowledge limits our ability to use mechanistic molecular information rationally to repurpose drugs, understand their side effects, and predict their interactions with other drugs. Here, we present DrugRouter, a novel method for generating drug-specific pathways of action by linking target genes, disease genes, and pharmacogenes using gene interaction networks. We construct pathways for more than a hundred drugs and show that the genes included in our pathways (i) co-occur with the query drug in the literature, (ii) significantly overlap or are adjacent to known drug-response pathways, and (iii) are adjacent to genes that are hits in genome-wide association studies assessing drug response. Finally, these computed pathways suggest novel drug-repositioning opportunities (e.g., statins for follicular thyroid cancer), gene-side effect associations, and gene-drug interactions. Thus, DrugRouter generates hypotheses about drug actions using systems biology data.

Introduction

Pathways form the basis of our understanding of how cellular processes occur, and provide a framework for inferring cellular phenotypes. Drug research and development has provided powerful medications over the last several decades (1). However, our understanding of the therapeutic effects of the drugs, their side effects and drug interactions is still limited by incomplete knowledge of the underlying cellular pathways through which drugs act. For many applications, including drug discovery, drug repurposing and the definition of pharmacogenomic modulators, we need a molecular-level understanding of drug effects and this is often either missing or incomplete.

We focus here on inferring the pathways of interacting biological macromolecules that modulate drug response. By generating drug-specific pathway hypotheses, we reduce the search space and enable researchers to focus their experimental efforts on the most promising directions. The primary challenge for building accurate pathways is our inadequate understanding of gene interactions, both their location and temporal dependencies. Thus, straightforward network algorithms applied to gene interaction data sets yield a very high rate of false positives when they are used to connect drug targets to the genes that produce end-phenotypes. Previous methods manage network noise (in the context of drug response) by using only curated cellular pathways from public databases (2–4) or by constructing pathways using only short paths prevalent in multiple drugs (5). These methods tend to ignore cross talk between pathways, or focus only on pathways that are common to multiple drugs.

Here, we borrow an analogy from roads and traffic in which gene interactions (protein-protein, metabolic and transcriptional) are roads and traversing the network is akin to finding the quickest route between points of interest. Network interactions that are part of a curated biological pathway have higher credibility than other gene interactions, and are considered “highways”. The less reliable and un-curated connections are viewed as “side roads”. Our method, DrugRouter, adopts a conservative strategy that assembles drug-specific pathways by which ‘highways’ are used preferentially and ‘side roads’ are used only when the highways do not connect the desired starting and ending points. The inputs to our method are genes and gene products (henceforth called genes for brevity) of three classes related to a particular drug of interest: (1) the drug’s target genes, (2) the drug’s pharmacogenes that are known to modulate its mechanism of action (i.e. genes whose variation influences drug response), and (3) the genes associated with the drug’s therapeutic effect or disease target. DrugRouter selects robust paths that connect these three sets of genes to one another; the genes that are visited during this ‘tour’ are then assumed to be relevant to the molecular drug response. We focus on the action of drugs (pharmacodynamics, PD) and not their metabolism (pharmacokinetics, PK-also an important area) by excluding pharmacokinetic genes before applying our algorithms. Figure 1 illustrates the steps of our method.

Figure 1

Illustration of the DrugRouter method. Method input, including drug targets, pharmacogenes and disease genes and the network constructed of highways and side-roads (A), building the raw pathway, connecting drug targets, pharmacogenes and disease genes (B) and pruning the raw pathway by keeping only paths traversed by more than one color (C).

We show that the pathways we construct are useful for four applications: (1) elucidating drug-specific PD pathways, (2) suggesting alternative indications for a drug (drug repositioning), (3) associating genes with drug side effects, and (4) associating genes with drug-drug interactions. We validate each of these applications independently.

Results

Drug-pharmacodynamic pathways as perturbed cellular pathways

A key assumption of our method is that drug-related pathways of action are chiefly drawn from existing knowledge of biology, and do not represent uncharacterized cellular interactions. Existing pathway databases reflect current knowledge of cellular mechanisms that are studied for many reasons, including their relevance to basic metabolism and disease processes. It is important for us to demonstrate that these databases contain useful and relevant knowledge for inferring the mechanism of action of drugs. Indeed, known PD pathway display significant overlap with other curated cellular pathways ( supplementary material, section 1).

Building drug-specific PD pathways

We constructed pathways for 113 drugs with at least one known drug target, pharmacogene, and disease-associated gene (Methods, sections 1–4, supplementary material, data file S1). The pathways display high variability in size (90±75 genes on average). Our pathway computations are robust: when we removed individual pharmacogenes or disease-associated genes, we found little overall change in the inferred pathways (see Methods, section 4 and Figure 2).

Figure 2

Robustness of inferred pathways to removal of pharmacogenes and disease-genes. The mean and standard deviation of the Jaccard score between the genes in the inferred pathway and in the leave-one-out pathway (intersection size divided by union size) as a function of the total number of pharmacogenes and disease genes per drug.

We clustered the drugs according to the overlap of the genes in their computed pathways (Methods, section 4). Figure 3 displays clusters of minimal size four, spanning 67 drugs. As expected, the clusters reflect similar or related therapeutic families. We provide an illustration and analysis of a small-sized pathway of tiludronate in the supplementary material (section 2 and figure S1).

Figure 3

Clustering of drugs based on drug-specific pathway similarity. The x-axis and y-axis show the same 67 clustered drugs (same ordering for both axes). Colors correspond to the relative overlap of the constructed pathways for each drug pair (Jaccard score). The number in parenthesis beside the drug names correspond to the cluster number. Labels on the left are “rough guides” based on the plurality of drugs in that cluster, but acknowledge that some drugs may be in unusual clusters and this represents interesting hypotheses about their pathway connections that may deserve follow-up.

We demonstrate the validity of our pathways with three observations made on the inferred genes of the pathway (genes in the drug-specific pathway, excluding the drug targets, pharmacogenes and disease genes): (1) genes in our drug-specific pathway have significant co-occurrence in the literature; (2) our drug-specific pathways have significant overlap with available gold-standard PD pathways, and (3) our pathways are enriched for genes that are hits in drug-based genome-wide association studies (GWAS) of warfarin, paclitaxel and gemcitabine. To assess the significance of these measures, we compared the drug-specific pathways against a set of randomly constructed control pathways (Methods, section 4).

For the support of literature co-occurrences, we queried PubMed and PubMed Central (PMC) for textual co-occurrence of drugs and genes (Methods, section 4). We found that inferred genes significantly co-occur with the corresponding drug compared to unrelated controlled genes (Wilcoxon rank sum test, 109/113 drugs below FDR of 0.05). The obtained p-values were also lower than p-values obtained compared to other pathways (88% of the pathways) and to random pathways (Wilcoxon rank sum test, p< e−120), where 87% of the individual pathways performed better than random pathways and the remaining fifteen pathways were all small (14±10 inferred pathway genes as compared to 75±59 for other drug-specific pathways, p<2e−7 - see also Discussion).

As a second verification, comparison to gold standard PD pathway, we focused on 43 drugs for which 21 PD pathways are available in PharmGKB. 23 of the drug-specific pathways (corresponding to 11 known PD pathways) were enriched for overlapping genes (hypergeometric enrichment, false discovery rate (FDR) <0.05). Moreover, when we measured the distance of inferred genes in the drug-specific pathways to the known PD pathways, the remaining 20 drug-specific networks were significantly closer (in terms of network distance) to the known PD pathways than the random control pathways (FDR<0.05).

The third evaluation, focused on GWAS for warfarin (6), paclitaxel (7) and gemcitabine (8), is described in the supplementary material (section 3). Additionally, we describe associations between genes modulating response to doxorubicin in yeast and the doxorubicin-specific pathway genes in the supplementary material (section 4).

Applying inferred pathways to suggest drug repositioning

We can use the drug-specific pathways to propose new indications for existing drugs. Using two independent criteria, we obtained 1,484 pairs of novel drug-disease pairs using the first and 1,348 pairs using the second (Methods, section 5). 195 of the pairs satisfied both criteria. These predictions associate 113 drugs to 139 diseases (Table S1).

We evaluated our repositioning predictions by computing enrichment of the predictions in (1) current clinical trials in phases I–III and (2) off-label uses extracted from electronic health records (Methods, section 5). Between 15% and 27% of the experimental (i.e. not yet approved) drug and disease associations being tested in clinical trials satisfy the first and second criteria, respectively (hypergeometric test, p<2e−7 and p< 2e−36, respectively). Similarly, between 19%–21% of off-label drug uses, extracted from electronic health records (9) satisfy one of the two criteria (hypergeometric test, p<5e−4 and p<0.02, respectively). Last, our predictions were enriched in a ‘silver-standard’ set of drug-disease associations (p<0.02 and p<0.05 for the two criteria, respectively) (Methods, section 5). Notably, our predictions were also enriched with the phenotype-driven repositioning prediction set published in (10) (p<0.04 and p<2e−7, respectively).

Another evaluation of our repositioning focused on repositioning of cancer drugs to other cancer types. We downloaded bioassays from Pubchem (11) for 34 cancer drugs and mapped the bioassays to 30 cancer diseases (supplementary methods, section 2). Aggregating the number of active and inactive bioassays across all studies and excluding inconclusive treatments resulted in 361 drug-disease pairs, where the drug was exclusively active or inactive. We found that only the second criteria received significant results (hypergeometric test, p<0.05. odds-ratio=1.57).

Among the predicted drug repurposing opportunities, we highlight five that were predicted using both criteria. First, conjugated estrogens and verapamil are predicted to affect type 2 diabetes mellitus. Indeed, conjugated estrogens improve glycemic control in postmenopausal women with type 2 diabetes (12) and short-term verapamil decreases fasting plasma glucose and glucose turn-over in non-insulin dependent diabetics, possibly by inhibition of gluconeogenesis (13). Second, thalidomide may be useful against Alzheimer’s disease (AD), in part through effect on the AD-associated genes - amyloid precursor protein (APP) and amyloid beta A4 precursor protein-binding family B member 2. Indeed, a recent study showed that long-term treatment of thalidomide may treat a mouse model of Alzheimer’s disease through inhibition of beta-site APP cleaving enzyme 1 (14). Third, NSAIDs, including aspirin and the withdrawn drug rofecoxib, may have a role in colon cancer. Indeed, the former was reported to reduce the risk of colon cancer (15) while the latter is chemo-preventive in a mouse model of colon cancer (16). The last two examples show evidence of a drug-disease association, but the type of association (treatment or causal) is less evident at this time: (1) valproic acid may be associated with acute leukemia. Indeed a study reported this association with proposed mechanism of inhibition of histone deacetylase (17), while a case report for three patients claimed that valproic acid causes leukemia through the same mechanism (18) and (2) statins on our list (Simvastatin, Lovastatin, Pravastatin, Atorvastatin and Fluvastatin) are associated with follicular thyroid cancer. While studies reported that lovastatin induces apoptosis in a different thyroid cancer - human anaplastic thyroid carcinoma cells (19, 20), a single case study reported development of thyroid follicular adenoma on simvastatin therapy (21) and another study reported the development of thyroid neoplasms at high dosages in rats (22). A nine-year follow-up reported inconsistent findings, where increased risk of thyroid cancer in men after 5 years of statin use was not supported by the 2-year lag results or by the findings in women (23).

Several drugs predicted to interact with diseases may not treat them, but instead cause them or increase their severity. Aspirin (acetylsalicylic acid) induces asthma (24) and calcium or paclitaxel may increase the risk of myocardial infarction (25, 26). These predictions, while not leading to new repositioning may be, however, useful in exposing the potential molecular mechanisms underlying these disorders.

Inferred pathways suggest novel associations between genes, side effects and drug interactions

In a similar manner to the drug repositioning, the existence of a gene associated with a side effect (SE) on a drug pathway may indicate that the drug induces the side effect. Using a literature curated list of gene-side effect associations, 26% of the predicted side effects are known drug side effects (hypergeometric test, p<3e−15) (supplementary material, section 6).

As the set of known genes associated to side effects are currently limited, our drug-specific pathways can suggest novel associations between genes and potential side effects. We computed the enrichment of genes within their pathways for 764 SEs associated with the 113 drugs, accounting for similar drugs. We obtained a final list of 135 gene-SE associations, spanning 33 genes and 50 SEs (supplementary material, section 6, Table S3 and Figure S2). Finally, our drug-specific pathways can associate genes and PD drug interactions. We found fifteen genes enriched for co-occurrence in drug-specific pathways of severely interacting drugs (supplementary material, section 7, Table S4 and Figure S3).

We evaluated our gene-SE predictions by querying PubMed and PMC for associations between 28 SEs from our prediction set and the genes in our network and our gene-DDI predictions by querying PubMed and PMC for the co-mentioning of the genes with all possible drug combinations. Nineteen SEs (out of 28) were significantly co-mentioned with the predicted associated genes and ten (out of fifteen) genes had more frequent co-mentioning with severely interacting drug pairs than with non-interacting drug pairs (FDR <0.05 for both). We provide detailed analysis and highlight examples in the supplementary material, section 6 and 7.

Discussion

In this paper, we introduce a novel method for inferring drug-specific pathways. We connect known drug associated genes (drug targets, pharmacogenes and genes associated with diseases treated by the drug) over protein, metabolic and transcriptional interaction networks while preferring high confidence interactions participating in curated cellular processes. In that sense, our method is conservative and unlikely to propose radically new pathways unless the high-throughput evidence is very strong. Upon evaluation of our constructed pathways, we were able to suggest novel drug repositioning, associate genes with side effects and suggest potential causes for drug-drug interactions.

In order to reduce the inherent noise in biological networks, our algorithm searches for the most confident paths, using preferentially high confidence “highways” and maintaining only highly traversed paths. Our conservative approach emphasizes precision to control for false positives. However, we might miss additional genes that reside either on less traversed pathways or on “side roads” with these requirements. This property may have caused the genes on fifteen small pathways to have fewer co-mentioning with the drug than the genes in the random pathways. In addition, most protein-protein interactions lack directionality as well as sign (activation/inhibition). Incorporating of such additional information into the model would enhance the pathway construction task.

Not all the drug-specific pathways were enriched in curated PD pathways, but were significantly closer to those curated PD pathways on the interaction network. The human curation requirements enforced by PharmGKB curators lead to a very high specificity of these curated pathways, but low sensitivity, potentially missing several lower evidence parts.

We assumed that discovery of new disease genes or drug target genes along a pathway propose a new drug repositioning opportunity. We disregarded, however, some factors that may prevent such an opportunity from materializing such as the actual role of those disease or drug target genes in the treatment (e.g. activation vs. inhibition) as well as the expression of those genes in a given tissue. Indeed some of our repositioning suggestions were found to induce or elevate the disease risk.

As noted by several authors (27, 28), drug discovery is fast moving from the single gene research paradigm to the systems biology analysis paradigm. DrugRouter represents a general-purpose tool to harness pathway information for multiple uses.

Methods

1. Data sets

Gene interaction network

Protein-protein interactions (PPIs) were assembled from the union of BioGrid ver. 3.1.94 (29), DIP (Aug 2012) (30), HPRD release 9 (31), IntAct (Oct 2012) (32), MINT (Oct 2012) (33), MIPS (34) and HIPPIE ver. 1.4 (35). Curated PPIs were extracted from KEGG (May 2012) (36) human signaling pathways. Metabolic interactions between enzymes were extracted from KEGG (36) human metabolic pathways. Protein-Gene (transcriptional) interactions were retrieved from the ChEA database (37). The network includes more than 223,000 interactions.

Pathways

Pathways were imported from the Pathway Interaction Database (38) which includes the NCI-Nature human curated pathways and selected pathways from Reactome (39) and BioCarta (40) (1331 Pathways). Pharmacodynamic pathways were downloaded from PharmGKB (41) (41 pathways).

Drug-specific genes

Drug targets were downloaded from DrugBank (42). Drug sensitivity variants were retrieved from PharmGKB (41) (data file S2).

Drug-disease genes

Genes associated with drug indications were assembled by mapping disease-associated genes retrieved from OMIM (Apr 2013) (43) to drug indications retrieved from (44). On average, 18.5±18.3 disease genes were assigned to a drug (data file S2)

Drug side effects and drug interactions

Drug-side effect associations were downloaded from the SIDER2 database (45). DDIs were retrieved from DrugBank (42) and the drugs.com website (46) as described in (10).

2. Building the network of highways and side roads

We constructed a network by integrating three types of interactions: (1) PPIs, (2) Metabolic, where enzymes are connected via a mutual metabolite, and (3) transcriptional, where a transcription factor is connected to a transcribed gene. Each network interaction was tagged as ‘highway’ if it is a (1) a PPI appearing in a KEGG signaling pathway, (2) a PPI where the two interacting genes appear in the same curated cellular pathway, or (3) a curated metabolic interaction from KEGG metabolic pathways. The remaining interactions were tagged as ‘side roads’. Transcriptional interactions and ‘highway’ interactions from KEGG are directed and the rest are undirected (data file S2).

3. The DrugRouter algorithm

DrugRouter constructs pathways in two consecutive stages: (1) A construction stage and (2) A pruning stage.

The first, construction, stage connects pairs of start and destination points (e.g. a drug target as a starting point and a pharmacogenes as the destination point). We included five different start-destination pair types: (1) Drug targets to pharmacogenes, (2) drug targets to disease genes, (3) pharmacogenes to pharmacogenes, (4) pharmacogenes to disease genes, and (5) disease genes to disease genes.

For each pair type, we connect all the start and destination pairs by applying three steps:

Locate all the nearest highway entry points (“on-ramps”), i.e minimal number of side-roads between the start point and the highway entry point. If the distance to the nearest “on - ramp” is farther than one standard deviation above the mean network path length (more than three interactions in our network), that start point is excluded from the analysis.

Locate all the highway exist points (‘off-ramps’) nearest to the destination point in a similar manner to step 1.

Find the shortest paths between the ‘on-ramps’ and ‘off-ramps’. If none of the ‘off-ramps’ is reachable from the ‘on-ramps’, we allow the use of side roads by weighing each side road as high as 10 highway interactions. We include all equidistant shortest paths.

Since the effect of a pharmacogene on the drug action is less likely to stem from transcriptional regulation of the pharmacogene, we ignored transcriptional regulation of the pharmacogene when it was a destination point. The effect of this decision was negligible in line with the robustness of the pathways (Methods section 4 and Figure 2). The result of the construction stage is the union of all the start-destination routes.

In the second, pruning, stage, we applied a conservative approach in which we retained only higher confidence interactions that are traversed by tours from more than one of the five start-destination pair types (e.g. tours that connect a drug target and a pharmacogene and a tour that connects a pharmacogene to a disease gene). If following this pruning, a pharmacogene or disease gene becomes disconnected, we also retain all the shortest routes, discovered in the construction stage, that connect the drug targets to that pharmacogene or disease gene.

4. Building drug-specific PD pathways

Some of the drugs in our set of 113 belong to the same drug family. However, no two drugs shared the exact set of inputs (drug targets, pharmacogenes and disease genes). Specifically, more than 90% of the drugs share less than 80% identity in their associated gene set and more than 90% of the drugs have a non-redundant chemical structure with Tanimoto coefficient lower than 0.7.

In order to simulate construction of de-novo pathways, when building pathways for drugs that have a gold-standard curated PD pathway, we converted the highways that are specific to that PD pathway to side roads.

In order to test the robustness of the drug-specific pathways, we systematically constructed the drug-specific pathway after removal of each of the pharmacogenes or disease genes. We observed high robustness in terms of inferred genes or inferred interactions. The Jaccard scores (47) between the sets of genes or interactions in the drug-specific pathways and the leave-one-out pathways, (computed as the size of the intersections between the two sets divided by the size of the union of the two sets), was 0.87±0.08 and 0.85±0.08, respectively. As anticipated, the greater the number of drug-associated pharmacogenes and disease genes, the more robust is the drug-specific pathway to their removal (Figure 2).

We performed bi-clustering of the drugs by computing the relative overlap of the genes in the inferred pathways of each pair of drugs (Jaccard score), using a modification of the spectral co-clustering algorithm of (48), whereby each bi-cluster is further clustered by single-linkage hierarchical clustering.

For evaluation purposes, we constructed 200 random pathways per drug by connecting the known drug targets with the same number of randomly picked genes from the interaction network as the set of known pharmacogenes and disease genes. We controlled for the randomly selected genes in three ways: (1) random shuffling of the known pharmacogenes and disease genes (100 pathways), (2) maintaining the same network degree distribution as the network degree distribution of the true pharmacogenes and disease genes of that drug (50 random pathways) and (3) maintaining the same distribution of network distances between the random pharmacogenes and disease genes as the network distances between the true ones (50 random pathways). We required the shuffled and randomly picked genes to be distant by at least the mean network path (more than two edges) from the true pharmacogenes and disease genes. The performance difference between each type of random pathways was negligible (1% difference in the number of significant drugs in PubMed test and no difference in the enrichment against curated PD pathways).

For the first support of literature co-occurrences, we queried PubMed abstracts and PubMed central whole articles using the PubMed E-Utilities interface for all pairwise combinations of one of the 113 drugs (using generic names) and one of the 19,176 genes in the gene network (official gene symbols). We manually excluded 63 ambiguous gene names (e.g. a known English word such as ‘rest’ or ‘tag’ or a prevalent abbreviation such as ‘ORF’ or ‘PDF’). Our statistics included rank sum comparison of the number of literature co-occurrences of the drug and the drug-specific inferred genes to (1) co-occurrences of all other genes (p-values for all the drugs below FDR of 0.05), (2) un-associated genes involved in the curated cellular pathways (109/113 drugs below FDR of 0.05), (3) all inferred genes from other pathways, or (4) the inferred genes from the random pathways constructed for that drug. The rank-sum p-values were lower than p-values obtained upon shuffling of gene assignments to drug-specific pathways.

For drug-specific pathways which were not enriched with the gold standard PD pathways, we measured the average distance between the PD pathway and the drug-specific pathway over the network (ignoring the distinction between highways and side-roads), compared to the distance between the PD pathway and the random pathways. As expected, switching roles between highways and side-roads resulted in only four enriched known PD pathway.

5. Applying inferred pathways to suggest drug repositioning

We considered two potential criteria for applying the drug-specific pathways for drug repositioning. A drug repositioning opportunity for drug A is found when: (1) a known target of a drug A appears along the path connecting the known target and disease gene of drug B (drug A may treat the disease that drug B treats); or (2) a known disease gene appears in the pathway computed for drug A (drug A may treat this newly implicated disease) (Figure 4).

Figure 4

Illustration of the drug repositioning prediction scheme. A repositioning of drug A is suggested if: (1) a drug target of an drug A is found along the path between the drug target and disease genes of drug B, or (2) a disease gene unrelated to drug A is found within the pathway of drug A.

For verification of our drug repositioning predictions, we downloaded phases I–III clinical trial up to June 28th, 2013 from the clinicaltrials.gov website. Drug names were matched to DrugBank generic, synonymous and brand names. Condition names were converted to OMIM disease names using the MetaMap tool (49) and filtering operations described in (44). When a drug-disease pair appears in more than one phase, we chose the highest phase. Overall, we obtained 410 unique drug-disease pairs that involve both drugs and diseases in our prediction set.

Off-label drug indications, as well as approved indications that do not appear in our strict gold standard were obtained from electronic health records from Stanford Hospital.

We found additional support comparing to a ‘silver-standard’ set of drug-disease associations. As described in (10), the authors constructed drug-disease associations from four independent sources. Two of the sources, based on extraction from textual indications, were noisier and required additional evidence to be included in the gold standard. However, they used the remaining set of single-evidence associations as a ‘silver standard’ for verification purposes.