A gene expression signature distinguishes innate response and resistance to proteasome inhibitors in multiple myeloma

Extensive interindividual variation in response to chemotherapy is a major stumbling block in achieving desirable efficacy in the treatment of cancers, including multiple myeloma (MM). In this study, our goal was to develop a gene expression signature that predicts response specific to proteasome inhibitor (PI) treatment in MM. Using a well-characterized panel of human myeloma cell lines (HMCLs) representing the biological and genetic heterogeneity of MM, we created an in vitro chemosensitivity profile in response to treatment with the four PIs bortezomib, carfilzomib, ixazomib and oprozomib as single agents. Gene expression profiling was performed using next-generation high-throughput RNA-sequencing. Applying machine learning-based computational approaches including the supervised ensemble learning methods Random forest and Random survival forest, we identified a 42-gene expression signature that could not only distinguish good and poor PI response in the HMCL panel, but could also be successfully applied to four different clinical data sets on MM patients undergoing PI-based chemotherapy to distinguish between extraordinary (good and poor) outcomes. Our results demonstrate the use of in vitro modeling and machine learning-based approaches to establish predictive biomarkers of response and resistance to drugs that may serve to better direct myeloma patient treatment options.


INTRODUCTION
Wide interindividual variation in response to chemotherapy is a major limitation in achieving consistent therapeutic effect in many cancers, including multiple myeloma (MM), the second-most common hematologic malignancy with an estimated 30 330 new cases (~2% of all new cancer cases) and 12 650 estimated deaths in 2016 (NCI-SEER (The Surveillance, Epidemiology, and End Results program of the National Cancer Institute) Cancer statistics). [1][2][3][4] Such heterogeneity in response to treatment is governed in large part by the underlying molecular characteristics of the tumor, including differences in the expression of genes involved in mechanisms of chemoresistance. [5][6][7] Deciphering key changes in gene expression levels underlying personalized sensitivity to chemotherapy is therefore essential to predict the efficacy of anticancer drugs and to prevent delay in the selection of more effective alternative strategies Proteasome inhibitors (PIs) are effective chemotherapeutic agents in the treatment of MM, used alone or in combination with other anticancer agents like alkylating agents, corticosteroids, immunomodulatory agents and histone deacetylase inhibitors. [7][8][9] Bortezomib (Bz/Btz/Velcade) was the first PI to be approved by the US Food and Drug Administration for clinical application in 2003 for the treatment of relapsed and refractory MM. [10][11][12] Other examples include second-generation PIs Carfilzomib (Cz/Cfz/ Kyprolis), Oprozomib (Opz) and Ixazomib (Ix/MLN9708/ Ninlaro). 7,10,13 However, MM still remains mostly an incurable disease with 5-year survival rate of 48.5% (NCI-SEER (The Surveillance, Epidemiology, and End Results program of the National Cancer Institute) Cancer statistics). Furthermore, most MM patients ultimately undergo relapse, including patients with good response to initial treatment who eventually develop resistance to the therapy. 7 Moreover, there are reports that patients who fail to respond to Bz may still respond to other PIs. 8,14 Most patients receive PIs in combination with other therapeutic agents; thus, the variability in PI response is difficult to assess. In addition, survival end points in clinical applications are measured in months to years, and hence developing prediction algorithms of response can be a long process.
Therefore, we utilized a collection of more than 50 human myeloma cell lines (HMCLs) generated through the immortalization of primary MM cells that represent a broad spectrum of the biological and genetic heterogeneity of MM 15 to create an in vitro chemosensitivity profile in response to treatment with the four PIs: Bz, Cz, Ix and Opz as single agents. Then, we used machine learning-based computational approaches to identify gene signatures that could distinguish sensitive and resistant responses in cell lines. When applied to gene expression profiling (GEP) data of MM patients from four different PI-based clinical trials, our GEP model of response/resistance to PIs successfully distinguished differences in disease progression and distinguished extraordinary (good and poor) responses. Thus, these results can provide a PI treatment-specific predictor of clinically relevant outcomes that could affect therapeutic choices.

