PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE/ACM Trans Comput Biol Bioinform. Author manuscript; available in PMC 2017 May 1.
Published in final edited form as:
PMCID: PMC4951188
NIHMSID: NIHMS794410

Efficient Drug-pathway Association Analysis via Integrative Penalized Matrix Decomposition

Abstract

Traditional drug discovery practice usually follows the “one drug – one target” approach, seeking to identify drug molecules that act on individual targets, which ignores the systemic nature of human diseases. Pathway-based drug discovery recently emerged as an appealing approach to overcome this limitation. An important first step of such pathway-based drug discovery is to identify associations between drug molecules and biological pathways. This task has been made feasible by the accumulating data from high-throughput transcription and drug sensitivity profiling. In this paper, we developed “iPaD”, an integrative Penalized Matrix Decomposition method to identify drug-pathway associations through jointly modeling of such high-throughput transcription and drug sensitivity data. A scalable bi-convex optimization algorithm was implemented and gave iPaD tremendous advantage in computational efficiency over current state-of-the-art method, which allows it to handle the ever-growing large-scale data sets that current method cannot afford to. On two widely used real data sets, iPaD also significantly outperformed the current method in terms of the number of validated drug-pathway associations that were identified. The Matlab code of our algorithm publicly available at http://licong-jason.github.io/iPaD/

Index Terms: drug discovery, high dimensional data, matrix decomposition, latent factor model, penalized method, template

I. INTRODUCTION

Identification of potential therapeutic drug molecules for a disease of interest, which is known as “drug-discovery”, is a critical step in new drug development. Traditional drug discovery practice in the pharmaceutical industry usually follows the “one drug – one target” approach, which seeks to identify drugs that act on individual molecules that are involved in the pathological process [1]. However, this approach ignores the systemic nature of human diseases and may account for a significant increase in the attrition rate of new candidate drugs due to low efficacy and serious side effects [2], [3]. Specifically, the pathological process of a disease usually involves extremely complex interplays between many functionally related biological molecules within certain disease-related pathways [4]. Targeting a single component in a disease-related pathway may not provide sufficient interference to the whole pathological process, which may result in unsatisfactory therapeutic efficacy. Moreover, it does not allow researchers to investigate the effects of drug molecules at a systems level, limiting its utility in evaluating drug safety and toxicity in early drug development stages [5].

Pathway-based drug discovery recently emerged as an appealing approach to overcome the limitations of the one drug – one target approach [6]–[8]. Instead of focusing on the interference of drug molecules on individual targets, pathway-based analysis seeks to understand the effect of a drug molecule on a whole biological pathway. The knowledge about the effects of drug molecules upon the disease-related pathways would allow researchers to identify more effective candidate drugs. However, currently not much is known about the complex relationships between drug molecules and biological pathways, and an important first step is to infer associations between them before delving into their complex relationships.

High-throughput transcription and drug sensitivity profiles have become routine products of pharmacogenomics studies [9]–[11]. For example, the NCI-60 project [12], [13] performed whole genome transcription profiling of 60 cancer cell lines and measured their sensitivities to several thousands of drug molecules. Assuming the transcription and drug sensitivity levels are both associated with the activity levels of certain biological pathways, joint modeling of these two types of data could provide valuable information regarding the pathway activities, drug-pathway and gene-pathway relationships. With this mindset, [14] proposed an integrative statistical framework named “iFad”, which is built on a Bayesian sparse factor analysis model, to infer drug-pathway associations by jointly analyzing the gene expression data and drug sensitivity data. They applied their method to the NCI-60 gene expression and drug sensitivity data set and identified significant drug associations for eight types of cancers separately. Despite its satisfactory power in identifying drug-pathway associations, iFad performs statistical inferences using Markov Chain Monte Carlo (MCMC), which is very computationally expensive. Moreover, its results might be sensitive to several prior parameters that need to be specified by the users.

Besides the NCI-60 project, a more recent effort to characterize the transcription and drug sensitivity profiles of various cancer cell lines is the Cancer Cell Line Encyclopedia (CCLE) project [15]. Compared with the NCI-60 data set, a noteworthy feature of the CCLE data set is its large sample size – more than 400 cancer cell lines have been profiled for their transcription levels and drug sensitivities. Analyzing such a large data set using iFad would incur very huge computational cost. In fact, in [14], iFad was only used to analyze each cancer type separately (each cancer type has only 6 ~ 9 samples). Analyzing all the cell lines together would cost tens to hundreds of hours depending on the number of MCMC iterations. Apparently, the ever-increasing data sizes of paired transcription and drug sensitivity data sets present an pressing need for more computationally efficient methods for drug-pathway association analysis. In this article, we introduce “iPaD”, an integrative Penalized Matrix Decomposition method for drug-pathway association analysis. Compared with iFad, our method has two unique features: 1) the association analysis is formulated into a bi-convex optimization problem which can be solved with a very efficient algorithm and hence greatly reduces the computational cost compared to the MCMC algorithm in iFad; 2) our method is almost tuning-free – the only tuning parameter can be easily determined through cross-validation. We demonstrated the practical utility of our method on two relatively large real data sets, the pooled NCI-60 data set with all eight types of cancers and the CCLE data set. We also compared power to detect drug-pathway associations between iPaD and iFad through comprehensive simulations. We organize the rest of this article as follows. We explain the iPaD model and its optimization algorithms in the next section, followed by real data applications and results for simulation studies. We summarize the conclusions in the final section.

