Gene expression profile alone is inadequate in predicting complete response in multiple myeloma.

With advent of several treatment options in multiple myeloma (MM), a selection of effective regimen has become an important issue. Use of gene expression profile (GEP) is considered an important tool in predicting outcome; however, it is unclear whether such genomic analysis alone can adequately predict therapeutic response. We evaluated the ability of GEP to predict complete response (CR) in MM. GEP from pretreatment MM cells from 136 uniformly treated MM patients with response data on an IFM, France led study were analyzed. To evaluate variability in predictive power due to microarray platform or treatment types, additional data sets from three different studies (n=511) were analyzed using same methods. We used several machine learning methods to derive a prediction model using training and test subsets of the original four data sets. Among all methods employed for GEP-based CR predictive capability, we got accuracy range of 56-78% in test data sets and no significant difference with regard to GEP platforms, treatment regimens or in newly diagnosed or relapsed patients. Importantly, permuted P-value showed no statistically significant CR predictive information in GEP data. This analysis suggests that GEP-based signature has limited power to predict CR in MM, highlighting the need to develop comprehensive predictive model using integrated genomics approach.

Introduction

With the introduction of novel agents, there has been significant improvement in outcome in multiple myeloma (MM).[1-5] Although median survival has improved from 3 years to over 7 years, the clinical course is highly variable and unpredictable. Initially high-dose therapy[6, 7] and subsequently novel agents[8-1112-14] have contributed to improvement in both response and survival[15]; however, MM remains uniformly fatal.[3, 4, 14, 16, 17] Each of these treatments individually and in combination is effective in only a portion of patients, and currently predictability of response to therapy remains unreliable. A number of disease, host and therapy-specific features including age, performance status, tumor burden, tumor proliferative index, serum β-2 microglobulin, albumin, creatinine, LDH, International Staging System (ISS), and cytogenetic abnormalities have been reported to provide prognostic information.[16] Although this and other systems help with prognosis for survival, they do not reliably predict probability of response and as of now, there is no pre-treatment response predictive biomarker to tailor patient-specific therapy.[18]

Moreover, with each new therapy attempted, patients often experience complications and toxicities requiring specialized intervention. With the improved understanding of oncogenomics in MM, a number of targeted agents have been investigated[19, 20]; three such novel agents, thalidomide, bortezomib and lenalidomide have already received FDA approval for use in MM, making selection of effective agents more difficult as well as economically taxing. As some of these novel agents were identified based on molecular identification of targets and the combinations were formulated based upon perturbation of expression profile by these agents, it is anticipated that integrative genomics will help provide more precise prognostic and predictive tools.

Gene expression profiling (GEP) is being widely used for tumor classification[21-25] and survival risk prediction.[26-31] However, it remains unclear whether the strategies used to define such prognostic genomic classifiers can be used to develop classifiers that predict response to specific therapy.[32, 33] Efforts assessing potential use of GEP-based signatures in therapeutic decision making in MM have been limited. Mulligan et al. have studied pharmacogenomics in 169 relapsed multiple myeloma patients samples receiving bortezomib[32] and another study by Kumar et al[34] have attempted response prediction to thalidomide - dexamethasone combination using GEP based signature in newly-diagnosed patients. However, results from these studies are limited either by small sample size[34] or response signature derived from GEP of relapsed and not from newly diagnosed patients with MM.[32] Here, we have analyzed GEP data from one large study of 136 newly diagnosed patients with MM through a rigorous statistical approach, and have independently investigated 3 additional studies of total 349 newly diagnosed, and 162 relapsed cases to determine ability and significance of gene expression signature in predicting complete response (CR) in patients with MM.

Patients and Methods

Datasets for Response Prediction

