Inferring drug-disease associations from integration of chemical, genomic and phenotype data using network propagation

Construct phenotype homo-network

A node in a phenotype homo-network as a disease phenotype is extracted from Online Mendelian Inheritance in Man (OMIM) database [18]. We use a scoring schema of phenotypic similarities as edges that quantitatively measures the phenotypic overlap of OMIM records constructed by van Driel [19] using text mining techniques. If the similarity score of two diseases falls in the range [0.6, 1], it means informative similarity which indicates potentially relevant phenotypic similarity. On the other hand, if the similarity score falls within [0, 0.3], it indicates non-informative similarity. Therefore, we apply a logistic function from [20] to convert the phenotypic similarity scores among diseases into a value as close to 1 as possible while the non-informative score into a value as close to 0 possible over all the entries in the phenotype similarity score matrix. The symmetric similarity matrix W _p (p _i ,p _j ) in phenotype homo-network denotes the phenotypic similarity score between phenotypes p _i and p _j .

Construct drug homo-network using chemical similarity

We extract the FDA-approved drugs and their canonical simplified molecular input line entry specification (SMILES) from DrugBank database [21][22]. We calculate the hashed fingerprints using Chemical Development Kit (CDK) [23]. The chemical similarities are calculated by two hashed fingerprints using Tanimoto coefficient [24]. It calculates the size of the common substructures over the union between two fingerprints of the drugs which is defined as sim(x, x')=|x∩x'|/|x∪x'| between two chemical structures of drug x and x'. The symmetric chemical structures similarity matrix as the edge weight in drug homo-network is denoted as W _d and each value falls in the range between zero (no bits in common) to unity (all bits the same)

Construct protein interaction homo-network using gene expression data

where W _g (g _i ,g _j ) denotes the weight function from gene g _i to gene g _j . $\left|R_{g_i{g}_j}^d\right|$ denotes the absolute value of Pearson correlation coefficient of the interaction between gene g _i and g _j from case data. $E_{g_i}^d$ and $E_{g_i}^n$ are the average gene expression values of gene g _i in case sample and control sample.

Integrated disease, protein interactions and chemical homo-networks

The gene-phenotype hetero-network shows the relationships between disease phenotype and disease-associated genes extracted from OMIM database [18]. The drug-protein hetero-network denotes the drug and its targets which is obtained from DrugBank database [21]. The asymmetric matrices W _pg , W _dg represent the adjacency matrices of link structures from phenotype-gene relationships and drug-target protein interactions, respectively. If drug d _i has a target g _j , then W _dg (d _i , g _j ) = 1, otherwise W _dg (d _i , g _j ) = 0. When a drug target or disease-associated gene has no link with other proteins in PIN, we set the probability of connection to any other protein as 1/(n-1), where n is the total number of proteins in PIN. Since n is usually very large, so the probability will be very small. The reason that we use small probability instead of zero probability is to prevent a node in the network becoming a "sink node" in PIN and allows the probability to be propagated through the node.

Network propagation in the integrated network

We identify the inferring drug-disease associations problem as probability propagation over a network which simulates a random walker stochastically move on query phenotype to its immediate neighbors in heterogeneous network [26, 27]. We adopt the idea from [28] which developed a label propagation algorithm for an integrated network.

where i denotes the node in a homo-network and j denotes the node in the other homo-network.

The first term denotes the random walker can "restarts" to the initial probability distribution among the nodes with the diffusion parameter 1-α. The second term denotes an iterative walk to reach the further nodes in the network based on the transition matrix S and a diffusion parameter α. The random walker will be trapped at initial nodes if α is zero. Let P ^t be a probability distribution where a node in the network holds the probability of finding itself in the iterative random walker process up to the step t. After certain steps, the probabilities will reach a steady state which the difference between P ^t and P ^t-1 measured by L2 norm falls below a very small value such as 10^-9.

Thus, network propagation is calculated with an enriched initialization from the other homo-networks through hetero-subnetworks and the proof of convergence is in [31].

Given a query phenotype, we first set the initial probability distribution over nodes where the probabilities to the query disease nodes set to one and other nodes in the other homo-networks to zero. Second, we apply our network propagation method iteratively until the probability converges on each homo-network. Finally, we use the coverged probabilities of the nodes in homo-networks as initail probability distribution and then repeated the network propagation processes until all homo-networks converge to a final probability distribution.

Evaluation of association specificity between drug and disease

Here, P _i (v) denote the probability of node v in homo-network i using genomic data calculated by our method. The functions avg and std denote the average and standard deviation for the set of reference probabilities P _i ^ref in homo-network i respectively.

Gene expression profile

We adopt microarray data taken from [33] that consists of 62 primary prostate tumors and 41 normal tissues from Stanford Microarray Database (SMD) [34]. We use genome-wide gene expression profiles from tissue samples of 18 healthy normal controls and 27 patients with colorectal cancer evaluated by HG-U133 Plus 2.0 platform microarrays (Affymetrix, Santa Clara) through from Gene Expression Omnibus (GEO) database (GSE4183 and GSE4107) [35, 36].

Protein interaction network

We successfully obtained 137,037 interactions among 13,388 genes by integrating five protein interaction network databases (HPRD, BIND, IntAct, MINT, and OPHID) and by mapping the UniProt protein ID to the human Entrez gene ID, erase the duplicated interaction pairs.

Phenotype network and Phenotype-genotype hetero-network

The OMIM database constructed the catalogue of genetic diseases in human and provides the phenotype-genotype association for 14,433 genes and 5,080 diseases [18]. The gene-phenotype hetero-network contains 275 disease phenotypes and 649 genes from 877 relations while mapping the genes in microarray data and PIN.

Drug network and drug-target hetero-network

We collect 1,571 FDA-approved drugs and 1,410 of them with available SMILES data in DrugBank database [21]. There are 4,456 relations between 1,215 drugs and 1,141 targets to be the drug-target hetero-network.

Benchmark of drug-disease associations

We extract 53 and 106 known associations between drug, prostate cancer and colorectal cancer extracted from CTD database [6] in May 2013 as our benchmark.

The performance of our method

We compare our method with the previous state-of-the-art Cmap project to evaluate the performance [10]. We first prepare the ranked lists of over-expressed genes and down-expressed genes in both prostate and colorectal cancers which are conducted with affymetrix HG-U133A platform for Cmap project. The number of over-expressed and under-expressed gene signature in prostate cancer and colorectal cancer microarray data with different fold changes are shown in Table 1. Given the disease-specific gene signature as input query, the drug compounds ranked lists obtained from Cmap which are scored based on how well they are correlated with the input query. The score of 1 represents the input query perfect matches the changes among drug expression profiles in the database, -1 represents perfect anti-correlation between them, and 0 represents a null match. Both negative and positive enrichment belong to the drug-disease associations. Therefore, we take the absolute values of scores and re-rank them and higher score represents the stronger drug-disease associations. Due to Cmap project is derived from the expression profile of a drug compound to the isolated cell lines, multiple instances may correspond to same drug and even with same dose in the obtained results. We calculate the maximum score of multiple instances that correspond to the same drug as the associated score of a specific drug for dealing with such cases. In prostate and colorectal cancer data, we use 601 drug compounds which are overlapped between DrugBank and Cmap results and the drug rank list have shown in the previous works [37]. We compute observed the area under the receiver operator characteristic (ROC) curve (AUC) for analyzing the quality of performance in Figures 2 and 3, respectively. The Cmap method obtains much lower AUC of 0.59 ± 0.04 and 0.59 ± 0.03 under the gene signature with different fold changes while our method obtains 0.94 and 0.89 in prostate cancer and colorectal cancer respectively. Previous study showed that only depending on the drug response expression profile is incapable to acquire the drug-disease associations accurately due to the profiles generated under different conditions [38]. There are several problems that limit the performance in Cmap project: First, the set of differentially expressed genes that constitute disease signatures or drug signatures were chosen empirically that cannot guarantee the biological relevance of the selected signatures. A bad selection of signatures may tend to capture similarities in the experimental settings rather than revealing the underlying mechanisms. Second, only the overlap among genes between disease state and drug treatments are not quantified as the overall effect of a drug and the values of differential expression are also not taken into account in Cmap project. Therefore, the feasibility and benefits of using network-based analyses of chemical, genomic and phenotype data provides a good chance to reveal drug-disease associations.

Cancer	Prostate cancer					Colorectal cancer
Fold change	1.0	1.1	1.2	1.3	1.4	1.2	1.4	1.6	1.8	2
# over-expressed genes	89	60	45	31	27	539	308	175	106	68
# down-expressed genes	232	156	115	86	58	72	42	27	12	5

The AUC with varying diffusion parameters

We investigate the systematic effect of the different diffusion parameters among networks in our approach. We apply the diffusion parameters α_i $\in${0.1, 0.3, 0.5, 0.7, 0.9} used in the experiments and 125 combinations of the parameters are tested. The ranking performances of our method with different combinations are measured by AUC values in Figure 4. The results show that the diffusion parameters set in drug, gene, and phenotype homo-network α_d, α_g, and α_p as 0.1, 0.7 and 0.3 can get a higher AUC in both prostate and colorectal cancer. A higher similarity score between the chemical structures of a pair of drugs may sometimes have different targets to affect different functional modules in PIN. The drugs with similar chemical structures without binding to similar enzymes would reduce the predictive accuracy for drug similarity [39]. On the other hand, the bit-comparison method of chemical structures using Tanimoto coefficient also has its limitations [40]: First, this kind of method does not include biochemical information at the atomic level of the representation. Second, sometimes it may yield low similarity values and it has an inherent bias related to the size of the molecules that are being sought. Furthermore, many physiological effects cannot be predicted accurately by chemical structure properties alone without more detailed information of metabolic and pharmacokinetic transformations of drugs [41]. Those may be the reason that why the diffusion parameter value in the drug network should be relatively smaller in our experiment. PIN tends to have much higher diffusion parameters to get higher AUC and it is reasonable to infer drug-disease associations not only by targeting the specific proteins directly but also by modulating the pathways involved in the pathological process [42].

The performance of our method with different data source

The drug is designed to target on the primary targets, but it also may interact with unexpected proteins (called off-targets) to trigger the associated biological processes and pathways. It may display undesired off-target toxicity or drug reposition for the disease phenotype. Therefore, we aim to add additional drug targets information to better understanding of dynamic processes under drug exposure. The database STITCH contains the relations between chemicals as chemical-chemical associations (CCA), and chemicals and their interacted proteins as drug-protein interactions (DPI) which are all curated by the evidence derived from experiments, publicly databases and the literature extraction [43]. However, the textual co-occurrence from text mining does not necessarily indicate meaningful relationships [44]. Therefore, we take 29,275 chemical-chemical interactions and 45,567 chemical-protein pairs from STITCH excluded the relations extracted from the text mining method. In the rank lists of 1,410 drugs, our method with the chemical structures similarity from CDK in drug homo-network and the chemical-protein interactions from STITCH in drug-protein hetero-network obtain highest AUC of 0.936 and 0.861 in prostate cancer and colorectal cancer dataset in Figure 5 and 6, respectively. It shows that our method with more potential off-targets from STITCH gets higher AUC than that with only primary drug targets from DrugBank in both prostate cancer and colorectal cancer data. Due to only a limited number of chemical-chemical associations have been identified in STITCH [45], its performance is worse than the associations constructed by chemical structures similarity from CDK in drug homo-network.

Case study: prostate cancer

Case study: colorectal cancer