II. METHODS

Model description

Consider N samples (usually cell lines) with both transcription and drug sensitivity profiles available. Let Y(1)N×G(1) and Y(2)N×G(2) denote the transcription and drug sensitivity profiles respectively, where G(1) is the number of genes that are measured and G(2) is the number of drug molecules that are assayed. We consider the activity levels of biological pathways as a set of common latent factors that underly the gene expression and drug sensitivity profiles. The activity levels of K pathways for the N cell lines are denoted as a matrix X [set membership] RN×K and we use B(1) and B(2) to denote the loading matrices for pathway-gene associations and pathway-drug associations. Then the gene expression and drug sensitivity profiles are modeled as follows:

Y(1)=XB(1)+E(1),Y(2)=XB(2)+E(2),
(1)

where E(1)N×G(1) and E(2)N×G(2) are residuals. Instead of adopting a Bayesian inference approach, we seek to estimate the latent factors X and the loading matrices B(1) and B(2) by minimizing the following sum of squared residuals,

Y(1)XB(1)F2+Y(2)XB(2)F2,

where ‖·‖F stands for the Frobenius norm, i.e. AF=ijAij2. We introduce an indicator matrix L(1){0,1}K×G(1) to represent the prior knowledge about gene-pathway relationships. Specifically, we let Li,j(1)=1 indicate that the j-th gene is involved in the i-th pathway and Li,j(1)=0 otherwise. Because gene-pathway relationships are relatively well defined in the public databases compared to drug-pathway relationships, we assume that the known gene-pathway relationships are complete and accurate and focus on identifying novel drug-pathway associations. Moreover, B(2) matrix is typically sparse because a drug molecule is usually only associated with a few pathways, and vice versa. Based on these two assumptions, we formulate the following optimization problem:

minimizeX,B(1),B(2)Y(1)XB(1)F2+Y(2)XB(2)F2+λB(2)1subjecttoiXi,j21,j=1,,K; Bi,j(1)=0,(i,j):Li,j(1)=0;
(2)

where ‖·‖1 stands for the [ell]1 norm, i.e. ‖A1 = Σi Σj |Ai j|. This [ell]1 norm penalty is an effective and popular device to achieve sparse solutions in high-dimensional regression which was first introduced in the so-called “lasso” (Least Absolute Shrinkage and Selection Operator) regression [16], which attempts to solve a high-dimensional regression problem while applying an [ell]1 norm constraint on its solution. With such a constraint, some coefficients in the solution may be shrunk to exact zeros. The constraint on B(1) is used to incorporate known gene-pathway relationships. The constraint on the columns of X ensures the identifiability of the model.

Optimization algorithm

We note that problem (2) is a bi-convex optimization problem, i.e., when X is fixed, optimizing B(1) and B(2) are both convex problems and when B(1) and B(2) are fixed, optimizing X is also a convex problem. This property naturally motivates the following alternating optimization algorithm (Algorithm).

Next, we describe the optimization algorithm for X, B(1) and B(2) separately.

Optimizing X

The optimization problem for X can be rewritten as follows:

minimizeXYXBF2subjecttoiXi,j21,j=1,,K
(3)

where Y = [Y(1) Y(2)] and B = [B(1) B(2)]. We used an iterative projected gradient descent algorithm to solve this problem. Specifically, at each iteration, the gradient of the objective function is calculated and we let X take a step along the negative direction of the gradient and projected it to the feasible region: {X:iXi,j21,j=1,,K}. We use the Nesterov’s method to accelerate its convergence [17]. In Nesterov’s method, after each gradient descent iteration, the solution is moved further in the direction and step size given by the the last two steps, achieving a convergence rate of order 1/t2 instead of order 1/t of the ordinary gradient descent method.

Algorithm I: The alternating optimization algorithm for iPaD

Data Input: Y(1), Y(2), L(1)
Parameter: λ
Initialization: Set B(1) = L(1) and set B(2) = 0

Optimization:
 1. Optimize X:
 X=argminXY(1)XB(1)F2+Y(2)XB(2)F2subject toiXi,j21,j=1,,K
 2. Optimize B(1):
 B(1)=argminB(1)Y(1)XB(1)F2subject toBi,j(1)=0,(i,j):Li,j(1)=0
 3. Optimize B(2):
 B(2)=argminB(2)Y(2)XB(2)F2+λB(2)1
 4. Repeat steps 1, 2 and 3 until convergence.

1) Optimizing B(1)

The optimization of B(1) is relatively trivial because we can decompose the problem into G(1) separate ordinary least squares (OLS) problems by observing that each column of B(1) can be optimized separately:

Forg{1,2,,G(1)},minimizeB:,g(1)Y:,g(1)X:,L:,g(1)BL:,g(1),g(1)22
(4)

where Y:,g(1) refers to the g-th column of Y(1), X:,L:,g(1) refers to the sub-matrix of X consisting of the columns corresponding to the non-zero entries of L:,g(1) and BL:,g(1),g(1) refers to a sub-vector of the g-th column vector of B(1) corresponding to the non-zero entries of L:,g(1).

2) Optimizing B(2)

Similar to B(1), we observe that each column of B(2) can also be optimized separately, which decomposes the original problem into G(2) separate lasso regressions:

Forg{1,2,,G(2)},minimizeB:,g(2)Y:,g(2)XB:,g(2)22+λB:,g(2)1
(5)

We use the readily available and very efficient coordinate descent algorithm to solve each lasso problem [18], [19].

