Accurate Outcome Prediction in Neuroblastoma across Independent Data Sets Using a Multigene Signature Free

One of the main challenges in clinical cancer research remains accurate prediction of outcome, enabling better choice of risk-related therapy. This is particularly true for neuroblastoma, a pediatric tumor of the sympathetic nervous system, which is characterized by a remarkably heterogeneous clinical course. Tumors that are found in infants frequently regress spontaneously or show differentiation features on treatment, whereas tumors diagnosed in children >1 year of age often metastasize, causing accelerated cancer-related death despite intensive therapies. Accordingly, different therapeutic schemes exist ranging from watch-and-see approaches to multimodal therapies. Four major risk stratification systems are currently being used in various parts of the world (Europe, United States, Japan, and Germany) based on a combination of clinicopathologic and genetic parameters, such as age at diagnosis, tumor stage, MYCN gene status, histopathologic classification, ploidy, and chromosome 1p and 11q status (1–8). Clinical experience within these systems indicates that the stratification is useful, but misclassifications occur, resulting in overtreatment or undertreatment. Identification of more specific and sensitive markers for response to therapy and outcome prediction is clearly required and is expected to result in better choice of risk-related therapy.

As differences in outcome are considered to reflect underlying genetic and biological characteristics that have their effect on mRNA gene expression profiles, several microarray expression profiling studies have been undertaken to predict patient outcome in different cancer entities.

An important limitation of many published gene expression profiling studies is the lack of statistical power to identify markers and lack of independent validation. Typically, around 30,000 to 40,000 transcripts are tested, generating hundreds of thousands of data points for a relatively small subset of tumors (between 20 and 100). When such a high number of genes are evaluated as prognostic markers, there is a substantial chance that a random association between a gene and the prognostic classes is observed (9, 10). Consequently, many published studies do not classify patients better than chance due to lack of internal validation by repeated random sampling of training sets or external validation on independent samples. As such, there are a few inherent but often overlooked statistical issues, such as data overfitting, unstable gene lists, and lack of study power (11).

In this study, we established a prognostic 42-gene classifier for children with neuroblastoma by reanalysis of four published gene expression studies from four different microarray platforms comprising 582 patients in total (12–15). To facilitate data comparison across different platforms, probe annotations were updated with respect to the original publications. When available, clinical follow-up information was updated. All these aspects critically contribute to the success of our multigene signature. Successful validation of the multigene signature in four independent unpublished data sets shows its robust performance and platform independence.

Four published studies were used for selecting the genes and building the prognostic classifier (phase 1 data sets), and four unpublished data sets were used as independent validation sets (phase 2 data sets).

The phase 1 data sets were downloaded either from the National Center for Biotechnology Information Gene Expression Omnibus (GSE2283 and GSE3960; refs. 14, 15) or from the European Bioinformatics Institute ArrayExpress database (E-TABM-38; ref. 13), or from the authors' Web site¹⁰

A trained multigene correlation signature was validated on the four independent phase 2 data sets from which the 42 genes (when present) were extracted and standardized (per gene, the median value across the samples was subtracted followed by division by the SD of the gene): (a) hgu95av2 Affymetrix gene expression data from 106 neuroblastoma patients (validation set 1; 40 genes present), (b) hgu133plus2 Affymetrix gene expression data from 53 neuroblastoma patients (validation set 2; 40 genes present), (c) data set for 91 neuroblastoma patients obtained using an 11K custom Agilent oligonucleotide microarray (validation set 3; 41 genes present), and (d) Human Exon 1.0 ST Affymetrix expression data from 101 neuroblastoma patients (validation set 4; 42 genes present; standardized data of the 42-gene selection as well as clinical data are available in Supplementary Tables S1 and S2; Fig. 1).

For the remainder of the article, we will label the data sets according to the first author for the published phase 1 studies [Oberthuer (13), Wang (15), Berwanger (12), and Ohira (14)] and as validation sets 1, 2, 3, and 4 for the unpublished phase 2 studies.

To make the data from the different microarray platforms maximally comparable, annotation information of the probes was updated using the MatchMiner tool (16) for the custom-made cDNA or oligonucleotide arrays (12–14) and using the latest version of the R packages hgu95av2 and hgu133plus2 for the Affymetrix array data (15). Probe identification numbers were converted into gene symbols to enable straightforward comparison of the gene lists between the different studies. Throughout the text, the number of unique gene symbols (represented by one or more array probes) in each study is indicated.

