Pharmacogenomics of bortezomib test-dosing identifies hyperexpression of proteasome genes, especially PSMD4, as novel high-risk feature in myeloma treated with Total Therapy 3

Much of the recent progress in the treatment of multiple myeloma has been attributed to the introduction of several novel agents which, when combined with standard cytotoxic agents and each other, have imparted a high frequency of clinical responses.^1-4  We have previously reported on the superior survival outcomes in Total Therapy 3A (TT3A), with added bortezomib (Bz),⁵  compared with Total Therapy 2 (TT2), which randomly assigned patients to a control arm or an experimental thalidomide arm of a multiple agent chemotherapy and tandem autograft-supported high-dose melphalan program.⁶  The major advance with TT3A was observed in the 85% of patients presenting with gene expression profiling (GEP)–defined low-risk myeloma^7,8  and, considering molecular subgroup designations,⁹  in the 13% exhibiting MMSET/FGFR3 translocations t(4;14) (the MS subgroup).¹⁰  As part of TT3A, pharmacogenomic investigations of Bz, using GEP analysis before and after 48-hour test-dosing (1.0 mg/m²), were performed in an attempt to further delineate low- versus high-risk disease. Strong leads obtained in 142 of 303 patients receiving TT3A motivated the conduct of a successor trial, Total Therapy 3B (TT3B), for validation purposes. The clinical results of both trials have been recently updated.¹⁰  Herein, we report on the Bz-associated gene expression alterations observed in the TT3 trials toward unraveling the mechanisms of action of this first-in-class proteasome inhibitor, the addition of which profoundly improved the outcome of patients with low-risk myeloma as defined by a previously published and validated 70-gene GEP model (GEP70).⁷ 

All patients met the criteria of symptomatic or progressive multiple myeloma according to International Myeloma Working Group 2006 criteria.¹¹  Standard prognostic parameters included β2-microglobulin, albumin, C-reactive protein, lactate dehydrogenase (LDH), M-protein quantifications in serum and urine for response assessment, BM examinations with metaphase karyotyping to determine the presence of cytogenetic abnormalities, and GEP.

Characteristics of patients with GEP data both at baseline and after Bz in the 2 TT3 trials were comparable, except for overrepresentation of low albumin and CD-1 and CD-2 subgroup classifications in patients in the TT3B trial (Table 1). Despite these less favorable features, clinical outcomes were comparable to those of patients enrolled in TT3A, whether examined among all 177 and 303 patients, respectively,⁷  or in the subset of 128 and 142 patients with GEP data both at baseline and after Bz.

Progression-free survival (PFS) and overall survival (OS) durations were measured from the time of initiation of protocol therapy; events included relapse or death from any cause in the former and death from any cause in the latter. OS and PFS curves were estimated by the product-limit method¹²  and were compared by the log-rank test.¹³  Univariate and multivariate regression analyses of covariates and time-to-event outcomes were performed by Cox regression.¹⁴ 

In keeping with institutional, federal, and Helsinki Declaration guidelines, all patients signed a written informed consent both for receiving protocol-based therapy and for undergoing repeated BM sampling for GEP analyses. Both protocols and their revisions had been reviewed and approved by the Institutional Review Board at the University of Arkansas for Medical Sciences, and annual progress reports on translational research investigations have been submitted. According to National Cancer Institute requirements, ∼ 80% of patient records have been audited by an external team of experienced independent auditors.

The details of the 2 TT regimens (TT3A, 2003-33; TT3B, 2006-66) have been published.^5,10  Briefly, each protocol called for 2 cycles of induction therapy with VTD-PACE, comprising Bz (Velcade), thalidomide, and dexamethasone and 4-day continuous infusions of cisplatin, doxorubicin, cyclophosphamide, and etoposide, before melphalan-based high-dose therapy with autologous hematopoietic progenitor support, followed by 2 consolidation cycles of dose-attenuated VTD-PACE. Maintenance regimens differed in that TT3A called for VTD in year 1 and TD in years 2 and 3, whereas TT3B used VD and lenalidomide (Revlimid) instead of thalidomide for all 3 years. TT2 data were used for validation purposes.^6,15 