We note that this setup does not allow incorporating known drug-pathway associations. When such prior knowledge is available and needs to be incorporated, we modify problem (5) as follows,

Forg{1,2,,G(2)},minimizeB:,g(2)Y:,g(2)XB:,g(2)22+λ(B(1L:,g(2)),g(2)1+BL:,g(2),g(2)2)
(6)

where L(2){0,1}K×G(2) is the indicator matrix representing the known drug-pathway associations similarly to L(1).

Handling missing values

In this section, we describe how missing values in the transcription data and drug sensitivity data are dealt with in iPaD. As mentioned earlier, when optimizing B(1) and B(2), each column of these two matrices can be solved separately. Therefore, missing values can be simply removed from each column of Y(1) and Y(2) while solving each column of B(1) and B(2).

Optimizing X is not straightforward because it needs to be solved as a whole matrix. Suppose that the observed values in Y are indexed by Ω{0,1}N×(G(1)+G(2)). We define an operator PΩ that projects matrix X onto the linear space supported by Ω,

PΩ(X)i,j={Xi,j,ifΩi,j=1.0,ifΩi,j=0.

Then when Y(1) or Y(2) contains missing values, the optimization problem for X can be written as:

minimizeXPΩ(Y)PΩ(XB)F2subjecttoiXi,j21,j=1,2,,K
(7)

Define Ω[perpendicular] = 1 – Ω as the index matrix for the missing values in Y and rewrite the objective function in problem (7):

PΩ(Y)PΩ(XB)F2=PΩ(Y)(XBPΩ(XB))F2=(PΩ(Y)+PΩ(XB))XBF2
(8)

This simple observation motivates the following iterative algorithm to solve problem (7), in which the estimate of X at each iteration is then plugged into PΩ(XB) for optimizing (8) in the next iteration (Algorithm 2). This is essentially the spirit of the soft-impute algorithm from the incomplete matrix learning problem [20].

We prove that Algorithm 2 indeed gives the global optimal solution for problem (7) in the Supplementary Materials.

Algorithm II: Optimizing X with missing values in Y

Data Input: Y, B, Ω
Initialization: Initialize X

Optimization:
 1. SetY=PΩ(Y)+PΩ(XB)
 2. Optimize X:
 X=argminXYXBF2subject toiXi,j21,j=1,2,,K
 3. Repeat steps 1 and 2 until convergence.

Parameter tuning and significance test

In our method, λ is a critical parameter that controls the sparsity of the drug-pathway association matrix, B(2). Similar to many other penalized methods such as lasso, one way to assess the relative importance of the coefficients in B(2) is to solve the problem for a decreasing sequence of λ values and record the order of the coefficients in which they became non-zero – the more important coefficients are supposed to become non-zeros earlier than the less important ones. However, this procedure does not allow evaluation of the significance of the coefficients. Alternatively, we can first find an appropriate value of λ and then perform permutation test to assess the significance of the coefficients.

We use ten-fold cross-validation to find an appropriate λ value. Specifically, the Y(2) matrix is evenly partitioned into n folds as illustrated in Figure 1. In each round of cross-validation, each of the ten folds is masked as missing values while the other entries are treated as observed values. Then we solve the iPaD problem for a sequence of λ values. The residual sum of squares (RSS) of the masked values for the ten folds are evaluated for each λ. The λ that gives the smallest RSS is chosen. Given the λ value, the significance of the coefficients in B(2) can be assessed via permutation test. In each permutation, the rows in Y(2) are shuffled whereas Y(1) is kept unchanged. After the estimates of the B(2) in the permuted data sets are obtained, the p-value of each coefficient in the B(2) matrix is calculated as follows,

Pi,j=t=1T𝟙(|Bi,j(2)(t)||B^i,j(2)|)T
(9)

where T is the total number of permutations, B(2)(t) is the B(2) estimate in the t-th permutation and B^(2) is the B(2) estimate in the original data.

Fig. 1
An illustration of data partitioning in the cross validation procedure. The black entries in each figure represent one fold of the data. In each round of the cross validation, the black entries are masked as missing values whereas the white entries are ...

III. Results and discussion

Real data applications

We applied iPaD to analyze the same NCI-60 data set that was analyzed in [14]. The detailed data pre-processing procedure can be found in their paper. The processed data set consists of 57 cell lines from eight different cancer types with gene expression levels of 1863 genes covering 58 KEGG pathways [21] and drug sensitivity profiles for 101 chemical compounds/drugs as measured by GI50 values, which is a commonly used measure for drug potency and is defined as the concentration required for 50% of maximum cell growth inhibition. The gene-pathway association matrix L(1) has a density of 3.95% whereas the drug-pathway association prior knowledge matrix L(2) has a density of 0.51%. [14] advocated analyzing the eight types of cancer cell lines separately by demonstrating the heterogeneity among the inferred drug-pathway association patterns across different cancer types. Therefore, we first used our method to analyze the eight types of cancer separately. However, we did not find much consistency between our results and the results reported in [14] (results not shown). A possible explanation is that the very limited sample size (only 6 ~ 9 cell lines per cancer type) for each cancer type probably could not provide reliable statistical inference. Then we pooled the data for the eight cancer types and analyzed them using both iFad and iPaD. For iFad, we ran 10,000 iterations with the first 7,000 iterations discarded as burn-in. MCMC samples were recorded every tenth iteration. Both η0 and η1 parameters were set to zero for gene-pathway associations. η1 was set to zero and η0 was set at 0.05 for drug-pathway associations. As for iPaD, we first used ten-fold cross-validation to choose the λ value and then performed 10,000 permutations to evaluate the significance levels of the identified associations. For the purpose of comparison, we also analyzed the data set with a Gene Set Enrichment Analysis (GSEA) method. Specifically, for each drug, we tested for its association with each gene using a simple linear regression and ranked the genes according to their significance levels of association. Then we performed a Kolmogorov-Smirnov test for the enrichment of significant genes in each pathway.

Due to the stochastic nature of the MCMC algorithm in iFad, we repeated the analysis for five times. On average, 123 drug-pathway pairs had posterior probabilities greater than or equal to 0.9. Among them, 25.2 associations were validated in the CancerResource database [22], which was also used in [14] to validate their inferred drug-pathway associations. Whereas for iPaD, although the optimization algorithm itself is deterministic, the permutation test is stochastic. Therefore we also ran iPaD for five times. On average, 277.2 associations had p-values less than or equal to 0.05, with 81 validated in the CancerResource database. If we set the p-value cutoff at 0.005, 59.4 associations were identified, with 20.2 associations validated in the CancerResource database. Among the top 100 associations identified by GSEA, 23 were validated. Apparently, the drug-pathway pairs identified by iPaD are more enriched for validated associations than those identified by the other two methods. To better compare the results from these two methods, we performed the following analysis. First, the list of all drug-pathway pairs were ranked based on their significance in both iFad, iPaD and GSEA results. Then we slid a bin that consists of 20 consecutive drug-pathway pairs from the top to the bottom of the list. For example, the first bin consists of the 20 most significant drug-pathway pairs and the next bin starts from the second most significant drug-pathway pair and runs through the 21st most significant pair, and so on. Then the number of validated drug-pathway pairs in each bin was counted. The results averaged for the five runs for the first 100 bins were shown in Figure 2a. Clearly, among the top drug-pathway pairs identified by the three methods, iPaD had higher enrichment for validated associations. The top 50 drug-pathway pairs for each run were listed in Table S1–S5 and Table S6–S10 for iFad and iPaD respectively. We also list the top 100 drug-pathway associations identified through GSEA in Table S11.

Fig. 2
The enrichment curves of the results from iFad, iPaD and GSEA on real data sets. All the drug-pathway pairs were ranked based on their posterior probabilities in iFad, p-values in iPaD and GSEA or orders of becoming non-zeros (iPad_seq). Then the numbers ...

We also applied both iFad and iPaD to the data set from the CCLE project. This data set contains 480 cancer cell lines with transcription profiles for 18,988 genes and drug sensitivities profiles for 24 chemical compounds/drugs. The drug sensitivities were measured by EC50, IC50, maximum activity (‘Amax’) and area over the dose-response curve (‘activity area’). Like GI50, both EC50 and IC50 are commonly used measures of drug potency. ‘Amax’, instead, is a measure of drug efficacy. In this paper, we chose the ‘activity area’ as our measure for drug sensitivity due to two considerations: 1) it simultaneously captures both drug potency and efficacy; 2) it had much fewer missing values than the other measures in the CCLE data set. We took log-transformation of both the transcription and drug sensitivity profile data before further analysis. Because this is also a cancer-specific data set, the genes in the CCLE data set were intersected with the 1,863 genes used in [14] and only 1,802 genes were left. The same 58 KEGG pathways used in [14] were also used in our analysis. The CCLE data set also provided the known targets of each of the 24 drugs. However, the target genes of two drugs were not annotated in any KEGG pathways and thus we were excluded them from our analysis. For the remaining 22 drugs, a pathway is defined as the target of a drug if any member of the pathway is a target of the drug. The known drug-pathway associations were not used as prior knowledge but as post-analysis validations. Because the CCLE data set is relatively large, running 10,000 iterations for iFad is infeasible in practice. Therefore we ran only 2,000 iterations, with the first 1,000 iterations discarded as burn-in. Correspondingly, only 2,000 permutations were performed in the iPaD analysis. However, the estimated p-values are of limited resolutions with only 2,000 permutations. To overcome this limitation, we estimated the p-values by approximating the null distribution of each coefficient as a mixture of a point mass at zero and a normal distribution.