Updated clinical information with regard to progression-free survival (PFS) and overall survival (OS) times was obtained from the authors (14, 15) or was publicly available (13). For the Berwanger and Ohira studies and validation set 1, only OS data were available.

Patients were divided in two clearly defined risk groups. The low-risk subgroup was defined by stage I, II, or IVS without MYCN amplification, and the high-risk subgroup comprised patients with age of diagnosis >1 y with stage IV tumors (irrespective of MYCN status) or with stage II and III tumors with MYCN amplification. To develop our classifier, as many patients as possible from the four phase 1 data sets were divided in the two risk groups with maximally divergent clinical course (Table 1), that is, low-risk patients with PFS time (or OS time for Berwanger and Ohira data sets) of at least 1,000 d and high-risk patients that died from the disease. The patients that did not belong to the above-mentioned low- or high-risk subgroups were used as independent test set.

Identification and validation of prognostic classifiers (for each single phase 1 data set) were done by prediction analysis of microarray (PAM) classification with 10-times repeated 10-fold cross-validation in the R statistical language using the Bioconductor package MCRestimate (Fig. 1A; refs. 13, 17). Forty-two genes were present in at least two of the four resulting gene lists.

A cross-platform gene signature was built using standardized expression data of the 42 genes (if present on the respective arrays, see Supplementary Data 2) from four published phase 1 studies. The correlation method was used to build and test a cross-platform prognostic signature (Fig. 1B). Log-transformed data were merged in one file (if more than one probe was present for a certain gene, the probe with the highest expression value was selected), and for each of the 42 genes, the mean expression value in low-risk neuroblastoma patients with PFS of at least 1,000 d was subtracted from the mean expression value in high-risk neuroblastoma patients that died of disease. For classification, the Pearson's correlation coefficient of the signature with the standardized expression values of independent test patients was calculated. Patients with a correlation coefficient below 0 were predicted to have good prognosis, whereas the other patients were predicted to have bad prognosis [according to Liu et al. (18)].

Kaplan-Meier survival analysis was done with the R survival package (R version 2.6.1). The area under the receiver operating characteristic curve (AUC) was used as a measure for the accuracy of the classifiers (ROCR R-package). Multivariate forward conditional logistic regression analysis was done using SPSS version 16.

A complete 10-times repeated 10-fold cross-validation using the PAM algorithm (13, 19) was done on the training patients belonging to one of the two clearly defined risk groups from the four published phase 1 studies separately to identify robust prognostic markers (Fig. 1). This process was accompanied by determination of the classification accuracy, providing a first estimation of the utility of the expression data to predict outcome (Table 1).

For each data set, we selected the probes that were included in at least 65 of the 100 cross-validation gene lists, as these genes are likely to be the ones with the highest prognostic value as determined by Oberthuer et al. (13). The resulting prognostic gene lists from the four studies showed significant overlap (Table 2; Supplementary Data 1). Two genes were in common between three lists (i.e., MYCN and NTRK1), whereas 40 genes were in common between two lists. Thirty-two were previously reported in at least 1 of 10 published prognostic gene lists, of which only 10 were found in 2 or more published prognostic lists (12–14, 20–26). The occurrence of the 42 genes in at least two of the four lists makes them robust, platform-independent, prognostic markers.

Next, we investigated whether the 42-gene list is able to predict prognosis across different data sets. The classification performance was estimated in the different phase 1 data sets using a complete 10-times repeated 10-fold cross-validation method using all patients from the two clearly defined risk groups. For this analysis, it is important to note that not all 42 genes are present on all platforms; hence, the performance test was inherently done with a different number of genes for the different data sets (Supplementary Data 2). As already indicated, the 10-times repeated 10-fold cross-validation provides a good estimate for the classification performance using the expression data of the selected gene list.