Gene expression profiles (n=647) from four independent clinical trials were analyzed to predict treatment response. Details about these studies are listed in Table 1. HOVON[35] and Mulligan et al.[32] datasets are publically available, while IFM datasets are now available at NCBI GEO website with accession ID, GSE39754 (IFM I) and GSE55145(IFM II). Patients were newly diagnosed and had not received any prior therapy except those in APEX/SUMMIT trials from Mulligan et al datasets[32] who were relapsed/refractory and had received 1–3 prior therapies. In two of four datasets [IFM I and HOVON], patients received three-drug regimen (VAD/PAD) followed by autologous stem cell transplantation (ASCT), and patients in the other two datasets received bortezomib with dexamethsone [IFM II and Mulligan et al.] prior to response evaluation using uniform European Group for Bone Marrow Transplantation criteria (EBMT).[36]

Table 1

Characteristics and details of studies used for gene expression profile-based response prediction.

	IFM I	IFM II	HOVON	Mulligan et al.
Study	IFM 2005#	IFM 2005#	HOVON 65MM / GMMGHD4$	APEX /SUMMIT
Number of Samples	136	67	282	162
Platform	Affymetrix Exon 1.0 ST array	Affymetrix Exon 1.0 ST array	Affymetrix U133 Plus 2.0 array	Affymetrix U133 Plus 2.0 array
Patient Population	Newly-diagnosed	Newly-diagnosed	Newly-diagnosed	Relapsed
Treatment Protocol	VAD, ASCT	Bortezomib, ASCT	VAD/PAD, ASCT	Bortezomib
Time of Response Measurement	Post-transplant	Post-Induction	Post-transplant	Post Salvage therapy
Complete response	44 (32%)	24 (36 % %)	76 (27 %)	73 (43%)∞

: NCBI GEO accession IDs: IFM I: GSE39754, IFM II: GSE55145.

: Broyl A, et al. Blood 2010

: Includes patients with partial response; 13 patients with CR and 60 with PR.

Data Preprocessing

GEP from samples within all four datasets was obtained on Affymetrix microarray platform. The IFM I and II datasets were obtained using Affymetrix Exon 1.0 ST gene array, while HOVON and Mulligan et al datasets were obtained using HG_U133_Plus_2 gene array and gene-level signal intensities were used for subsequent analysis. Pre-processing and normalization of datasets was carried out using dChip and R package - aroma.affymetrix.

Class Prediction Methods

Class prediction analysis was carried out in each of four datasets using a series of steps, as summarized in Figure 1. Specifically, each dataset was split into training and test datasets, with 2/3 samples in training dataset. Various feature selection approaches were used to eliminate redundant and non-varying gene probe sets, including gene filtering by variance, variable selection using LASSO and Ridge regression, among others. [Suppl. Table 1s] Patients having CR were labeled as positive class for further analysis. Then, classifier models were built and trained within a training set based on gene expression pattern of selected genes (features). Class prediction methods from six major groups of machine learning techniques were used in building a classifier model, including decision tree, support vector machines (SVM), prediction analysis of microarray (PAM), K-nearest neighbors (KNN), Bayesian additive regression trees (BART), and artificial neural networks (ANN). [Suppl. Table 1s] Several classifier models were built using different feature sets and tuning of method-specific criteria to minimize misclassification error. Also, response labels were randomly permuted for 1000 times to build a robust model. A final classifier model was obtained following leave-one-out cross-validation (LOOCV) and/or K-fold cross-validation. The entire model-building process, including gene selection step, was repeated during each cross-validation cycle. Finally, the best-performing cross-validated model was applied to the test dataset for an unbiased assessment of response prediction, and model performance was evaluated by measuring model-specific sensitivity, specificity, positive and negative predictive value (PPV and NPV), and overall model accuracy. Additionally, receiving-operating characteristic (ROC) curves were drawn and compared among different class prediction methods. Machine learning packages in R language and BRB-ArrayTools developed by Dr. Richard Simon and BRB-ArrayTools Development Team were used for the entire analysis.[37] Specific details for each of classification methods were previously described.[38]

Figure 1

Flow-chart showing major steps to develop response (CR) prediction model.

Results