Again, the analysis was repeated five times for iPaD and iFad. Moreover, we also analyzed the CCLE data set using the GSEA approach described earlier. On average, iFad identified 39.4 drug-pathway pairs that had posterior probabilities no smaller than 0.9, with 4.8 of them being validated associations. Whereas for iPaD, 82 drug-pathway pairs had p-values no greater than 0.05 and 23 pairs were among the validated associations. Again, iPaD outperformed iFad in identifying the validated associations. As for GSEA results, among the top 100 associations, only 1 was validated. As a more comprehensive comparison of the results of the iFad and iPaD, we also performed similar analysis as we did on the NCI-60 results. The enrichment curves for the first 100 bins for the two methods are shown in Figure 2b. The top 10 associations identified by iPaD were listed in Table 1. In the iPaD results, two EGFR inhibitors, Erlotinib and Lapatinib are associated with the ErbB signaling pathway. It is well known that EGFR is a key component of the ErbB signaling pathway and the latter plays a key role in the regulation of cell proliferation, migration, differentiation, apoptosis. Mis-regulations of these activities are hallmarks most cancer cells. In addition, two Raf kinase inhibitors, Sorafenib and PLX4720 are associated with the Chronic myeloid leukemia pathway. In the Chronic myeloid leukemia pathway, Raf kinase is a critical player that regulates the MAPK signaling pathway and hence controls cell proliferation. We also note that Panobinostat, a histone deacetylase (HDAC) inhibitor, was found to be associated with the regulation of actin cytoskeleton pathway. Although HDAC and cytoskeleton are seemingly unrelated, we found that a study published in 2012 [23] suggested that Panobinostat could also induce hyperacetylation of tublin, a major component of the microtubule cytoskeleton, and affects microtubule dynamics and organization. As for the iFad results, there was not much consistency between the top associations from the five repeats. The identified drug-pathway pairs for each run are listed in Table S12–S16 and Table S17–S21 for iFad and iPaD, respectively. We also list the top 100 drug-pathway associations identified through GSEA in Table S22.