As a reference, the 35-, 330-, 81-, and 82-gene lists obtained through single PAM analysis of each of the four phase 1 data sets were evaluated in the same way as the 42-gene list. The classification performance was also tested for a subset of 11 genes (from the 42-gene list) that were present on all four platforms. This analysis showed that all performance parameters for the 42-gene list are best or second best for all studies compared with the other gene lists, whereby the overall accuracy is highest for the 42-gene list subset (AUC = 0.935; Supplementary Data 2). This analysis also shows that the performance of a classifier built for a given data set is not always best, which indicates the power and utility of our meta-analysis for the identification of a prognostic gene list by using expression data of 250 training samples (170 low risk and 80 high risk). When only 11 genes of the 42-gene list were selected that are present on all four platforms, the overall accuracy was lower due to loss in sensitivity and positive predictive value. The 42-gene classifier was also compared with two published classifiers (13, 27) and showed that the 42-gene classifier performs best.

A major disadvantage of the PAM classification method is the need for a training set of samples that are analyzed on the same gene expression measurement platform as the one used to evaluate the test samples. We therefore applied an alternative method to build a classifier based on the 42-gene list that can be used for completely independent data sets even on other platforms.

The prognostic signature is determined using 250 training samples from the four phase 1 studies. A 42-gene classification vector was created and tested using the correlation method (see Materials and Methods; Fig. 1).

First, the correlation signature was tested on the 129 test samples (patients not belonging to the low- and high-risk subgroup) from the four phase 1 studies and revealed a very high predictive power for OS (log-rank P = 2.41E−4) and PFS (log-rank P = 3.40E−7; Fig. 2).

Next, this correlation signature was evaluated on the four independent phase 2 data sets (351 patients), whereby the patients could be clearly separated into groups with significant differences in OS (log-rank P = 2.17E−23) and PFS (log-rank P = 2.03E−21; Fig. 3A). Kaplan-Meier analysis of patients stratified using known risk factors (i.e., age, stage, and MYCN gene status) showed that the correlation signature outperforms these risk factors (P < 0.001, except for MYCN-amplified samples; Supplementary Fig. S2). This was confirmed using multivariate logistic regression analysis evaluating age, stage, MYCN status, and the gene classifier, indicating that the 42-gene signature is an independent predictor for PFS and OS in the four phase 2 data sets as well as in the test samples of the phase 1 data sets (Table 3). Of note, whereas phase 2 data sets are representative of the general neuroblastoma population, test samples from the phase 1 data sets only represent intermediate risk patients.

As the different validation data sets include patients stratified using different risk stratification systems (Europe, United States, and Germany), we defined a common low- and high-risk group (Supplementary Data 3). As there was only 1 patient of 50 that died of disease within the common low-risk group of patients, we did not do Kaplan-Meier analysis. However, we could show that this single patient was classified in the high-molecular risk group using our classifier. Most interestingly, the correlation signature could partition patients within the common high-risk subgroup into groups with significant differences in OS and PFS (Fig. 3B) and was an independent prognostic marker (odds ratios, >4; Supplementary Table S4). To exclude that the significant survival differences in high-risk tumors is solely due to the effect of the MYCN amplification and related downstream MYCN signaling, we also tested the survival in high-risk tumors without MYCN amplification and could show that the classifier also significantly discriminates these patients with respect to outcome (Fig. 3C; Supplementary Table S4). In line with this, inspection of the 42-gene list indicated that not all 42-genes are related to MYCN amplification (Supplementary Data 4).

In this study, we developed and validated a 42-gene prognostic classifier for children with neuroblastoma through a reanalysis strategy of published data complemented with gene expression data from 351 patients from four unpublished data sets (Fig. 1). To accomplish this, four published microarray studies comprising >500 neuroblastoma patients were reanalyzed generating four new prognostic gene lists with a high overlap of genes between them. Comparison of the genes in the classifiers showed that 42 unique genes were present in at least two of the four lists. Not surprisingly, this set of 42 predictor genes contains numerous genes that have been reported in the context of neuroblastoma (e.g., MYCN, NTRK1, NME1, CADM1, FYN, ODC1, and WSB1). The finding of these genes in at least two independent studies indicates their robustness as prognostic markers. Comparison of the performance of the 42-gene list with the lists that were generated on the individual phase 1 studies and with two published prognostic gene lists (13, 27) showed that the classifier based on the 42-gene list has the highest overall accuracy while using the lowest number of genes. However, we have to keep in mind that this observed superiority of the 42-gene set might in part be due to the fact that, for some of the other gene lists, a large proportion of genes were not present on the platform (Supplementary Table S3).