GEP sample procurement and processing, as well as analyses for the derivation of GEP70-based high-risk designation,⁷  molecular subgroup classification,⁹  and delTP53 assignment,¹⁶  have been reported previously. For Bz test-dosing, consenting patients received a bolus dose of Bz 1.0 mg/m² intravenously and, 48 hours later, a further BM sample was obtained under local anesthesia from the posterior iliac crest. Plasma cell (PC) isolation from mononuclear cell fractions was performed by immunomagnetic bead selection with monoclonal mouse anti–human CD138 antibody with the use of the AutoMACS automated separation system (Miltenyi-Biotec). PC purity of > 95% homogeneity was confirmed by 2-color flow cytometry with the use of CD138⁺/CD45⁻ criteria (Becton Dickinson), immunocytochemistry for cytoplasmic light-chain immunoglobulin, and morphology by Wright-Giemsa staining. RNA was extracted with the QIAGEN RNeasy kit. cDNA was prepared and biotinylated with the Affymetrix GeneChip HT 3′ IVT Express Kit. Samples were hybridized to an Affymetrix U133Plus2.0 microarray according to the manufacturer's recommendations and then read on a GeneChip Scanner 3000 (Affymetrix). The data files were deposited in the ArrayExpress archive (http://www.ebi.ac.uk/arrayexpress) under the accession number E-TABM-1138.

Tricolor interphase FISH analysis was performed so that we could correlate GEP data with FISH-derived 1q21 copy number.¹⁷  BM aspirates were first subjected to Ficoll-Hypaque gradient-centrifugation separation to remove erythrocytes. The cells were immediately fixed in 100% ethanol and stored at −20°C. FISH probes were labeled by a nick translation reaction. Briefly, 1 μg of BAC DNA (RP11-180N18 for the AHCYL1 locus [1p13.3], RP11-498A2 for the CKS1B locus [1q21.3]) was mixed with 10μM fluorescent Green/Red-dUTP, 5μM dTTP, and 15μM each of dATP, dCTP, and dGTP in a nick translation buffer system including NT enzymes (Abbot Molecular Inc). After coprecipitation with 10 μg of human Cot-1 DNA (Invitrogen Corp), the labeled DNA probes were suspended in 100 μL of DNA in an in situ hybridization solution (Dako Corp). The 2 probes were mixed 1:1 before the hybridization. At least 100 tumor cells each were counted for the copies of 1p13.3 (green) and 1q21.3 (red). The threshold of significance abnormality (gain or loss) of each probe was set at ≥ 90%, as we previously described.¹⁷ 

Data were analyzed with the Ingenuity Pathways Analysis online tool (Ingenuity Systems, www.ingenuity.com). All 80 genes from the GEP80 model were considered for analysis. Canonical pathways analysis identified the pathways from the Ingenuity Pathways Analysis library of canonical pathways that were most significant to the dataset. The significance of the association between the dataset and the canonical pathway was measured in 2 ways. First, a ratio was calculated by dividing the number of molecules from the dataset that mapped to the canonical pathway by the total number of molecules that mapped to the pathway. Second, the Fisher exact test was used to calculate a P value.

Affymetrix U133Plus2.0 microarrays were preprocessed and normalized with the MAS5.0 method. With the use of TT3A protocol data as a training set (n = 142 patients with GEP data both at baseline and after Bz), we identified 1051 genes whose expression levels significantly changed 48 hours after the Bz test dose relative to baseline by applying a paired t test on each gene and a .005 P-value cutoff with a false discovery rate (q value) of 0.22.¹⁸  We then correlated each of these 1051 genes with PFS in a Cox regression model and ranked the genes by P values. To explore the various potential survival implications, 2 sets of covariates were considered in the Cox models for each gene: model a (percentage of change at 48 hours compared with baseline expression) and model b (48-hour expression).

Fitting model a on each of the 1051 genes produced a minimum q-value of 0.94,¹⁸  which means that, even if we selected only the most significant gene, the chance of it being a false positive was 94%. Fitting model b produced much smaller q-values than model a, with a minimum of 0.005, indicating that the 48-hour expression, compared with the percentage of change, was more strongly associated with PFS. Accordingly, we chose model b over model a and ranked genes by the P values of the 48-hour expression levels.