Table 1
The top 10 associations identified by iPaD on CCLE data set

We note that although pooling together data from different cancer types allows us to increase the sample size and better identify the pan-cancer common signals, it may also dilute the signals that are specific to certain cancer types. Therefore we also performed an analysis on the lung cancer cell lines of the CCLE data set. We chose lung cancer between it has the largest sample size (N=91) among all the cancer types. Similarly, iFad, iPaD and GSEA were separately used to analyze this data set. Similar to previous results, iPaD achieved the highest enrichment for validated associations, followed by iFad, and then GSEA. The enrichment curves for the three methods are shown in Figure 2c. The 15 associations with p-values smaller than 0.01 in iPaD were listed in Table 2. Similar to previous results, EGFR inhibitors such as Lapatinib and Erlotinib continue to appear in validated association. In addition, multi-kinase inhibitor Sorafenib was found to be associated with the Focal adhesion pathway. Focal adhesions are sub-cellular structures that link cells to the extracellular matrix. Several studies have already reported the up-regulation or activation of focal adhesion kinases in human lung cancer cells [24], [25]. In the focal adhesion pathway, the focal adhesion kinase is regulated by several target kinases of Sorafenib such as KDR, FLT4, PDGFRB and so on. A more comprehensive list of the identified associations on this lung cancer data set is in the supplementary materials (Table S23–S33).

Table 2
The top 15 associations identified by iPaD on the CCLE Lung Cancer data set

The permutation test allows us to evaluate the significance of the identified drug-pathway associations at the expense of substantially increased computational cost. However, we note that such cost is still much more affordable than the MCMC in iFad. For example, on the NCI-60 data set, 10,000 MCMC iterations cost around 230 hours for iFad. Whereas for iPaD, even with 10,000 permutations, it took only about 30 hours. For the CCLE data set, 2,000 MCMC iterations cost about 150 hours. Yet it took only 6 hours on average to run 2,000 permutations with iPaD. All the computation times in this article were estimated on a standard laptop computer (2.4GHz dual core CPU with 4G memory running the Mac OS X 10.9 operating system). In addition, if the permutation test is still not affordable, iPaD could offer an alternative way to evaluate the relative importance of the identified associations without performing the costly permutation test. When a decreasing sequence of λ values are solved, the order in which the drug-pathway association coefficients become non-zero can be used as indications of their relative importances, i.e., the coefficients that becomes non-zeros at larger λ values can be considered as more important. In the real data applications, we also plotted the enrichment curve obtained through this procedure in Figure 2 (denoted as ‘iPaD_seq’). Compared with their rankings of p-values, the drug-pathway pairs may be ranked differently using this procedure. But it is noteworthy that it still achieved better overall enrichment for validated associations than iFad on the two real data sets. Moreover, the computational cost was further greatly reduced – it took only about 15 and 30 minutes to obtain the rankings of the first 50% of the drug-pathway pairs for the NCI-60 data set and the CCLE data set, respectively.

Simulation studies

We first generated a pathway activity level matrix X [set membership] RN×K with each entry simulated from a standard normal distribution independently. To better mimic real applications, we generated an indicator matrix L(1){0,1}K×G(1) by randomly choosing K pathways from the 58 KEGG pathways used in the previous section, where G(1) is the number of unique genes that are involved in any of the K pathways. Given a pre-defined η value that controls the sparsity of B(2), we generated another indicator matrix L(2){0,1}K×G(2) by randomly choosing ηKG(2) entries to be ones. After these two indicator matrices were generated, we simulated B(1)N×G(1) and B(2)N×G(2) by simulating their entries corresponding to the ones in their indicator matrices from independent standard normal distributions while keeping other entries as zeros. Finally, we simulated noise matrices E(1)N×G(1) and E(2)N×G(2) based on a pre-defined signal-to-noise ratio (SNR). Specifically, we first obtained Y(1)=XB(1) and Y(2)=XB(2) and calculated the variance of each column of E(1) and E(2) as VAR(E:,g(1))=VAR(Y:,g(1))/SNR and VAR(E:,g(2))=VAR(Y:,g(2))/SNR. Then each column of the two noise matrices was simulated from a normal distribution with zero mean and its calculated variance. The simulated transcription data and drug sensitivity data were obtained by Y(1)=Y(1)+E(1) and Y(2)=Y(2)+E(2). When a whole column of Y(2) is zero (this will happen when a drug has no target pathway), we use the sample variance of all the entries in the non-zero columns as its variance and simulate its entries from a mean-zero normal distribution.