We analyzed the ability of gene expression profile to predict response to therapy in multiple myeloma. We focused on the ability to predict complete response (CR), as the majority of patients (>90%) receiving newer therapeutic modalities and/or high-dose therapy and transplant (HDT) achieve at least a partial response (PR). Moreover, attaining CR is critical to achieving a curative outcome. With this consideration in mind, we first utilized the gene expression signature as a tool to predict CR. Our primary dataset (IFM I) included 136 patients with MM receiving standard dose induction followed by HDT. Forty-four (32%) patients achieved CR in this group. Using random 2:1 sampling division, the training set of 91 patients included 30 achieving CR and the test set of the remaining 45 patients included 14 attaining CR. We applied SVM & KNN class prediction approach to the training set, as detailed in Methods and Figure 1, and observed an overall CR prediction accuracy of 60% with 24 % sensitivity, 77% specificity, 33% PPV, and 69% NPV. Applying the trained model to the test dataset gave accuracy of 62% with 27% sensitivity, 80% specificity, 40% PPV, and 69% NPV. To improve upon these results and to overcome any technical and/or statistical issues that may affect predictability, we reanalyzed the data using a number of different supervised machine learning methods. As seen in Figure 2, none of the analytical methods used significantly improved overall predictive accuracy, which remained between 56 – 78 % in both training and test sets.

Figure 2

Summary of CR prediction performance

Performance of various class prediction methods in predicting CR in the validation dataset using IFM I dataset (a) and HOVON dataset (b). CCP: Compound covariate predictor; LDA: Linear discriminant analysis; KNN: K-nearest neighbor; NC: Nearest centroid; SVM: Support vector machine.

Different microarray platforms do not improve CR prediction accuracy

To evaluate any confounding effect of microarray platforms used on response predictability, we analyzed an independent dataset generated using the Affymetrix HG_U133_Plus_2 platform from 282 patients (HOVON trial) treated with a standard-dose induction regimen followed by HDT. Actual CR rate was 27% (76/282) in this trial. All patients within IFM I and HOVON datasets received uniform treatment, but their MM cell expression profile was performed using different microarray platforms, namely, Affymetrix Exon 1.0 ST and Affymetrix HG_U133_Plus_2 arrays, respectively. We used an identical set of prediction methods as was used for IFM I dataset described above. As noted in table 2, the maximal achievable accuracy was 78% in the HOVON test set.

Table 2

Summary of CR Prediction Analysis

Maximum achievable performance (MAP) in both training and validation sets for all four datasets using seven different class prediction methods. MAP is defined as highest accuracy achieved by one of seven classifiers methods used for a given dataset.

a: Class Prediction Results – Maximum Achievable Accuracy – Training Set
	PredictionMethod	Sensitivity	Specificity	PPV	NPV	Accuracy
IFM I	KNN	24	77	33	69	60
IFM II	KNN	72	56	65	64	65
HOVON	SVM	17	92	48	71	70
Mulligan et al.	KNN	71	67	68	71	69

SVM: Support Vector Machines; KNN: K-nearest neighbors (n=1)

CR prediction accuracy does not improve even if patient selection includes extremes of responses

Next, we evaluated whether gene expression signature using patients achieving CR/nCR versus those with no response (NR) or progressive disease(PD) could improve predictive power. In this case, no patients with PR or MR were included in the analysis. The IFM II dataset had 67 patients, 34 with CR/nCR and 33 with NR/PD. Using K-nearest neighbor (K=1) method, we reached maximum sensitivity of 72% with overall accuracy of 65% in the training data set. However, the same 38-gene signature failed to enrich CR prediction in the test dataset (see below), with reduced sensitivity and accuracy of 67% and 58%, respectively.

To evaluate if early-onset CR has better predictability, we performed CR modeling based on time-points when CR was achieved. Specifically, HOVON study determined CR status at the end of 3rd and 8th cycle of induction therapy (early-onset CR) and at the end of protocol (CR20). Also, we evaluated whether CR can be better predicted in those who have sustained CR, This group was classified under sustained CR group. Using similar analysis as above we did not observe significant improvement in CR prediction in these newly regrouped subsets. [Suppl. Table 2s]. Furthermore, we assessed performance of CR prediction separately in high and low GEP risk groups as defined by proliferation index (PI) and cytogenetic abnormalities.[39-42] In these groups also, our results failed to show significant improvement in CR prediction. [Suppl. File 4 - Appendix]