Once genes were ranked, PFS was predicted by computing a summary score on the basis of the top-ranked genes for each patient as described previously.⁷  The score for a particular patient was defined as the mean differential expression of prognosis-unfavorable genes and prognosis-favorable genes, whereby the latter were those associated with a hazard ratio (HR) < 1 (“good” genes) and the former with a HR ≥ 1 (“bad” genes). To produce a robust model, we used a 10-fold cross-validation approach (supplemental Table 1, available on the Blood Web site; see the Supplemental Materials link at the top of the online article) to simultaneously select the number of genes to calculate the gene score and an optimal cut point to dichotomize that score. In brief, we examined a range of top genes (from the top 5 to the top 150) to calculate a score and a range of percentiles (from the 50th to the 90th) to dichotomize that score. For each top gene set and percentile, a log-rank test was conducted comparing the predicted high-/low-risk groups by cross validation. The optimal score with the lowest log-rank P value used 80 genes to compute the score and dichotomized at the 82nd percentile. After cross validation, the 80 genes and 82nd percentile were applied to the entire training set to find the optimal cut point of 2.48. Patients having scores greater than the cut point were deemed to have high-risk disease and patients with scores below the cut point low-risk disease.

Once model selection was completed, we then examined the independent prognostic power of the post-Bz GEP80 [GEP80(PB)] in multivariate analysis of PFS and OS, along with other prognostic factors. This score was then applied to the test set of the TT3B protocol (n = 128 patients with GEP data both at baseline and after Bz, the latter of which had also been preprocessed and normalized with the Affymetrix MAS5.0 algorithm, followed by log2 transformation). The risk score was again calculated as the mean expression of the prognosis-unfavorable genes minus the mean expression of prognosis-favorable genes (Table 2), as with the training set. Each patient of the test set was then classified as having high-risk disease if the score was ≥ 2.48 and low-risk disease otherwise.

Using the training set of 142 patients in the TT3A trial, we identified 80 differentially expressed genes (Table 2) whose 48-hour expression levels were highly associated with PFS. The heat maps of these 80 genes are portrayed in Figure 1A, in which patients are grouped according to ascending GEP80 scores at 48 hours (upper panel) and at baseline (lower panel). The heat map of the 48-hour expression shows 2 major gene clusters (separated by the yellow horizontal line) with the upper gene cluster comprising many genes coding for subunits of the proteasome. The 26 patients in the high-risk group showed a concerted up-regulation of 42 bad genes (including genes coding for the proteasome) and down-regulation of 38 good genes. The clinical course for these patients was characterized by highly significantly inferior PFS and OS (Figure 1B). Results were validated in the test set of 128 patients receiving TT3B, showing similar gene expression patterns (Figure 1C) and survival differences (Figure 1D).

We next examined whether application of the GEP80 model to baseline data [GEP80(BL)] in both TT3 trials had outcome-discriminatory power. Results showed significant survival differences both in the TT3A training set (Figure 2A) and in the TT3B test set (Figure 2B). In the training set, 88% and 87% of low-risk patients, respectively, were alive and progression free at 2 years compared with 46% and 23% of the patients with high-risk disease. Applied to the test set, 2-year OS rates were 92% and 60% and PFS estimates were 87% and 53%, respectively, in low- and high-risk groups. Applying the GEP80 model to the combined TT3 protocols showed no significant differences in OS whether the model was applied to samples after Bz or at baseline (data not shown).

The profound predictive power of the GEP80 model in TT3 led us to investigate its discriminatory power in TT2, which randomly assigned 668 patients between a control arm and an experimental arm with thalidomide, both without Bz. Results showed significant differences in both OS and PFS when both arms of TT2 were considered together (Figure 2C), as well as separately (Figure 2D).

Given the profound outcome differences according to the GEP80 risk model, we next examined its performance in the context of the well-validated GEP70 risk model. Applying the GEP80 model to GEP70-defined low- and high-risk disease showed significant segregation in survival outcomes (Figure 3). Thus, among the 221 patients with GEP70-defined low-risk disease, the GEP80 model, when applied to samples after Bz, identified 19 patients (9%) with a grave prognosis. The 2-year estimates of PFS and OS were 63% and 68%, respectively, in these patients and 89% and 92%, respectively, in the 202 patients with GEP80-defined low-risk status. Conversely, among the 49 patients with GEP70-defined high-risk myeloma, the GEP80 model identified 20 patients (41%) with a favorable prognosis. These patients had 2-year PFS and OS estimates of 84% and 83%, respectively, compared with 45% and 59%, respectively, in the 29 patients in whom GEP80-defined risk was high (Figure 3A).