After the simulation data were generated, we used both iPaD and iFad to analyze the data and measured their performances by the AUC (area under the receiver operating characteristic curve). For both methods, the gene-pathway relationships in L(1) matrix was treated as prior knowledge whereas L(2) was treated as unknown. For iPaD, we solved a decreasing sequence of λ values until 70% of the coefficients in B(2) were non-zero and compared the order in which the coefficients in B(2) became non-zero against their true association status in L(2) to obtain the ROC (receiver operating characteristic) curve and then calculated the area under this curve (AUC). For iFad, both η0 and η1 parameters were set to zero for gene-pathway associations so that their association statuses were fixed during the MCMC iterations. The η1 parameter was set to 0 for drug-pathway associations so that the known associations were never “turned off” whereas the η0 for drug-pathway associations were set to 0.1 to allow novel associations to be inferred. The default values were used for the other parameters. The total number of MCMC iterations was set to 500 (burn-in = 250). MCMC samples were recorded every tenth iteration. The AUCs for iFad were calculated by comparing the association posterior probabilities for the drug-pathway pairs and their true association statuses. We only ran 500 iterations for two reasons: 1) the AUC almost plateaued at 500 iterations; 2) the computational cost is too high. Because the performance of iFad also depends on the number of MCMC iterations, we also evaluated its performance with 10 and 50 iterations. For the purpose of comparison, we also analyzed the simulated data set with GSEA.

The performances in terms of AUC and comprehensive time cost profiles of three methods are shown in Table 3 and and44 (the time cost for GSEA is not shown because its AUC is not competitive). We let the sample size N = 480, the number of drugs G(2) = 22, the number of pathways K = 58, η = 0.1, and varied the SNR from 0.1 to 1. We also evaluated the effect of sample size by fixing SNR = 0.5 and varied N from 120 to 480. As expected, the AUC increases as the sample size and SNR increase. We also simulated data with unequal SNRs for transcription levels and drug sensitivities and with different η values. The AUCs of iFad and iPaD are always similar to each other regardless of the sample size and SNR when we run enough number of iterations for iFad. Both methods substantially outperformed GSEA. Although the performances of the two methods were similar, iPaD cost much less time than iFad. For example, when N = 480, it took about 5 minutes to run iPaD on the simulated data set. Whereas 500 MCMC iterations would cost more than one day for iFad. Even only 10 iterations would cost as much as half an hour for iFad.

Table 3
The simulation settings and time costs of the iFad and iPaD
Table 4
The performances of the three methods on the simulated data sets

We note that although the two methods had similar performances on the simulated data sets, iPaD surprisingly achieved higher enrichment for validated associations in the real data applications. One possible explanation is that the MCMC algorithm in iFad did not converge with the number of MCMC iterations used in our analyses. However, further increasing the number of iterations will make the computation time even more unaffordable in practice. Another explanation is that the default prior parameters of iFad might not be appropriate for the real data sets in our analyses. If this is the case, the results might indicate that iFad is not very robust with respect to the choice of its prior parameters in practice. It could also be explained by the potential difference between the simulation settings and the unknown nature of the real data sets. Finally, it might also be the case that the known drug-pathway associations are not 100% complete and accurate. We admit that these associations are probably not complete, which may limit the interpretability of the results to some extent. But we think their accuracy should be relatively reliable since the drug targets are either well established in the literature or carefully curated in the specific databases by experts. Therefore we argue that our results still provide compelling evidence that iPaD performs better than iFad on the two real data sets.

IV. CONCLUSIONS

Various high-throughput experiment data have become increasingly available to facilitate drug target identification, which have also posed a major challenge for effectively integrating these different types of data. We have developed a method called iPaD, which is based on penalized matrix decomposition, to jointly analyze transcription and drug sensitivity data to identify drug-pathway associations. Thanks to its scalable bi-convex optimization algorithm, iPaD has tremendous advantage in computational efficiency over the existing state-of-the-art method called iFad, which has been demonstrated through both simulation studies and real data applications. Another notable advantage of iPaD is that only one parameter needs to be tuned and we have derived an efficient cross-validation procedure to tune this parameter. Yet iFad has quite a few prior parameters that need to be specified by users. In terms of the accuracy of the identified associations, iFad and iPaD had similar performances on simulated data sets. In real data applications, iPaD seemed to outperform iFad, although the assessment might be complicated by the completeness and accuracy of the known associations that were used for validation. We also note for the two data sets we analyzed, the compounds in the CCLE data set are better characterized than those in the NCI-60 data set in general. Therefore the known associations are more reliable for the CCLE data set and interpretation of its results is easier than that of the NCI-60 data set. Therefore in the real data analysis, we focused on the interpretation of the results from the CCLE data set.

We note that statistical inference/test is notoriously difficult for problems with lasso-type penalty and therefore we had to use a permutation test to evaluate the significance levels of the identified associations. Although the permutation test substantially increased the computational cost, iPaD still takes much less time than iFad in the real data applications. In addition, compared with MCMC, permutation is more amenable to parallel computing, and therefore the computation time of iPaD can be further reduced if the permutation procedure was parallelized. Moreover, as demonstrated in both simulation studies and real data applications, a much more computationally affordable and still useful way to evaluate the relative importance the identified drug-pathway association is to rank them based on the order in which their coefficients become non-zeros. Indeed, this strategy cannot provide rigorous statistical inference as in iFad. But it should be kept in mind that the posterior probabilities provided by iFad need to be interpreted with caution as they can be affected by the several pre-specified prior parameters, especially the prior probability of each drug-pathway association coefficient being non-zero.