Finally, we evaluated whether predictive accuracy changes if patients received therapy in the relapsed setting. We analyzed the Mulligan et al dataset using similar methods as above, except that we used PR as a response endpoint since not many CRs were achieved in this relapsed patient population. We achieved an accuracy of 44% in test set. Using all the additional methods described above, we did not improve upon these results.

Permutation to assess the prediction power

Finally, we compared the actual CR achieved by the patients (real CR) versus the enriched CR or positive predictive value from the classifier model giving maximum accuracy in our test set. As seen in table 3, we do not observe significant enrichment of CR compared to actual CR rate. Moreover, we performed a response permutation by randomly assigning the response labels of patients and analyzing the ability to predict. We performed 1,000 such permutations to predict CR. The permuted p-value is the proportion of permutations that give predictive ability higher than the one obtained using the actual response labels. As seen in table 3, none of the data sets have permuted p-value of < 0.05, suggesting that the data from gene expression profile is not adequately informative to predict CR outcome.

Table 3

Permuting class labels to assess the power of predicting CR

During the training of classifiers, the treatment response labels were permuted 1000 times to measure the statistical significance of prediction performance. P-values indicate the proportion of total 1000 permutations giving better prediction performance than the original analysis using real treatment response labels. Hence, higher p-value suggests lower confidence in prediction analysis results using gene expression data.

	Real CR %	Predicted CR:Positive Predictive Value % (p-value– 1000 permutations)
IFM I	32 (44/136)	25 (0.36)
IFM II	36 (24/67)	27 (0.11)
HOVON	27 (76/282)	45 (0.13)
Mulligan et al.	43 (73/162)	43 (0.49)

Predicted CR or Positive Predictive Value is derived from the test or validation dataset analysis using criteria giving the best possible overall model accuracy.

Discussion

In this study, we show that the ability of gene expression profiling (GEP) to predict CR in patients with MM is very limited. We have used uniformly treated patient population, and treatment responses were uniformly measured across all four studies using EBMT Blade Criteria.[36] In our primary dataset, newly-diagnosed patients with MM in IFM I, we found the best accuracy of predicting CR at less than 67% in the test dataset. To confirm our initial observation, we have analyzed 3 different datasets using 2 different microarray platforms, as well as different treatment protocols. Among them, the Mulligan et al. study involved patients with relapsed MM who were refractory to 1–3 previous treatments. We used a set of common feature selection and supervised machine learning methods to build a robust response prediction signature in training set for each study, and evaluated the performance in a test dataset from the same study. In this thorough analysis, we have performed class prediction analysis within each of the four studies to define independent classifier gene signatures to ensure the best predictability within each dataset, and to avoid batch effects when merging different datasets. Despite these efforts, as seen in Figure 2 and Table 2, our response predictability remains low in all the analyses.

To uncover potential information that may reside in the expression data that may allow response predication, we performed permuted prediction analysis. In this approach, we randomly shuffled patients’ response labels and analyzed the ability to predict CR. If the data has some predictive power, then the prediction performance achieved with such random assignment should have significantly lower ability to predict CR than the performance achieved with the real response labels. Following 1,000 such permutations, as seen in Table 3, we observed predictive p-value to be greater than 0.05 in all the four data sets. This indicated limited inherent information in the expression profiles to predict CR.

Our study clearly shows that GEP alone has inherent limitations in CR-prediction in MM. Despite using six major classification methods with various feature selection procedures in building more than 240 models [Figure 2 and Suppl. Figure 2s], the overall prediction accuracy remained low. Furthermore, our analysis reveals inadequate response prediction using GEP regardless of the treatment type used, i.e., patients in IFM I and HOVON trials received VAD followed by ASCT therapy, while patients in IFM II and Mulligan et al. studies received bortezomib therapy. It is possible that GEP of pre and post-treatment may help in feature selection by identifying treatment-specific genes and thereby improving prediction power. However, so far such datasets are not available in public domain. Also, we do not observe significant improvement in CR prediction by stratifying data based on time of achieving CR [Suppl. Table 2s] and GEP based risk-status [Suppl. File 4 - Appenidx].