Such survival differences were also noted when the GEP80 model was applied to baseline samples (Figure 3B). The GEP80 model identified 6 patients (3%) with grave prognosis among the 221 patients in the GEP70 low-risk group. Thus, 2-year PFS and OS were 88% and 91%, respectively, in the 215 patients with GEP80(BL) low-risk disease, compared with 33% and 50%, respectively, in the 6 patients with GEP80(BL) high-risk status. However, among the 49 patients with GEP70 high-risk disease, 27 (55%) had a favorable prognosis on the basis of the GEP80 model. These patients had 2-year PFS and OS estimates of 77% and 81%, respectively, as opposed to 41% and 55%, respectively, among the 22 patients with GEP80(BL) high-risk disease. Such discriminatory power of the GEP80 model was not observed when applied to GEP70-defined subgroups treated with TT2, whether applied to all patients or examined in the context of the 2 arms of the trial (data not shown).

To identify predominantly affected pathways we first used the Ingenuity Pathway Analysis online tool to further explore the 80 selected genes of the GEP80 model. This analysis showed significant enrichment for genes of the protein ubiquitination pathway (Figure 4). We also computed overlaps of our 80 genes with curated gene sets with the use of the Molecular Signature Database of the Gene Set Enrichment Analysis online tool.¹⁹  Computing overlaps with a total of 5411 gene sets, the top 4 gene sets thus identified all involved the proteasome^20-23  (Table 3). Further analysis showed that the Bz-induced gene expression alterations were unique to this proteasome inhibitor and were not observed in similar pharmacogenomic investigations after test-dosing with dexamethasone, thalidomide, lenalidomide, and melphalan (Figure 5A; supplemental Tables 2.1 and 2.2). Comparing the GEP80 and GEP70 risk models we found only 3 genes, 2 unfavorable and 1 favorable, common to both (Figure 5B). One of these genes, PSMD4, is one of the non-ATPase subunits of the proteosomal 19S regulator and maps to chromosome 1q21, the amplification of which had been linked to inferior clinical outcomes.²⁴  The other 2 genes, BIRC5 and KIAA1754, map to chromosomes 17q25 and 10q25, respectively. The former is a member of the inhibitor of apoptosis gene family whereas KIAA1754, the only favorable gene found in common, encodes an inositol triphosphate receptor interacting protein.

Correlating the GEP80 model with FISH-derived 1q21 copy number showed that PSMD4 expression levels increased significantly with FISH-derived 1q21 copy number (Figure 6A), as did GEP80(BL)–defined risk scores (Figure 6B). The latter finding also persisted in the context of combined GEP70 and GEP80 risk designations (low/low, low/high, high/low, and high/high), whether applied to samples at baseline (Figure 6C) or after Bz (Figure 6D). We then identified optimal cut points in PSMD4 expression levels, such that patients with 4 copies were designated as having high PSMD4 expression, patients with 3 copies medium expression, and those with 2 copies low expression. The cross-validated error rate was estimated to be 16%. In TT3, both PFS and OS progressively shortened with an increase in copy number–linked PSMD4 expression (Figure 6E). In the case of TT2, PFS was equally inferior in patients with 3 and 4 copy numbers compared with the 2-copy number group; a more graded effect, as in TT3, was noted for OS (Figure 6F).

When applied to both TT3 protocols, low albumin and hemoglobin concentrations; high levels of β2-microglobulin, LDH, C-reactive protein, and creatinine; and the presence of cytogenetic abnormalities conferred inferior PFS and OS on univariate analysis (Table 4). High-risk designations by GEP70 or GEP80 models (whether applied to samples at baseline or after Bz), GEP-derived delTP53, and high Proliferation Index represented highly adverse features for PFS and OS. FISH 1q21 copy number–discretized PSMD4 levels imparted inferior PFS and OS, whereas the FISH-derived CKS1B 1q21 copy number had borderline deleterious effects on PFS and significantly affected OS. In the multivariate stepwise Cox regression analysis, GEP80(BL)–defined high-risk status was the only surviving adverse genetic variable, along with high LDH and low albumin, that predicted inferior OS and PFS.