The high prognostic classification performance of the 42-gene list is undoubtedly due to our unique reanalysis approach. First, annotations of the probes on the different platforms were updated according to the latest genome build. Second, a uniform risk definition was applied to select training patients across the different studies. Only patients with maximally divergent courses were used for training. Third, the same powerful algorithm with built-in cross-validation was used for identification of prognostic genes in four major published data sets, enabling the generation of relatively stable prognostic gene lists with high overlap.

This list of 42 prognostic genes was used to build a cross-platform classification signature. As the PAM algorithm is not suitable for cross-platform classification, we used a more intuitive, alternative method for building a 42-gene classifier. In this study, we generated a prognostic correlation signature based on expression data of the 42 genes in 250 training samples of the four phase 1 data sets. The signature was subsequently applied on independent test samples from the phase 1 data sets and on four independent and unpublished phase 2 data sets, generated on different expression profiling platforms, totaling 480 patients. The excellent prognostic performance of the 42-gene list (Table 3) further shows the validity of our meta-analysis approach and the utility of the recognized prognostic markers for neuroblastoma. The classifier predicts overall (OS) and PFS for the patients from the four phase 2 studies as well as for the test patients from the phase 1 studies (which could not be unequivocally classified in the low- or high-risk subgroups using known risk factors) with high sensitivity and specificity (Table 3). Importantly, the classifier was shown to be an independent predictor for both PFS and OS when stratifying for known risk factors such as age, stage, and MYCN status. Indeed, the 42-gene list does not only contain MYCN-regulated genes and, thus, not only reflects the MYCN copy number status of the samples. This is further substantiated by the excellent performance of the classifier in the high-risk neuroblastoma patients without MYCN amplification.

Thus far, this is the largest prognostic meta-analysis study in neuroblastoma, totaling >900 patients, including 351 patients from four independent and unpublished validation data sets. In contrast to other studies on neuroblastoma gene expression classifiers (13, 14, 21, 25, 27, 28), we could show an excellent performance of our classifier on these four independent data sets involving patients from different risk protocols from Germany, Europe, and United States by using a smaller gene set and a more intuitive classification method.

This survival classifier will definitely help to identify patients with increased risk in the current risk groups and to make a better choice of risk-related therapy. For example, low-risk patients with high molecular risk might benefit from more aggressive treatment protocols, whereas more intensive follow-up and new experimental therapies might be considered for high-risk patients with high molecular risk.

In conclusion, we applied a unique meta-analysis strategy for the identification of a robust set of 42 prognostic genes for outcome prediction in neuroblastoma. Furthermore, we propose a prognostic gene signature that is significantly associated with outcome prediction in neuroblastoma samples from independent studies using different technological platforms, making it a useful and practical classifier for risk stratification in neuroblastoma patients. The signature remains to be tested in a prospective clinical validation. The low number of genes makes this signature very well suited for cost-effective and fast PCR-based analysis, requiring only minimal amounts of tumor material, as exemplified by a recently published quantitative PCR study in which a 59-gene classifier containing the 42 genes from this study was trained, tested, and independently validated on a large cohort of patients (29). The outlined strategy for robust selection of prognostic markers and the use of a cross-platform correlation signature have wide application potential in other cancer entities.

In the search of an optimal prognostic classifier, it could prove useful to do an integrated analysis to determine the combined prognostic power of a mRNA gene expression signature along with gene copy number levels, microRNA gene expression patterns, and epigenetic modifications.

No potential conflicts of interest were disclosed.

We thank Vincent Detours for critical reading of the manuscript and helpful suggestions and the Children's Oncology group.

Grant Support: European Community under the FP6 (project: STREP: EET-pipeline, number: 037260). This article presents research results of the Belgian program of Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister's Office, Science Policy Programming. K. De Preter is a postdoctoral researcher with the Fund for Scientific Research-Flanders. J. Vermeulen is a Ph.D. student with the Belgian Kid's Fund and the Fondation pour la recherche Nuovo-Soldati. This study was sponsored by the ‘Kinderkankerfonds,’ the ‘Stichting tegen Kanker,’ the Fund for Scientific Research and BOF-UGent, and the fondation Fournier Majoie pour l'innovation.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.