Among the number of classification methods and data sets used, the highest sensitivity was 67% in IFM II dataset, representing accurate prediction of CR amongst those achieving real CR. False negative rate (1-specificity) remained between 36% – 80%, which is the percentage of patients misclassified as non-CR who might have responded to the treatment. Finally, the sensitivity and overall accuracy dropped significantly below 60% when classifier gene signatures from one dataset was used to predict CR in a different dataset. Among various methods from six major classification groups we used, SVM and KNN performed better than the other methods. [Figure 2]

Based on our analysis, we recognize that the expression profile alone has limited ability in predicting treatment response, especially when actual response rate is high, i.e., prediction using expression profile fails to significantly enrich actual response patients using the data available from the published clinical study. [Table 3] Although we do report sensitivity of more than 70% in some prediction models (IFM II and Mulligan et al.), these values are associated with poor accuracy (65–69%); and importantly, the signature trained from the training dataset performs poorly in the validation dataset of the same study. Previous efforts at GEP-based response prediction in MM have met with very limited success. Mulligan et al.[32] showed overall 71% accuracy in predicting bortezomib response using 169 patients with relapsed MM. A report from Mayo clinic researchers showed inadequate prediction power in a small number of patients receiving thalidomide-dexamethasone combination therapy[34]. Such limited ability of expression profiling to predict response has also been described in breast cancer, leukemia, and other cancer studies.[33, 43, 44] For example, two of seven distant recurrence risk-prediction signatures in breast cancer have recently been commercialized, with test accuracy of more than 80%; however, these signatures are prognostic rather than true predictors of response. Also, risk prediction is restricted to particular subtypes of breast cancer (ER +, Node −). Similarly, a number of risk signatures have been described in myeloma, but they are developed for event-free or overall survival, not correlated with predicting response to initial therapy.[27, 45, 46] In the era of personalized medicine, a response predictive signature will be required to select correct combinations of agents that provide best response with reduced toxicity.

What are the possible biological reasons for not being able to predict response using a large number of gene expression datasets? This analysis highlights the limitation of gene expression by itself to influence the eventual cellular behavior. The transcriptome and proteome modifiers play a significant role in cell growth and survival, especially following therapy.[47, 48] It is well described that alternate splicing and miRNA affect the eventual translation of the expressed gene, and post-translational modifications further control cellular effects of the transcribed proteins.[49-51] Variability in single nucleotide polymorphism also affects response to therapy.[52, 53] Finally, intra-tumoral and inter-tumoral genetic heterogeneity[54, 55] can significantly influence both, individualized response to therapy as well as a major analytical issue for development of GEP-based response signatures.

While the strength of our analysis was in having larger patient cohorts from comparable datasets and use of rigorous and comprehensive statistical classification methods, we acknowledge potential limitations in response prediction imposed by two fundamental statistical challenges in high-dimensional data analysis, e.g., curse of dimensionality (large p for genes) and curse of sparsity (small n for sample size).[44, 56, 57] Nonetheless, we addressed these limits to the best of current knowledge and yet fail to significantly improve response prediction. Also, during the time of this analysis, we have not been able to assess some of GEP datasets with survival information but without response information, which is required for our analysis.

In conclusion, we report that GEP based signature have very limited ability to predict CR in MM which is in accordance with recent consensus report on risk stratification.[40] This analysis highlights the need to develop more comprehensive signatures that can incorporate various genomic elements in final predictive models. Failure in predicting response as an end-point (binary) also warrants need to concurrently measure durability of response (continuous variable), which may improve prediction performance similar to available survival and prognostic GEP based signatures. Finally, integration of next generation sequencing data with available gene expression and miRNA profiling information, along with GEP, may allow development of better predictive models for potential application in clinics.