Pharmacogenomic investigations after test-dose applications of Bz have shown that high expression levels of proteasome genes are linked to inferior prognosis in both of our Bz-based TT3 protocols. The GEP80 model also discriminated outcomes whether applied to samples at baseline or after Bz. These findings, observed in TT3A, were validated in TT3B and also applied to both arms of TT2 (control and thalidomide). With the use of sample before Bz and at baseline in TT3, the GEP80 model provided segregation of high-risk subsets of 9% and 3%, respectively, in the GEP70 low-risk group and low-risk subsets of 41% and 55%, respectively, in the GEP70 high-risk setting. In TT2, the GEP80 model failed to discern a low-risk subset among the patients with GEP70 high-risk disease, and the GEP80 model discriminated only one patient with high-risk status among patients identified to have GEP70 low-risk disease.

To predict the risk category of a patient with newly diagnosed disease who will be treated with a Bz-based regimen, one can apply the GEP80 model to the patient's log2-transformed gene expression profiles that have been preprocessed and normalized with the Affymetrix MAS5.0 algorithm. The risk score is then calculated as the mean log2 expression of those unfavorable probe sets minus the mean log2 expression of the favorable probe sets listed in Table 2. The patient should be classified as having high-risk disease if the score is ≥ 2.48; otherwise, the patient should be considered to have low-risk disease.

PSMD4 was only 1 of 3 genes common to both GEP70 and GEP80 models. PSMD4 and other proteasome genes were uniquely up-regulated by Bz but not by dexamethasone, immunomodulatory agents, and melphalan. Correlative FISH analysis of CKS1B copy number showed that PSMD4 expression levels were highly sensitive to copy number. This enabled FISH-guided distinction of 3 PSMD4 expression levels (high, medium, low) that had prognostic implications in both TT3 and TT2, mainly by segregating a superior subset in case of low-tertile expression levels (Figure 6E-F). Prognostically, GEP80(BL)–defined high-risk status was the sole genetic parameter that survived multivariate PFS and OS models, along with low albumin and high LDH levels, although all the other standard and both cytogenetic and molecular genetic variables also contributed when examined alone.

We have previously reported that both the copy number and the percentage of cells with amp1q21 invariably increase comparing samples obtained at diagnosis and at relapse.²⁴  Because genes residing on chromosome 1q21 contribute critically to the high-risk designation in the GEP70 model, an increase in risk score from baseline to relapse portends poor survival after relapse.²⁵  These data suggest that clonal evolution in myeloma can be traced to a copy number–dependent increased expression of genes within the 1q21 amplicon. We found a significant up-regulation in the expression of proteasome genes after a 48-hour Bz test dose (Figure 5A; supplemental Tables 2.1 and 2.2). These changes, which ranged from 1.1- to 1.2-fold, could reflect the active up-regulation of these genes in the CD138-selected cells in response to suppression of the proteasome. However, an alternative explanation, and one that we support, is that these changes reflect the differential sensitivity of PCs to Bz-induced cell death within 48 hours of Bz dosing.

CD138 selection does not discriminate between normal and malignant PCs, nor does it distinguish clonal myeloma cells with 2, 3, or ≥ 4 copies of 1q21.²⁴  We have reported that myeloma cell samples often exhibit an intraclonal heterogeneity with respect to amp1q21 and that a higher copy number at diagnosis portends a higher risk of death.²⁴  More importantly, we have noted that the proportion of cells with amp1q21 always increases at relapse, suggesting that cells with more copies of 1q21 have a resistant phenotype.

Consistent with this finding on serial samples at diagnosis and relapse, we have also noted that the GEP70 score inevitably increases at relapse and that a shift from low to high risk imparts a significantly shorter survival after relapse.²⁵  Recent murine studies showed that Bz treatment reduces normal PCs in the spleen and BM by > 95%.²⁶  In the current study, we demonstrated that PSMD4 expression is copy number dependent, with normal PCs from healthy donors having the lowest levels and > 90% of clonal cells from myeloma patients exhibiting > 2 copies of 1q21. Taken together, we propose that the minor yet significant up-regulation of PSMD4 and other proteasome genes 48 hours after Bz dosing is not because of a drug-induced up-regulation of the gene but rather to the preferential killing of normal PCs and myeloma cells with 2 copies of 1q21 (Figure 7A).

We have been able to show that high-risk myeloma in the Bz era is associated with high expression levels of proteasome genes. If proteasome activity can be completely extinguished by Bz, how might high-risk disease be related to Bz resistance? Given the rapid recycling of the proteasome,²⁷  high levels of mRNA may allow for the rapid rescue of proteasome activity, so that higher doses or shorter intervals between doses of Bz may be required for patients with high-risk disease to benefit. Alternatively, protein translation inhibitors may prevent the creation of said proteasome subunits. In addition, the extra-proteasomal effects of PSMD4 may counter the activity of Bz.²⁸  Thus, as a polyubiquitin receptor, PSMD4 could bind to and interfere with the function of tumor suppressor genes, causing them to be destined for degradation; however, in the absence of 1q21 gains, Bz would be able to counteract these effects.