From a broader perspective, integration of distinct high-throughput data sets has also become a trending theme in biomedical research. For example, in the cancer genomics community, much effort has been spent on integrating high-throughput genomic data generated from various platforms (e.g. somatic mutations data, gene expression data, epigenomics data, etc). A commonly used strategy is also to treat various types of observed data as manifestations of a set of shared underlying latent factors [26], [27]. As these data sets become larger and more heterogeneous, a major task in the future is to develop useful and computationally efficient methods to perform joint analysis upon them. Compared with a previous method for drug-pathway association analysis, our method has greatly increased the computational efficiency without much, if any, sacrifice in statistical power. However, as more types of high-dimensional data sets become available in the future, their integrative analysis will call for more novel methods to address the forthcoming statistical and computational challenges.

Supplementary Material

tcbb-l-2462344-mm.zip

Acknowledgments

This study was supported by NIH grants GM59507 and CA154295 and NSF grant DMS 1106738. We also thank Yale University High Performance Computing Center (funded by NIH RR19895) for the computation resource and data storage.

Biographies

An external file that holds a picture, illustration, etc.
Object name is nihms794410b1.gif

Cong Li received his bachelor’s degree in Bioinformatics from Huazhong University of Science and Technology, Wuhan, China in 2010. He also just received his Ph.D. degree in Computational Biology and Bioinformatics from Yale University, New Haven, Connecticut in 2015. His current research interests include statistical genetics and genomics, bioinformatics, pharmacogenomics, high-dimensional statistics and data integration.

An external file that holds a picture, illustration, etc.
Object name is nihms794410b2.gif

Can Yang received the bachelor’s and master’s degrees in Automatics from Zhejiang University, China, in 2003 and 2006, respectively. He received the Ph.D. degree in electronic and computer engineering from the Hong Kong University of Science and Technology in 2011. He worked as a postdoctoral associate (2011–2012) and an associate research scientist (2012–2014) at Yale University, New Haven, Connecticut. Now he is working as an assistant professor at department of mathematics, Hong Kong Baptist University. His research interests include statistical learning and bioinformatics.

An external file that holds a picture, illustration, etc.
Object name is nihms794410b3.gif

Greg Hather holds a M.A. in Physics and a Ph.D. in Statistics from the University of California, Berkeley. Since graduating in 2008, he has worked for the healthcare and pharmaceutical industry, primarily supporting drug discovery and manufacturing. His most recent employers include Takeda Pharmaceuticals and Pfizer. His current research interests cover optimal experimental design, empirical Bayes methods, and genomic analysis.

An external file that holds a picture, illustration, etc.
Object name is nihms794410b4.gif

Ray Liu is the head of Statistical Innovation and Consultation group at Takeda Pharmaceutical Company. His group is responsible for providing statistical consultation and project support to various functional areas in pharmaceutical R&D, including Discovery, CMC, Translational Research, and Outcome Research. Ray received his PhD degree from Columbia University. He is a member of ASA, PSTC Stats working group, IQ consortium, and NCBLF. His current research interests include analysis of pharmacogenomics data, semantic analysis and integrated analysis.

An external file that holds a picture, illustration, etc.
Object name is nihms794410b5.gif

Hongyu Zhao is the Ira V. Hiscock Professor of Biostatistics and Professor of Statistics and Genetics, Chair of the Biostatistics Department, and the Co-Director of Graduate Studies of the Inter-Departmental Program in Computational Biology and Bioinformatics at Yale University. He received his B.S. in Probability and Statistics from Peking University in 1990 and Ph.D. in Statistics from the University of California at Berkeley in 1995. His research interests are the developments and applications of statistical methods in molecular biology, genetics, drug developments, and personalized medicine. Some of his recent projects include large scale genome wide studies to identify genetic variants underlying complex diseases (Crohns disease, schizophrenia, bipolar, autism, and substance abuse), biological network modeling and analysis, disease biomarker identifications, proteomics, genome annotations, microbiome analysis, and systems biology study of herbal medicine. He has published over 370 articles on statistics, human genetics, bioinformatics, and proteomics, and edited two books on human genetics analysis and statistical genomics. Dr. Zhao is a Co-Editor of Statistics in Biosciences, and serves on the editorial boards of several leading statistical and genetics journals. He was the recipient of the Mortimer Spiegelman Award for a top statistician in health statistics under the age of 40 awarded by the American Public Health Association. His research has also been recognized by the Evelyn Fix Memorial Medal and Citation by UC Berkeley, a Basil O’Connor Starter Scholar Award by the March of Dimes Foundation, election to the fellowship of the American Association for the Advancement of Science, the American Statistical Association, and the Institute of Mathematical Statistics.

Contributor Information

Cong Li, Program in Computational Biology and Bioinformatics, Yale University, New Haven, CT, USA.

Can Yang, Department of Biostatistics, Yale School of Public Health, New Haven, CT, USA. Department of Mathematics, Hong Kong Baptist University, Hong Kong, China.

Greg Hather, Takeda Pharmaceuticals International Co., Cambridge, MA, USA.

Ray Liu, Takeda Pharmaceuticals International Co., Cambridge, MA, USA.

Hongyu Zhao, Program in Computational Biology and Bioinformatics, Yale University, New Haven, CT, USA. Department of Biostatistics, Yale School of Public Health, New Haven, CT, USA.

References