Cell lines
Fifty HMCLs were procured from various institutions, established and characterized, and maintained in HMCL media with interleukin-6. 16 Supplementary Table S1 provides the cytogenetic characteristics of the HMCLs.
In vitro chemosensitivity assays Cell cytotoxicity assays were performed on the HMCLs to create a drug sensitivity profile in response to treatment with increasing concentrations of Bz, Cz, Ix and Opz as single agents and half-maximal inhibitory concentration (IC 50 ) values and area under the survival curve (AUSC) were calculated as described earlier. 16 Briefly, cells were counted using Countess automated cell counter (Invitrogen, Carlsbad, CA, USA) and seeded in 96well plates at a concentration of 4 × 10 5 cells per ml. After 24 h, cells were treated with increasing concentrations of Bz, Cz, Ix and Opz as single agents. Cell viability assays were performed 48 h post treatment using CellTiter-Glo luminescent cell viability assay (Promega, Madison, WI, USA) according to the manufacturer's instructions using Synergy 2 Microplate Reader (BioTek, Winooski, VT, USA) to generate survival curves. Percent survival values were normalized to untreated controls and IC 50 values were determined by calculating the nonlinear regression using sigmoidal doseresponse equation (variable slope). AUSC was calculated using trapezoidal rule and log 2 -transformed for further statistical analysis. 16 Caspase-3/7 activity was evaluated using Caspase-Glo 3/7 Assay kit (Promega) on Synergy 2 Microplate Reader.
Gene expression profiling RNA was isolated from six most Ix-sensitive and six most Ix-resistant cell lines and RNA sequencing was performed on llumina HiSeq 2000 (Illumina, San Diego, CA, USA) using 50 bp paired-end protocol with depth of 420 million reads per sample. RNA-sequencing (RNA-seq) data from CD138selected plasma cells were generated from MM patients enrolled in an ongoing phase-2 Ixazomib clinical trial at Mayo Clinic (Mayo-Ix; NCT01415882) that enrolled patients with relapsed myeloma who had less than six cycles of prior treatment with a Bz-based regimen and were not refractory to Bz 17 (see Supplementary Methods for details).
Gene expression, treatment arm and outcome data on newly diagnosed myeloma patients enrolled in HOVON65/GMMG-HD4 trial (ISRCTN64455289; n = 290) were downloaded from Gene expression omnibus (GEO) (GSE19784). 18 The APEX data set (GSE9782; n = 264) consists of bortezomib-based phase-2 and phase-3 relapsed and/or refractory myeloma clinical trials (The APEX phase-3 trial (039), a companion study (040), the SUMMIT (025) and CREST phase-2 trials (024)). 19 For both these trials, Affymetrix, Santa Clara, CA, USA gene probe set analysis data (U133-Plus2.0 for HOVON65/GMMG-HD4 and HG-U133A/B for APEX) were available. Pretreatment RNA-seq data, treatment arm and clinical outcome information on CoMMpass (Relating Clinical Outcomes in MM to Personal Assessment of Genetic Profile) trial patients were downloaded from the MMRF (Multiple Myeloma Research Foundation) researcher gateway portal (IA7c release; https://research.themmrf.org). The CoMMpass Trial (NCT0145429), sponsored by MMRF, is a non-registrational, longitudinal study of 1000 newly diagnosed MM patients followed over the course of their disease, up to 8 years. 20 Bioinformatics and statistical analysis All statistical analyses were performed using R software environment (https://www.r-project.org/) version 3.3.1 for statistical computing and graphics, and GraphPad Prism v7.0 (GraphPad Software, Inc., La Jolla, CA, USA). Spearman's rank-order and Pearson's product-moment correlation analyses were performed to compare the PI responses. All tests were two sided and P o0.05 was considered statistically significant. Gene expression data were preprocessed, log 2 -transformed and analyzed using Galaxy and Partek Genomics Suite v6.6 (Partek Inc., St Louis, MO, USA) to perform differential expression testing to identify gene expression signatures PI response (details in Supplementary Methods). Analysis of variance model or two-sided paired sample t-test was used to evaluate whether each gene is differentially expressed. Heatmaps were generated using unsupervised hierarchical clustering analysis based on the differentially expressed genes.
Supervised machine learning approaches construct algorithms that learn from training data, build models based on properties of training inputs and thus make learned predictions/decisions on new/test samples. 16,21 Random forest, a supervised ensemble machine learning algorithm, was used to establish the top differentially expressed genes as predictive GEP signatures of PI response. 21 GEP data on top Ix-sensitive vs top Ix-resistant HMCLs (n = 12) were used as 'training data set' to build random forest classification models (decision trees) and predict PI resistance of HMCLs (n = 44) in the 'test data set', the mRNA-seq data obtained from the Keatslab repository (http://www.keatslab.org/data-repository). The average bootstrap prediction error was generated using repeated bootstrapping of the train data set with a k-fold cross-validation (k = 100). The crossvalidation error rate was used to evaluate the accuracy of the method. 16,21 The predictive GEP signature was then applied to the four different PIbased MM clinical trials and random survival forest estimation method for right censored data (randomForestSRC), another supervised machine learning decision-tree based algorithm, was used to predict probability of progression/event (0 = censored; 1 = progression) within the first 3 years for each myeloma patient. 22 The predicted probability values derived using machine learning approaches (random forest and randomForestSRC) were rank-ordered and the predictions for the top (Q3) and bottom quantiles (Q1) were compared with observed PI response using Somers' D xy rank correlation 23 (see details in Supplementary Methods).
Unsupervised K-means clustering was performed using the algorithm of Hartigan and Wong 24 to identify clinically important K-subgroups based on our PI response GEP signature such that N/K~30 (N = total number of subjects in data set). Kaplan-Meier curves for survival were generated for the extraordinary PI response (good vs poor) clusters by computing progression-free survival (PFS) over time. 25 The Kaplan-Meier survival curves were compared statistically using log-rank test and Cox proportion hazard test. 26 Clusters with no 10 were combined for PFS comparisons.
Odds ratios (ORs) between observed clinical responses (available for APEX and CoMMpass data sets) vs extraordinary PI-response K-means clusters were computed using binomial logistic regression analysis (logit model). 27

Ingenuity pathway analysis
The differentially expressed genes were analyzed using Ingenuity Pathway Analysis (IPA) to identify canonical pathways, downstream effects, upstream regulators and causal networks and to perform predictive toxicology analysis using toxicogenomics approaches (IPA-Tox). 28

RESULTS
Wide variability in response to PI treatment To model and characterize the interindividual variation in therapeutic response, we assessed the in vitro PI sensitivity of a large panel of HMCLs. Results of our cytotoxicity assays in 50 HMCLs showed a wide range of response to treatment with the four PIs (Bz, Cz, Ix and Opz) identifying some lines that are highly sensitive and some lines relatively refractory to PIs, as represented by the IC 50 and AUSC plots (Figures 1a and b). In some cases, the AUSC provides a comparative measure of response through the higher test doses, especially if an IC 50 value was not achieved. The summary table (Table 1) depicts the extensive variation in response to each of the PIs among the panel of HMCLs. Our caspase-3/7 activity data corroborated with cytotoxicity data (data not shown). Figure 1c shows a correlation matrix representing comparative analysis of IC 50 values of four PIs across the HMCLs. Statistically significant (adjusted P-values of o0.001; Holm's method) positive correlation was observed between IC 50 values across the four PIs. Thus, in general, sensitivity was highly In (c), the scatterplot matrix is shown as a pairwise correlation of the natural log (Ln) of IC 50 and AUSC values for the response to each PI drug. Scatterplot matrix was generated using the R graphing package ggplot2. correlated among all four PIs, especially among the cell lines representing extraordinarily/exceptionally high and low responses. Therefore, we evaluated bortezomib and ixazomib as model proteasome inhibitor drugs for the gene expression signature-based prediction of extraordinary response in PIinhibitor-based therapies.
However, we did find some examples of HMCLs that were highly sensitive to one PI, and comparably less sensitive to other PIs. This demonstrates tumor heterogeneity even in response to inhibitors of the same class, and further demonstrates some tumors refractory to one PI may still respond to another.
Deriving a GEP-based signature profile of PI response We performed differential gene expression analysis between 5 (top 10%) most Ix-sensitive and 5 (bottom 10%) most Ix-resistant cell lines. Notably, these cell lines showed the same relative sensitivity and resistance to all four proteasome inhibitors. RNAseq data were prefiltered, normalized and used for further analysis (see Supplementary Methods). Differential gene expression analysis was performed using the remaining 7682 genes to identify GEP signatures that distinguish the highly sensitive from the highly resistant HMCLs. Results showed 506 genes differed significantly between the sensitive and the resistant groups (P o 0.05; fold difference ≠1). In all, 141 genes showed |fold difference| 42 and P o 0.05, whereas 42 genes, listed in Table 2, had P o 0.01 (|fold difference| 42) (Figure 2). Subsequent analyses used the more stringent highly variable/differentially expressed 42-gene list. The purpose of this initial selection was to generate a list of genes that were highly variable/differentially expressed between PI-sensitive and PI-resistant cell lines and could then be used in further steps to derive the prediction model and for validation.
GEP signature of PI response is predictive of in vitro PI chemosensitivity in HMCLs and progression in MM clinical trials Our approach was to develop a novel methodology to predict extraordinary/exceptional (very high and very low) PI responses in patients using a myeloma cell line-based model system. The training data set comprised gene expression profiles of HMCLs that represent extremities of PI responses, whereas the initial validation/test data set was derived from an independent gene expression study on HMCLs (n = 44) performed in the Keats laboratory at TGen laboratories (http://www.keatslab.org/datarepository). The final validation data set of human myeloma patients comprised four independent clinical trials (HOVON65/ GMMG-HD4 (n = 290); APEX (n = 264); CoMMpass (n = 765); and the Mayo-Ix trial (n = 22)).
First, the random forest algorithm for classification was 'trained' using mRNA-seq data on the 42-gene signature from 6 most PIsensitive and 6 most PI-resistant cell lines. The out-of-bag estimate, which is a measure of mean prediction error of the training model calculated by subsampling the training data set, was computed as 0%. Concurrently, the classification error rate derived from confusion matrix was also 0, thus validating the   Figure S1 provides the gene importance plot representing the most important genes within the 42-gene expression signature influencing PI response based on variable importance (VIMP) measure derived from random forest analysis. The probability scores of PI resistance were then calculated for each of the HMCLs obtained from the Keatslab data repository ('test' data set). The predicted probabilities of PI resistance were then rank-ordered and Somers' QUOTE rank correlation analysis was performed between the top quantile (Q3) and bottom quantile (Q1) resistance probability values and observed PI chemosensitivity as a binary outcome (sensitive = 0 vs resistance = 1). Results revealed high positive Somers QUOTE rank correlation for the PI drug cytotoxicity parameters (IC 50 and AUSC) ( Although the training set showed good prediction capabilities in an independent set of HMCLs, we were particularly interested to determine whether it is able to stratify patient outcomes in multiple clinical trials. Four different PI-treating trials were examined (treatment details provided in Materials and Methods). Microarray gene expression data from clinical trials was meancentered and scaled before analysis. The standardized transcriptomic profiling data from APEX trial was used as a training data set to perform random survival forest analysis to predict the percent probability of a progression/event within 3 years in the HOVON-GMMG-HD4, CoMMpass and Mayo-Ix clinical trials (test data sets). Somers' QUOTE rank correlation analysis between the top and bottom quantiles of predicted percentage values from random survival forest model on test data and the progression index of the test data sets showed consistent positive values for the HOVON-  (Table 4). Somers' c for the training data set (APEX) was = 0.852. As we were particularly interested in our signature's drug response performance, we chose to look at progression as one measure of response/nonresponse as well as associations that distinguish clinical definition of response, complete response versus nonresponse. 29 K-means clustering was performed on 188 MM patients from the Bz treatment arm of the APEX data set 19 to partition the samples into clusters/subgroups based on the expression of the 42 genes comprising the GEP signature of PI response. Results show significant differences in PFS between the signature-derived poor PI-response and good PI-response groups (hazard ratio (HR) = 2.346; P = 0.0076; Figure 3a). Conversely, no difference in PFS is observed between the K-means clusters when the 42-gene model was applied to the expression data of 76 patients in the dexamethasone arm of APEX phase-3 trial that compared single-agent Bz with high-dose dexamethasone (HR = 1.1; P = 0.732; Figure 3a), showing our GEP signature is drug Figure 2. Heatmap representing differential gene expression between PI-sensitive and PI-resistant myeloma cell lines. Gene expression was z-score normalized (standardized: shifted to mean of 0 and scaled to s.d. of 1) and compared among the five most Ix-responsive and five least Ix-responsive cell lines. Heatmap was generated using the top 42 differentially expressed genes (|fold difference| 42; P o0.01). Columns are ordered by Ix IC 50 of cell; genes are ordered by fold difference. Color indicates fold change between Ix-resistant and Ix-sensitive cell lines.  specific. Concurrently, statistically significant association is observed between the K-means clusters and clinical response in the Bz arm (OR responder vs nonresponder = 5.813; 95% confidence interval = 1.833-20.007; P OR = 0.0036) but not in the Dex arm (OR responder vs nonresponder = 2.139; 95% confidence interval = 0.753-6.326; P OR = 0.158) of the APEX data set.
When applied to the gene expression data from HOVON-GMMG-HD4 clinical trial (n = 290) that implemented a Bz-based drug regimen, 18 the 42-gene signature shows statistically significant differences in PFS among K-means clusters representing good vs poor PI response in the Bz-treated PAD (bortezomib, doxorubicin and dexamethasone) arm (HR = 2.161; P = 0.024), whereas no difference in PFS is observed in the VAD (vincristine, doxorubicin and dexamethasone) arm (HR = 1.282; P = 0.437; Figure 3b). What is particularly striking is that within the first 1.5 years the predicted good response group in the PI-containing therapy had no progression events, whereas 35% of the predicted poor response group had progression events.
Clinical data were available on 765 patients from the CoMMpass study. A total of 253 patients used Bz as first-line therapy, alone or in combinations. Of these, 128 had RNA-seq and PFS data available. K-means clustering based on the expression of the 42 genes distinguishes between good and poor PFS subgroups of myeloma patients (HR = 2.556; P = 0.0277; Figure 3c). In addition, the K-means clusters are also found associated with clinical response (OR = 5.20; 95% confidence interval = 1.22 36.077; P OR = 0.0453). In contrast, among the 96 patients using Lenalidomide (an immuno-modulatory drug) as first-line therapy, no association was observed between the good and poor K-means clusters and PFS (P = 0.49; Figure 3c). Our results thus further validate the PI response specificity of our method.
Finally, we applied our algorithm based on the our 42-gene PI response signature to an ongoing clinical trial at Mayo Clinic that uses Ix-containing drug regimen. Interestingly, even with small numbers of patients available, our GEP-based PI response classifier distinguishes between the top (medium PFS = 22.42 months) and bottom responders (medium PFS = 12 155 months) in this Ix trial (risk ratio = 2.5) (Figure 3d). Notably, all 10 patients in the good performance group were alive at the latest point of the interim analysis, whereas 4 out of 9 patients in the poor PI response cluster died (inset of Figure 3d).

Ingenuity pathway analysis
Genes in the 42-gene GEP model were then analyzed using IPA to assign them to different functional networks. Out of the 42 genes, 12 overlapped with the top IPA network obtained from the analysis of direct and indirect relationships (score = 26), as represented in Supplementary Figure S2a. Genes/molecules in this network include the 26s proteasome complex. Furthermore, when we used the IPA upstream regulator analysis, TFEB and CCND1 were identified as the top hits/factors that may control the genes and pathways highlighted by network analysis (Supplementary Figure S2b). The functional significance is discussed further in the Discussion.

DISCUSSION
Resistance to chemotherapy in cancers is a major limitation in achieving complete and sustained therapeutic effect, and may in fact lead to unwanted exposure to ineffective antitumor agents, thereby increasing the risk of negative side effects. 1 MM is a mostly incurable disease with a median overall survival of 4.6 years despite significant improvement in treatment strategies over the past 12 years. [2][3][4]30,31 Differential GEP studies have provided GEP signatures as a useful prognostic indicator of low-versus high-risk MM; however, none of these gene signatures were treatment specific. 15,18,32 Furthermore, no study thus far has made use of the vast array of HMCLs as model system to generate drug chemosensitivity profile as a representation of the response variation in patient subtypes that may be used to derive a PIspecific GEP signature predictive of resistance and treatment outcomes.
In the current study, using a panel of 50 HMCLs, we successfully developed an in vitro chemosensitivity profile to four PI treatments. We focused on common features that distinguish the most and least responsive HMCLs, noting that these responses are very similar across all four PIs. However, several outliers (HMCLs that were highly resistant to one PI and highly sensitive to another) were noted. Interestingly, several cell lines with intermediate response to Bz or Ix were found highly responsive to Cz treatment. However, in general the ranked correlation with response was very similar among all four PIs, particularly those lines that were collectively the most responsive and least responsive to all four PIs (extraordinary responders). Further characterization of the outliers may reveal important features that better direct which PI is most effective.
Using transcriptomic profiling of PI-sensitive vs PI-resistant HMCLs followed by extensive machine learning-based response prediction, we demonstrate here a gene expression signature that distinguishes strong sensitivity versus strong resistance to the four PIs when used in MM. Importantly, this signature was not applied across the full range of responses, but was very effective in identifying extraordinary (good versus poor) response. We also observed that this 42-gene signature of innate PI response is different from a signature of acquired PI response that we generated by comparing GEP between PI-sensitive and subsequent resistant outgrowths of Bz-resistant MM cell lines (unpublished data).
Despite the small numbers of cell lines used to generate the 42gene signature, when applied to clinical trials on MM patients undergoing PI-based therapies, our GEP signature was remarkably correlated with extraordinary PI response/resistance, could predict early or late progression events and efficiently clustered patient responses into favorable (response or complete response) vs unfavorable (nonresponse) outcome subgroups. Notably, the GEP signature did not stratify response to non-PIs in the APEX, HOVON and CoMMpass trials, suggesting PI treatment specificity, not simply high-versus low-risk tumor biology irrespective of drug treatment regimen as has been previously reported. 32 Several genes included in our GEP signature have been previously implicated with prognosis and disease progression in myeloma. For example, the cell cycle regulator CCND1 (cyclin-D1) has been reported as a strong prognostic predictor in multiple studies. 15,16 HSPA1B is an NRF2-mediated oxidative stress response gene that was also found significant (P o 0.001; |fold change| 42) in our earlier study that compared differential responsiveness with Bz in vitro. 33 Abnormal levels of serum bone Gla protein have been shown to play important role in the inhibition of bone formation in the occurrence of bone lesions in MM that may have prognostic implications. 34 Wolfram syndrome 1 encodes an endoplasmic reticulum (ER) transmembrane protein that plays a role in regulation of ER stress signaling and in the pathogenesis of diseases involving chronic, unresolvable ER stress. 35 Furthermore, Wolfram syndrome 1 is a downstream transcriptional target of the gene X-box binding protein (XBP1), a bZIP transcription factor that regulates unfolded protein response and ER stress-induced apoptosis. 35,36 Higher expression of XBP1 pathway genes is associated with better outcome/response in PItreated MM patients, whereas loss of XBP1 has been shown to induce PI resistance. 36,37 Interestingly, we found that the unspliced transcript of XBP1 (Xbp1u) had~1.5 times lower expression in Ixresistant cell lines when compared with Ix-sensitive cell lines. The transcription factor EB (TFEB) is a master gene for lysosomal biogenesis that drives expression of a large subset of autophagy genes leading to activation of autophagic cell death of myeloma cells owing to metabolic cell stress. 38 RNF170 encodes an ER membrane ubiquitin ligase that mediates ubiquitination and degradation and plays a key role in cell signaling. 39 Given previous reports on the association of these genes to myeloma biology and outcomes, we considered their collective inclusion in our unsupervised clustering further validation of their role in PI response. Although individually associated with MM biology, our GEP signature suggests a combination of effectors influence response to PIs. Additional studies will be required to functionally validate the combined effects of these genes towards response to PI-based therapy.
Earlier, we demonstrated that gene expression signatures may be used to identify secondary therapies in PI-resistant MM using in silico predictions that were confirmed in vitro. 33,40 In 2014, the NCI initiated the Exceptional Responders Initiative to understand the molecular basis of exceptional response to chemotherapy in cancer patients enrolled in clinical trials. 41 The primary goal of the study is to identify molecular features in malignant tissue that may aid to predict response to same or similar drugs. On similar lines, our work should serve as resource to use machine learning-based approaches for the personalized prediction of exceptional chemoresistance and to eventually identify signatures of drug combination regimens that may effectively reverse drug resistance by predicting drugs for various subpopulations/subclones of tumors based on pharmacogenomic signature profiles.