Another possibility is that this protein can act as a chaperone-like molecule during the Bz-induced unfolded protein response (UPR) and prevent the subsequent UPR-induced apoptosis of myeloma cells.²⁹  Given the secretory nature of PCs, another intriguing possibility is that PSMD4 fragments can suppress protein secretion. PSMD4 can undergo proteolytic cleavage, and an internal fragment, identified as antisecretory factor 1, has been shown to be a potent inhibitor of mucus secretion in the gut.³⁰  It is conceivable that antisecretory factor 1 exerts its antisecretory function at the transcriptional level and, therefore, could possibly suppress protein production by myeloma cells, thus reducing the load on the endoplasmic reticulum and proteasome. There is a growing body of literature to indicate that both differentiation and death of PCs is regulated by the load and capacity of the proteasome and UPR.³¹ 

Our study has clearly shown that myeloma cells with high proteasome mRNA activity, as evidenced by high PSMD4 expression, portend an inferior prognosis. Because protein production is often reduced or even lost in end-stage disease, we may face, in high-risk myeloma, a scenario of very high proteasome activity and, hence, high protein-degradation capacity but low protein load, resulting in reduced UPR-induced apoptosis by Bz. This mechanism of Bz resistance can potentially be overcome by targeting the UPR with other novel agents in combination with proteasome inhibitors, such as inhibitors of heat-shock protein 90³²  and histone deacetylase.³³  We have also reported that there is global up-regulation of miRNAs in high-risk myeloma, which in turn is associated with decreased protein translation,³⁴  so that strategies targeting miRNAs in high-risk multiple myeloma may restore Bz sensitivity by increasing protein translation and, consequently, the UPR (Figure 7B).

In summary, investigations of the short-term effects of Bz in TT3 have shown insights into novel mechanisms of resistance to this agent. As we have embarked on GEP-defined risk-based assignments of therapy in TT4 and TT5,³⁵  we are also performing pharmacogenomic analyses after test-dosing of melphalan (10 mg/m²) toward identifying the basis for synergistic interaction of this genotoxic drug with novel agents.³⁶ 

The authors thank the M. J. Murdock Charitable Trust for providing support for a genomics server at Cancer Research And Biostatistics (CRAB), Seattle, where part of the statistical analysis for this research was conducted.

This work was supported by the National Cancer Institute (grant CA 55 813).

National Institutes of Health

Contribution: J.D.S. and B.B. conceptualized the work; B.B., S.U., C.J.H., J.D.S., and J.E. wrote the paper; B.B., F.v.R., S.U., S.W., B.N., Y.A., and E.A. contributed patients and performed clinical research; O.S. and R.W. performed gene expression profiling; E.T., P.Q., and I.H. supervised and discussed gene expression profiling analyses; and J.C., P.Q., Q.Z., and Y.Z. performed statistical analyses.

Conflict-of-interest disclosure: J.D.S., B.B., P.Q., and Y.Z. have filed patents related to technology that is described in the present study. J.D.S. is a founder and has an ownership stake in Signal Genetics LLC, a biotechnology company that has licensed said technology from the University of Arkansas for purposes of commercial development purposes; holds patents, or has submitted patent applications, on the use of GEP in cancer medicine; receives royalties related to patent licenses from Genzyme Novartis and Signal Genetics; has received research funding from Celgene, Millennium, and Novartis; has advised Celgene, Genzyme, Millennium, and Novartis; and has received speaking honoraria from Celgene, ArrayBioPharma, Centocor Ortho Biotech, Genzyme, Millennium, and Novartis. B.B. has received research funding from Celgene and Novartis, is a consultant to Celgene and Genzyme, has received speaking honoraria from Celgene and Millennium, and is a co-inventor on patents and patent applications related to use of GEP in cancer medicine. The remaining authors declare no competing financial interests.

Correspondence: John D. Shaughnessy Jr, 4301 W Markham, Little Rock, AR 72205; e-mail: shaughnessyjohndjr@uams.edu.