1. Hopkins AL. Network pharmacology: the next paradigm in drug discovery. Nat Chem Biol. 2008;4(11):682–690. [PubMed]
2. Kola I, Landis J. Can the pharmaceutical industry reduce attrition rates? Nat Rev Drug Discov. 2004;3(8):711–716. [PubMed]
3. Schadt EE, Friend SH, Shaywitz DA. A network view of disease and compound screening. Nat Rev Drug Discov. 2009;8(4):286–295. [PubMed]
4. Lindsay MA. Finding new drug targets in the 21st century. Drug Discov Today. 2005;10(23):1683–1687. [PubMed]
5. Iskar M, Zeller G, Zhao X-M, van Noort V, Bork P. Drug discovery in the age of systems biology: the rise of computational approaches for data integration. Curr Opin Biotechnol. 2012;23(4):609–616. [PubMed]
6. Ganter B, Giroux CN. Emerging applications of network and pathway analysis in drug discovery and development. Curr Opin Drug Discov Devel. 2008;11(1):86–94. [PubMed]
7. Azmi AS. Adopting network pharmacology for cancer drug discovery. Curr Drug Discov Technol. 2013;10(2):95–105. [PubMed]
8. Ma H, Zhao H. Drug target inference through pathway analysis of genomics data. Adv Drug Deliv Rev. 2013;65(7):966–972. [PMC free article] [PubMed]
9. Giuliano KA, Haskins JR, Taylor DL. Advances in high content screening for drug discovery. Assay Drug Dev Technol. 2003;1(4):565–577. [PubMed]
10. Hughes JE. Genomic technologies in drug discovery and development. Drug Discov Today. 1999;4(1):6. [PubMed]
11. Ulrich R, Friend SH. Toxicogenomics and drug discovery: will new technologies help us produce better drugs? Nat Rev Drug Discov. 2002;1(1):84–88. [PubMed]
12. Shoemaker RH. The nci60 human tumour cell line anticancer drug screen. Nat Rev Cancer. 2006;6(10):813–823. [PubMed]
13. Shankavaram UT, Varma S, Kane D, Sunshine M, Chary KK, Reinhold WC, Pommier Y, Weinstein JN. Cellminer: a relational database and query tool for the nci-60 cancer cell lines. BMC Genomics. 2009;10(1):277. [PMC free article] [PubMed]
14. Ma H, Zhao H. ifad: an integrative factor analysis model for drug-pathway association inference. Bioinformatics. 2012;28(14):1911–1918. [PMC free article] [PubMed]
15. Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S, Wilson CJ, Lehár J, Kryukov GV, Sonkin Dea. The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature. 2012;483(7391):603–607. [PMC free article] [PubMed]
16. Tibshirani R. Regression shrinkage and selection via the lasso. J Roy Stat Soc B Met. 1996:267–288.
17. Nesterov Y. A method of solving a convex programming problem with convergence rate o (1/k2),” in. Sov Math Dokl. 1983;27(2):372–376.
18. Friedman J, Hastie T, Höfling H, Tibshirani Rea. Pathwise coordinate optimization. Ann Stat. 2007;1(2):302–332.
19. Tibshirani R, Bien J, Friedman J, Hastie T, Simon N, Taylor J, Tibshirani RJ. Strong rules for discarding predictors in lasso-type problems. J of the Roy Stat Soc B. 2012;74(2):245–266. [PMC free article] [PubMed]
20. Mazumder R, Hastie T, Tibshirani R. Spectral regularization algorithms for learning large incomplete matrices. J of Mach Learn Res. 2010;11:2287–2322. [PMC free article] [PubMed]
21. Kanehisa M, Goto S, Furumichi M, Tanabe M, Hirakawa M. Kegg for representation and analysis of molecular networks involving diseases and drugs. Nucleic Acids Res. 2010;38(suppl 1):D355–D360. [PMC free article] [PubMed]
22. Ahmed J, Meinel T, Dunkel M, Murgueitio MS, Adams R, Blasse C, Eckert A, Preissner S, Preissner R. Cancerresource: a comprehensive database of cancer-relevant proteins and compound interactions supported by experimental knowledge. Nucleic Acids Res. 2011;39(Suppl 1):D960–D967. [PMC free article] [PubMed]
23. Iancu-Rubin C, Gajzer D, Mosoyan G, Feller F, Mascarenhas J, Hoffman R. Panobinostat (lbh589)-induced acetylation of tubulin impairs megakaryocyte maturation and platelet formation. Experimental hematology. 2012;40(7):564–574. [PMC free article] [PubMed]
24. Carelli S, Zadra G, Vaira V, Falleni M, Bottiglieri L, Nosotti M, Di Giulio AM, Gorio A, Bosari S. Up-regulation of focal adhesion kinase in non-small cell lung cancer. Lung Cancer. 2006;53(3):263–271. [PubMed]
25. Mukhopadhyay NK, Gordon GJ, Chen C-J, Bueno R, Sugarbaker DJ, Jaklitsch MT. Activation of focal adhesion kinase in human lung cancer cells involves multiple and potentially parallel signaling events. Journal of cellular and molecular medicine. 2005;9(2):387–397. [PubMed]
26. Shen R, Olshen AB, Ladanyi M. Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis. Bioinformatics. 2009;25(22):2906–2912. [PMC free article] [PubMed]
27. Mo Q, Wang S, Seshan VE, Olshen AB, Schultz N, Sander C, Powers RS, Ladanyi M, Shen R. Pattern discovery and cancer gene identification in integrated cancer genomic data. Proc Natl Acad Sci U S A. 2013;110(11):4245–4250. [PubMed]