Predicting survival times for neuroblastoma patients using RNA-seq expression profiles

Predicting survival times for neuroblastoma patients using RNA-seq expression profiles Background: Neuroblastoma is the most common tumor of early childhood and is notorious for its high variability in clinical presentation. Accurate prognosis has remained a challenge for many patients. In this study, expression profiles from RNA-sequencing are used to predict survival times directly. Several models are investigated using various annotation levels of expression profiles (genes, transcripts, and introns), and an ensemble predictor is proposed as a heuristic for combining these different profiles. Results: The use of RNA-seq data is shown to improve accuracy in comparison to using clinical data alone for predicting overall survival times. Furthermore, clinically high-risk patients can be subclassified based on their predicted overall survival times. In this effort, the best performing model was the elastic net using both transcripts and introns together. This model separated patients into two groups with 2-year overall survival rates of 0.40 ± 0.11 (n = 22) versus 0.80 ± 0.05 (n = 68). The ensemble approach gave similar results, with groups 0.42 ± 0.10 (n = 25) versus 0.82 ± 0.05 (n = 65). This suggests that the ensemble is able to effectively combine the individual RNA-seq datasets. Conclusions: Using predicted survival times based on RNA-seq data can provide improved prognosis by subclassifying clinically high-risk neuroblastoma patients. Reviewers: This article was reviewed by Subharup Guha and Isabel Nepomuceno. Keywords: Accelerated failure time, Sparse PLS, Lasso, Elastic net, Data imputation Background been one of the most important markers for stratify- Neuroblastoma is the most frequently diagnosed cancer ing patients. Genome-wide association studies have found in the first year of life and the most common extracra- many other SNPs associated with an increased risk of nial solid tumor in children. It accounts for 5% of all neuroblastoma. However, while aberrations of these genes pediatric cancer diagnoses and 10% of all pediatric oncol- indicate an increased susceptibility to the disease, these ogy deaths [1]. These numbers have improved over the markers are less useful for stratifying patients into risk past decade, but accurate prognosis for the disease has groups after diagnosis. remained a challenge [1]. The difficulty is due to the highly The Children’s Oncology Group stratifies patients into heterogeneous nature of neuroblastoma; cases can range three risk groups using the International Neuroblastoma from tumors that spontaneously regress on their own, to Staging System (INSS) and various prognostic markers aggressive tumors that spread unabated by treatment. including age at diagnosis, tumor histology, MYCN ampli- In 1984, the MYCN oncogene was identified as a fication, and DNA ploidy. According to the American biomarker for clinically aggressive tumors [2]. It has since Cancer Society, the 5-year survival rate for these low-risk, intermediate-risk, and high-risk groups are > 95%, 90% - 95%, and < 50%, respectively. The high-risk group typi- cally consists of patients older than 18 months with INSS *Correspondence: somnath.datta@ufl.edu stage 4 or patients of any age with MYCN amplification. Department of Biostatistics, University of Florida, 2004 Mowry Rd, 32611, Gainesville, USA © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Grimes et al. Biology Direct (2018) 13:11 Page 2 of 15 Predicting survival outcomes using gene expression data PLS is implemented using the “spls” R package [6]. The has been explored with promising results [3, 4]. These number of latent factors v is a tuning parameter, which studies use gene expression profiles with classification is determined from 10-fold cross validation. The optimal methods to stratify patients into risk groups. However, value of v is searched over v = 1, ... , 10. patients that are clinically labeled as high-risk pose a SPLS particular challenge, and classifiers tend to struggle sep- Like PLS, the sparse partial least squares (SPLS) also con- arating those patients into subgroups. In this paper, we structs latent factors, but it incorporates L regularization take the approach of modeling survival time directly using in the process [7]. This induces sparsity in each linear RNA-seq data. This leads to two objectives: the first is to combination of the original covariates that make up the evaluate the accuracy of the model in predicting exact sur- latent factors. There are two tuning parameters, the num- vival times. The second is to determine whether the pre- ber of latent factors v < n and the shrinkage parameter dicted times can be used to subclassify high-risk patients η ∈ (0, 1) for the regularization. Both of these are deter- into distinct groups. mined from 10-fold cross validation using the “spls” R package [6]. The optimal values of v and η are searched Methods over the grid of points with v =1, ..,10 and η = 0.1, ...,0.9. Accelerated failure time (AFT) model Note, to implement PLS the shrinkage parameter, η,is Theacceleratedfailuretime(AFT)modelrelates thelog set to zero. survival times to a linear combination of the predictors. Lasso log(y) = Xβ + ,(1) The least absolute shrinkage and selection operator (lasso) fits the model using least squares subject to an L con- where y ∈ R denotes the vector of n observed sur- straint on the parameters |β |≤ λ,where λ> 0is j=1 vival times, X the n × p matrix with columns containing a tuning parameter that affects the amount of shrinkage the predictor variables for each observation, β ∈ R the [8]. This constraint induces sparsity in the estimated coef- vector of regression coefficients, and  ∈ R avector ficients, setting many coefficients to zero and shrinking of independent random errors with an unspecified dis- others. tribution that is assumed to be independent of X.The The model is fit using the “glmnet” R package [9], which predictors X are centered and scaled so that each col- performs 10-fold cross validation to select λ. umn X ,for i = 1, ... , p, has zero mean and unit variance There are two challenges to fitting this model: the high Elastic net dimensionality of X and the right censoring of y.Since The elastic net (elnet) uses a similar approach as the lasso. It combines both L and L penalties; the estimator p > n, ordinary least squares (OLS) should not be used as 1 2 minimizes the convex function it will simply overfit on the data. Instead, four approaches for dimension reduction are considered, which include 1 1 2 2 ||Y − Xβ|| + λ (1 − α)||β|| + α||β|| ,(2) both latent factor and regularization techniques. To han- 2 2 2 2 dle right censoring, a nonparametric, iterative imputation where λ> 0and α ∈ [0, 1] are two tuning parameters procedure is proposed, which allows the model to be fit as [10]. When α = 1, this reduces to the lasso. By includ- though complete data were available. ing some component of the L penalty, groups of strongly Each of the dimension reduction techniques require the correlated variables tend to be included or excluded in the selection of one or more tuning parameters. These param- model together. The “glmnet” R package [9]isusedtofit eters are determined by 10-fold cross validation, which the model and determine both tuning parameters. is implemented in R using two packages discussed in the following sections. Imputation for right censoring PLS Let (y , δ , X )|i = 1, ... , n denote the set of observed i i i With partialleast squares(PLS),acollectionof v < n survival times, indicators for death from disease, and the orthogonal latent factors are computed as linear combi- p-dimensional vector of covariates for the n patients in the nations of the original covariates. The construction of the dataset. Let T denote thetruesurvivaltimesforpatient latent factors takes into account both X and y;thisisin i = 1, ... , n.Ifthe ith patient’s survival time is censored contrast to principal component analysis (PCA), which (i.e. δ = 0) then we only observe y < T .Thatis, T is i i i i only considers X. A review of PLS and its application to unobserved. genomic data can be found in [5]. Once the v latent factors To deal with this right censoring, the dataset imputa- are computed, the AFT model is fit using OLS with these tion procedure from [11] is used. This procedure is briefly (0) new variables. summarized here. To begin, an initial estimate β is Grimes et al. Biology Direct (2018) 13:11 Page 3 of 15 obtained by fitting the AFT model using only the uncen- censored data is imputed once at the start of the ensemble sored data. Then, in each of k = 1, ... , n iterations, do training; the censored survival times are replaced with the the following. predicted times from the single best model (TI-4). (k) 1 Calculate the Kaplan-Meier estimate S (e) of the Classification: LPS vs. non-LPS { } distribution of model error using (e , δ )|i = 1, ... , n i i The second goal is to subclassify clinically high-risk T (k−1) where e = log(y ) − X β . i i patients. A new dichotomous variable is created to classify 2Impute n new datasets by replacing each censored patients: If the predicted survival time is less than t > 0 T (k−1) ∗ ∗ log(y ) with X β + e ,where e is a sampled i i i years, we say the patient has low predicted survival (LPS). model residual from the conditional distribution Otherwise, the patient is non-LPS. For patient i = 1, ..., n (k) S (e|e > e ). This condition ensures that the with predicted survival time y ˆ ,let imputed observation will be larger than the observed 1if y ˆ ≤ t right-censored time. i LPS = .(3) i,t 0otherwise 3 Use the new datasets to compute n new estimates (k) β for j = 1, ... , n . Two cutoffs were considered with t = 2and t = 5years. 4 Average the n estimates to obtain a final estimate For clinically high-risk patients, the t = 2cutoffisuse- n (k) 1 D (k) ˆ ˜ β = β . ful for identifying those with a significantly lower survival j=1 j rate. In the general population of neuroblastoma patients, The process is repeated for n iterations, and the final the t = 5 cutoff is useful as an alternative way to iden- (n ) ˆ K estimate β is returned. tify high-risk patients, but it cannot tease out the more To balance between computation time and simulation extreme cases. variability, we chose to run n = 5 iterations, imputing n = 5datasetsineach. Performance measures Performance is evaluated on the testing dataset by four Ensemble method different measures. The ensemble method incorporates bagging with rank The first involves the prediction error of survival times. aggregation over each performance measure. The 12 This is measured by the root mean squared error, adjusted models using genes, transcripts, and introns each with to account for the censoring by reweighting each error by PLS, SPLS, lasso, and elnet are considered, along with the the inverse probability of censoring [13]. This is given by, clinical data only model. These 13 models are combined 1/2 n 2 using the ensemble method presented in [12], which is δ y −ˆy i i i RMSE = ,(4) briefly summarized here. n ˆ S T − i=1 i For i = 1, ... , B iterations, do the following where n isthesamplesizeofthe testingdataset, δ is 1 1 From the original training dataset, resample n if the ith patient is uncensored and 0 otherwise, y is the observations with replacement. This set is referred to observed survival time for patient i, y ˆ is the predicted as the bag and will be used to train the ensemble. The survival time, and S is the survival function of censor- out of bag (OOB) samples consist of those not ing. Note that S can be estimated by the Kaplan-Meier chosen for the bag and are used to test the ensemble. estimator with δ replaced by 1 − δ. 2Eachofthe M = 13 models are fit on the bag samples. A reviewer suggested Harrell’s c-index as an alternative 3 Compute K performance measures for each model measure to RMSE. The c-index measures the concordance using the OOB samples. of predicted survival times with true survival times. It is 4 The models are ordered R ,for j = 1, ... , M,by (j) computed as rank aggregation of the K measures. The best model δ I y ˆ < y ˆ I y < y R is collected. i i j i j i=j (1) = .(5) δ I(y < y ) i i j i=j This process results in a collection of B models. The ensemble method uses the average of the predicted sur- In contrast to RMSE, the c-index only considers the rel- vival times from each of these B models. ative ordering of the predicted times. The c-index ranges In this study, we consider K = 3 different measures: from 0 to 1, with values close to 1 indicating strong the RMSE and two logrank test statistics described below. performance. A total of B = 20 iterations are conducted, which keeps The final two measures are based on the LPS classifica- the computational burden at a minimum while maintain- tion of patients using cutoffs t = 2and t = 5. A model is ing desirable results. In addition, to avoid repeating the considered to peform well if it is able to separate patients imputation procedure for each model at each iteration, the into two groups having distinctly different survival curves. Grimes et al. Biology Direct (2018) 13:11 Page 4 of 15 To measure this property, the logrank test [14]isused, transcripts, and 340,414 introns, respectively. A hierarchi- which compares the estimated survival curves for each cal version of the transcript annotation was also available group (LPS versus non-LPS). The test statistic is given by butwas notused. Normalization of the RNA-seq data was performed by O − E g g [16]. The gene counts were normalized as the log of the ,(6) Var O − E g g number of bases aligned in the gene, divided by the num- ber of terabases aligned in known genes and by the length where O − E = d − d (n /n ) is the sum g g g,f f g,f f f ∈F of the gene, with several corrections. The same normal- of observed minus expected deaths in group g = 1, 2, ization is used for the transcript counts. The expressions where F is the set of all observed survival times, d is g,f for the introns are computed as the number of deaths in group g at time f, n is the g,f number of patients at risk in group g at time f,and n (1 + number of supporting reads) ∗ 10 log . is the total number at risk at time f. The survdiff func- number of reads supporting an intron in this data tion in the “survival” R package [15]isusedtocompute The RNA-seq data are filtered prior to model fitting. this statistic. Under the null hypothesis of no difference Genes and transcripts without an NCBI ID are removed. between survival curves, the logrank test statistic has an Any variables with over 80% zero counts in the training asymptotically χ distribution with 1 degree of freedom. dataset are also omitted. A database of 3681 genes related The performance measures for each model are shown in to neuroblastoma was obtained from the GeneCards Suite Figs. 1 and 2. For RMSE and the logrank tests, smaller val- [18]. This dataset is used to subset the remaining genes ues correspond to better performance. For c-index, values and transcripts, resulting in 3389 genes and 47276 tran- close to 1 are better. The error bars are 95% confidence scripts. For the introns, their predictive ability for survival intervals obtained by bootstraping on the testing dataset; time is ranked by fitting each intron in a Cox proportional observations are resampled with replacement and each hazards model [19, 20]. This is repeated for both OS and measure is recomputed. The process is repeated B = 1000 EFS times of patients in the training set. The Cox model times. The 2.5th and 97.5th percentiles are used for the is fit using the “survival” R package [15]. The top 2000 lower and upper confidence limits, respectively. introns with the smallest p-values (testing that the coeffi- cient is zero) are used. This ranking is also performed on Datasets the remaining genes and transcripts; the top 2,000 of each The datasets can be accessed from the GEO database are retained. with accession number GSE49711 [16, 17].Thedataare comprised of tumor samples from 498 neuroblastoma Results patients from seven countries: Belgium (n = 1),Germany Eighteen models are considered in total. Each model is (n = 420),Israel (n = 11),Italy (n = 5),Spain (n = 14), used to estimate overall survival (OS) and event-free sur- United Kingdom (n = 5),andUnited States (n = 42). vival (EFS). For a baseline of comparison, a “null” model Several clinical variables are available for each patient, is fit using clinical covariates alone. Models are then con- along with the RNA-sequencing information from their structed by first selecting a set of predictors: genes, tran- tumor sample. In [16], thedatawererandomlysepa- scripts, introns, or both transcripts and introns (labeled rated into a training set and testing set; this partition was G, T, I, and TI, respectively); and then choosing one of the recorded with the clinical data and is used here. four dimension reduction techniques: PLS, SPLS, lasso, or elastic net (labeled 1-4, respectively). This gives 16 Clinical data possible combinations. Finally we consider an ensemble The clinical data consist of 11 variables. In this study, three model, which pools together the null model and individual of these variables are used as clinical covariates: sex, age, models containing genes, transcripts, or introns. and MYCN status. There are two outcomes of interest: overall survival and Predicting survival times directly event-free survival. Overall survival is calculated as the The models using RNA-seq data tend to perform bet- time from diagnosis to the time of death from disease or ter than the null model in predicting survival times. A the last follow-up date, if the patient survived. Event-free 95% confidence interval (CI) for the adjusted root mean survival is calculated as the time from diagnosis to the squared error (RMSE) of each model is estimated via boot- time of tumor progression, relapse, or death from disease, strapping on the testing set; these are shown in Figs. 1 or to the last follow-up date if no event occurred. and 2. RNA-seq data ForOS,theestimated95% CIforRMSEofthe null The RNA-seq data provide annotations at three feature model is (2.66, 7.61). Every other model besides G-1, G-3, levels, giving datasets comprised of 60,776 genes, 263,544 and G-4 (genes using PLS, lasso, and elnet, respectively) Grimes et al. Biology Direct (2018) 13:11 Page 5 of 15 Fig. 1 Performance measures for overall survival. Each of the 18 models are assessed using the testing dataset. Four measures of performance are considered: the adjusted root mean squared prediction error (RMSE); the logrank test statistic from using the predicted survival time as a classifieron high-risk patients, thresholded at 2 years (LPS2) and 5 years (LPS5); and Harrell’s c-index. 95% confidence intervals are obtained by bootstraping on the testing dataset. This is done by resampling observations with replacement and recomputing each measure. The process is repeated for B = 1000 times, and the middle 95% of measures are used for the confidence interval have smaller RMSE estimates than the null model. How- which overlaps slightly with the null model’s. The I-1 and ever, only the TI-2 model (transcripts and introns using I-2 models (introns using PLS and SPLS) have confidence SPLS) has a confidence interval bounded below the null intervals bounded below the null model’s (Fig. 7). model’s, with an estimated 95% CI of (1.23, 2.60) (Fig. 6). Overall, the performance of predicting exact survival For EFS, the improvements of the RNA-seq models over times is not completely satisfactory. For a patient with the null model appear to be less substantial. The estimated high predicted survival, say 20 years or more, an RMSE 95% CI for RMSE of the null model is (4.37, 5.52). Only five of 1-2 years is acceptable; we can reliably conclude that of the 16 RNA-seq models have lower RMSE estimates this is a low-risk patient who won’t require intensive treat- than the null model. The TI-2 model still performed well ment. However, a clinically high-risk patient may have a in comparison with a 95% CI for RMSE of (2.02, 4.49), predicted survival time of 5 years or less, in which case an Grimes et al. Biology Direct (2018) 13:11 Page 6 of 15 Fig. 2 Performance measures for event-free survival. Each of the 18 models are assessed using the testing dataset. Four measures of performance are considered: the adjusted root mean squared prediction error (RMSE); the logrank test statistic from using the predicted survival time as a classifier on high-risk patients, thresholded at 2 years (LPS2) and 5 years (LPS5); and Harrell’s c-index. 95% confidence intervals are obtained by bootstraping on the testing dataset. This is done by resampling observations with replacement and recomputing each measure. The process is repeated for B = 1000 times, and the middle 95% of measures are used for the confidence interval. Note, the upper limit of RMSE for T-2 is not visible in the plot RMSE of 1-2 years is troublesome; it is unclear whether or Classification of high-risk patients not an agressive course of treatment should be used. These models can be used as a classifier by compar- A reviewer suggested the use of Harrell’s c-index as ing the predicted survival times to a chosen threshold. an alternative measure to RMSE. This measure consid- Since the clinically high-risk group is notorious for hav- ers the relative ordering of predicted survival times with ing poor prognosis, our goal is focused on subclassifying the observed times [21]. We find that models provide pre- these patients. A threshold of 2 years is used. If a patient dicted times that are strongly concordant with observed has a predicted survival time less than 2 years, they are times (Figs. 1 and 2), which indicates an accurate rela- labeled as LPS (low predicted survival). Otherwise, they tive ordering of patients. These results suggests that the are non-LPS. A classifier is considered successful if the models may be useful as a classifier. two resulting groups (LPS versus non-LPS) have distinct Grimes et al. Biology Direct (2018) 13:11 Page 7 of 15 survival curves. The Kaplan-Meier estimates [22]ofthese between LPS and non-LPS groups than is found with OS curves for each RNA-seq model are shown in Figs. 3, 4, (Figs. 3, 4, 5 and 6). The T-1 model provides the largest 5 and 6, and the null model and ensemble are shown in distinction in 2-year EFS rates: 0.29 ± 0.06 versus 0.56 ± Fig. 7. 0.10 (Table 1). Using OS as the outcome, almost every RNA-seq model In general, subclassification is more successful with OS is able to partition high-risk patients into two distinct than with EFS. The ensemble approach (Fig. 7)reflectsthe groups, providing a substantial improvement over the null overall performance in both cases: the LPS and non-LPS model. The TI-4 model produces groups with the largest groups are well separated by the ensemble in OS (0.42 ± difference in 2-year OS rates: 0.40±0.11 versus 0.80±0.05 0.10 versus 0.82 ±0.05) but not for EFS (0.36 ±0.06 versus (Table 1). With EFS as the outcome, there is less separation 0.39 ± 0.09) (Table 1). Fig. 3 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the gene annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test Grimes et al. Biology Direct (2018) 13:11 Page 8 of 15 Fig. 4 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the transcripts annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test Pathway analysis fit using elnet are considered, as these contain an amount Pathway enrichment analysis provides a biological sum- of sparsity that is appropriate for pathway analysis. Two mary of the genes selected by the AFT model. Gene sets gene sets are constructed, one associated with OS and the are constructed by collecting the predictors with nonzero other with EFS. Pathway enrichment analysis (on KEGG coefficients in the fitted G-4, T-4 and TI-4 models. The I- pathways) is performed using DAVID 6.8 [23] and sum- 4 model with introns only is not considered, since introns marized in Tables 2 and 3. cannot easily be interpreted in the pathway analysis. The When predicting OS, a total of 354 unique genes are PLS and SPLS methods gave each predictor some weight given nonzero coefficients by one of the three models. in the AFT model, while the predictors selected by lasso Of these genes, 186 are annotated in KEGG pathways. are a subset of those selected by elnet. Hence, only models DAVID uses a modified fisher exact test to compute Grimes et al. Biology Direct (2018) 13:11 Page 9 of 15 Fig. 5 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the introns annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test p-values for enrichment, and the Benjamini-Hochberg genes annoted in KEGG pathways. However, the RNA-seq correction is applied to account for multiple testing [24]. data used in this study are filtered based on the GeneCards Two pathways are found to be significantly enriched: Path- database. Hence, the pathway enrichment may be more ways in Cancer and ErbB signaling pathway (Table 2). appropriately conducted using those GeneCard genes For EFS, 246 unique genes have nonzero coefficients, of as the background. The GeneCards database contained which 135 are annoted in KEGG pathways. However, no 3512 genes related to neuroblastoma, of which 2044 pathways are enriched for EFS at the 0.05 significance are annoted in KEGG pathways. Relative to this back- level. ground, three pathways are enriched for OS: ErbB sig- The preceeding enrichment analysis uses the entire naling pathway, Salivary secretion, and Inflammatory human genome as a background, which contains 6910 mediator regulation of TRP channels (Table 3). Five Grimes et al. Biology Direct (2018) 13:11 Page 10 of 15 Fig. 6 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using both the transcript and intron annotations from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test pathways are enriched for EFS: Terpenoid backbone procedure, to predict overall survival (OS) and event-free biosynthesis; Metabolic pathways; Valine, leucine and survival (EFS) times of neuroblastoma patients. Three fea- isoleucine degradation; Biosynthesis of antibiotics; and ture levels of an RNA-seq dataset were considered, includ- Fatty acid metabolism (Table 3). These pathways have p- ing genes, transcripts, and introns. Models were fit using valuesbelowthe0.05significancelevel,butarenonsignif- the three features independently and with transcripts and icant after applying the Benjamini-Hochberg correction. introns together. In terms of RMSE, the predictive performance of OS Discussion is greatly improved in the RNA-seq models over the null In this study we used the AFT model, fit using various model, but this improvement is curtailed when predict- dimension reduction techniques and a dataset imputation ing EFS. The high rate of censoring that is found in this Grimes et al. Biology Direct (2018) 13:11 Page 11 of 15 Fig. 7 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the null model (first row) and the ensemble approach (second row). The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test data will be a hinderance for any nonparametric model. OS, two cancer-related pathways are enriched. This anal- Alternative approaches can be considered: One possibil- ysis may be biased, however, since the RNA-seq data are ity is to switch to semiparametric estimation, but this initially filtered using the GeneCards database. If the back- ground is altered to reflect this filtering, we find that approach will be computationally intensive in this high- dimensional setting. A more practical solution may be to one of the two cancer-related pathways remains relatively employ a boosting algorithm (see [25] for example). These enriched. This alteration also reveals additional enriched alternatives were not explored in detail in this paper. pathways for the OS and EFS gene sets, but their relevance The second goal is to subclassify clinically high-risk to neuroblastoma is questionable. Since the prediction (HR) patients. In this venture the AFT model produces of EFS had limited success, it is no surprise that the very promising results. High-risk patients with low sur- genes selected for EFS appear to have limited biological vival times are more sensitive to the amount of error relevance. remaining in predicted times, but the estimates tend to The predictive accuracy and pathway enrichment for OS be in the right direction. That is, the relative ordering suggests that the AFT model with elastic net is able to pick of the patients by their predicted survival times is accu- out biologically meaningful genes. A future study pursu- rate. A reviewer suggested the use of Harrell’s c-index ing this kind of interpretation will need to consider the [21] to measure this effect. The c-index is above 0.8 for stochastic nature of the fitting procedure and determine a each model when predicting OS, indicating strong con- stable set of genes selected by the model. As suggested by a cordance between predicted OS time and true OS times reviewer, we can also explore relationships between these (Fig. 1). The concordance is less strong when predicting genes and those excluded by the initial filtering process. EFS (Fig. 2). Such an investigation may produce biological insights into Using a cutoff of 2 years, each model is converted to a the subgroups of high-risk patients. classifier. The TI-4 model provides the best results for OS. An ensemble of models was considered, which incorpo- For EFS, the I-4 model appears to be the best. A classi- rates bagging with rank aggregation of three performance fier using 5 years as a cutoff is also considered, but the measures. The performance of the ensemble method is performance is not as good; setting the threshold to a comparable to that of the best individual model. This suggests that the ensemble method is able to effectively value below 5 years seems to be necessary in order to iden- tify those patients who are at the highest risk in the HR combine models fit on separate datasets. If additional group. datasets are incorporated, such as copy number varia- A pathway analysis of the gene sets selected by the elas- tion or other -omics data, the AFT model can be fit tic net when predicting OS and EFS is performed. With by simply concatenating the datasets together, but the Grimes et al. Biology Direct (2018) 13:11 Page 12 of 15 Table 1 Summary of Kaplan-Meier estimates for 2-year OS and 2-year EFS for clinically high-risk patients using each of the 18 proposed models LPS non-LPS Outcome Prob SE N Prob SE N Data Model P-value OS 0.57 0.09 31 0.79 0.05 59 Null model 0.574 OS 0.60 0.07 57 0.91 0.05 33 G 1 0.081 OS 0.60 0.06 60 0.93 0.05 30 G 2 0.082 OS 0.55 0.09 29 0.79 0.05 61 G 3 0.166 OS 0.50 0.10 24 0.79 0.05 66 G 4 0.164 OS 0.57 0.07 55 0.94 0.04 35 T 1 0.128 OS 0.56 0.07 47 0.88 0.05 43 T 2 0.002 OS 0.41 0.12 17 0.78 0.05 73 T 3 0.048 OS 0.57 0.09 33 0.79 0.05 57 T 4 0.081 OS 0.61 0.07 56 0.88 0.06 34 I 1 0.124 OS 0.51 0.08 36 0.85 0.05 54 I 2 0.027 OS 0.47 0.10 25 0.81 0.05 65 I 3 0.137 OS 0.43 0.10 25 0.82 0.05 65 I 4 0.017 OS 0.57 0.07 56 0.94 0.04 34 TI 1 0.092 OS 0.54 0.07 49 0.90 0.05 41 TI 2 0.001 OS 0.40 0.11 21 0.81 0.05 69 TI 3 0.036 OS 0.40 0.11 22 0.80 0.05 68 TI 4 0.009 OS 0.42 0.10 25 0.82 0.05 65 Ensemble 0.016 EFS 0.39 0.08 37 0.36 0.07 53 Null model 0.616 EFS 0.35 0.07 53 0.39 0.08 37 G 1 0.299 EFS 0.34 0.07 52 0.41 0.08 38 G 2 0.266 EFS 0.30 0.06 58 0.50 0.09 32 G 3 0.125 EFS 0.31 0.06 59 0.48 0.09 31 G 4 0.125 EFS 0.29 0.06 62 0.56 0.10 28 T 1 0.036 EFS 0.35 0.07 51 0.40 0.08 39 T 2 0.946 EFS 0.29 0.07 38 0.43 0.07 52 T 3 0.065 EFS 0.31 0.07 47 0.43 0.08 43 T 4 0.150 EFS 0.35 0.06 66 0.43 0.10 24 I 1 0.256 EFS 0.35 0.06 66 0.43 0.10 24 I 2 0.256 EFS 0.29 0.07 43 0.45 0.07 47 I 3 0.052 EFS 0.29 0.07 42 0.44 0.07 48 I 4 0.082 EFS 0.31 0.06 57 0.47 0.09 33 TI 1 0.084 EFS 0.33 0.06 68 0.50 0.11 22 TI 2 0.183 EFS 0.29 0.07 42 0.45 0.07 48 TI 3 0.062 EFS 0.30 0.07 44 0.44 0.08 46 TI 4 0.085 EFS 0.36 0.06 58 0.39 0.09 32 Ensemble 0.599 Patients with predicted survival of less than 2 years are labeled as Low Predicted Survival (LPS), and otherwise are non-LPS. Columns labeled “Prob.”, “SE”, and “N” correspond to the estimated probability of 2-year survival, the standard error of the estimate, and the number of patients in the given cohort. The P-values are for the logrank test comparing LPS to non-LPS survival. The “Data” column refers to the type of RNA-seq data used, and the “Model” column refers to the dimension reduction technique used computational requirement quickly becomes too bur- that this heuristic works well in combining different anno- densome. The ensemble approach may provide a useful tations of RNA-seq data, but further investigation is heuristic for combining several datasets. We have shown needed to verify the performance with disparate datasets. Grimes et al. Biology Direct (2018) 13:11 Page 13 of 15 Table 2 Pathway enrichment analysis of genes selected by the is a very well-written paper. Overall, the analysis is scien- G-4, T-4, and TI-4 models when predicting OS (no pathways were tifically compelling and relies on creative applications of significantly enriched for EFS) sound statistical techniques. The classifier comparing the Outcome Pathway Count Size P-value BH predicted survival times to a 2-year threshold is success- OS Pathways in cancer 26 393 < 0.001 0.010 ful when it is based on transcript and intron annotations. The ensemble method and its potential application to fit- OS ErbB signaling pathway 11 87 < 0.001 0.012 ting disparate datasets holds much promise for future The columns include KEGG pathway name, the number of genes in the gene set work. that are in the pathway, the total number of genes annotated for the pathway, the p-value from a modified fisher’s exact test, and the Benjamini-Hochberg corrected Reviewer comment: As a suggestion for future p-value research, but entirely unrelated to the current paper which is more than satisfactory, I have the following suggestion. Conclusion From the second paragraph of the Discussion, it appears In this study, we explored the performance of the AFT that it may be helpful to explore Harrell’s C-index as an model in predicting survival times for neuroblastoma alternative measure of accuracy. This may be a better mea- patients. A classifier was constructed by comparing pre- sure than RMSE for the parametric models, especially dicted survival times to a 2-year threshold. Using both because they appear to get the relative ordering of the transcript and intron annotations in the model gave the survival times right rather than the actual magnitudes. best performance. We are able to subclassify clinically Author’s response: We thank Dr. Guha for this sug- high-risk patients into two distinct groups, one with a 40% gestion. The performance of each model using Harrell’s 2-year overall survival rate and the other at 80%. This c-index has been added to the revised manuscript. suggests that the AFT model is useful in subclassifying Reviewer comment: On Line 7 of page 2, should the high-risk patients, which can help clinicians in choosing comma following INSS be deleted? 2. On Line 7 of page 6, effective treatment plans. Only RNA-seq data was con- what is K? sidered in this study, but other types of data can be used Author’s response: Grammatical corrections have been as well. The ensemble method is a useful heuristic for made to the manuscript. For the latter point, there are K = combining several high-dimensional datasets under this 3 performance measures in this study. This is now clarified framework, and it has been shown capable of maintaining in the text. optimal performance. Reviewers’ comments Reviewer’s report 2: Isabel Nepomuceno, Universidad de Reviewer’s report 1: Subharup Guha, University of Florida, Sevilla, Seville, Spain Gainesville, USA In this paper, authors used the accelerated failure time (AFT) model with four dimension reduction techniques The authors explore the performance of the AFT model in and a dataset imputation scheme to predict overall sur- predicting survival times for neuroblastoma patients. This vival and event-free survival times of neuroblastoma Table 3 Pathway enrichment analysis of genes selected by the patients. Three feature levels of and RNA-Seq dataset G-4, T-4, and TI-4 models were considered. Authors shown that the use of RNA-Seq Outcome Pathway Count Size P-value BH data improves accuracy in comparison to using clinical data alone. In general the paper is appropriate to the OS ErbB signaling pathway 11 60 0.029 0.999 journal. The analysis presented in this paper is very inter- OS Salivary secretion 6 23 0.042 0.995 esting. I have several suggestions and comments to be OS Inflammatory mediator 9 48 0.049 0.983 revised: regulation of TRP channels Reviewer comment: The Method section is written in a clear manner but is difficult to reproduce. Authors EFS Terpenoid backbone 4 8 0.010 0.906 biosynthesis mentioned the R package used but they don’t provide the EFS Metabolic pathways 29 304 0.016 0.847 Rcodeofthe study. Author’s response: We thank Dr. Nepomuceno for her EFS Valine, leucine and 5 20 0.032 0.911 isoleucine degradation comments and suggestions. All R code and output is available from GitHub at https://github.com/tgrimes/ EFS Biosynthesis of 12 98 0.037 0.882 antibiotics CAMDA-2017-Neuroblastoma. The session info is also EFS Fatty acid metabolism 5 21 0.037 0.820 reported, which includes the R version, computer specifica- tions, and a list of the packages used during the analysis. In this analysis, the GeneCards genes are used at the background. The columns include survival outcome (OS or EFS), KEGG pathway name, the number of genes in Reviewer comment: The Ensemble Method subsec- the gene set that are in the pathway, the total number of genes annotated for the tion, authors use bagging with rank aggregation over each pathway, the p-value from a modified fisher’s exact test, and the performance measure and set B to 20. Why this parameter Benjamini-Hochberg corrected p-value Grimes et al. Biology Direct (2018) 13:11 Page 14 of 15 is fixed to 20 should be explained. And authors should the GeneCards database to subset on genes, which would explain why the use bagging instead of cross validation. bias our selection to genes in cancer-related pathways. In Author’s response: The choice of 20 iterations for bag- response, we have modified this section and now conduct ging is a compromise between computation time and model a pathway enrichment analysis. However, a question is performance. We also considered B = 50 but did not find a raised regarding the choice of background: should our gene sets be compared to all genes in the genome (as is usu- substantial change in performance. Reviewer comment: The description of the RNA-Seq ally done) or to the GeneCards genes that we subset on? Data, authors reduce the "raw data" with 60776 genes With the former, there is a concern that the analysis may be into 3401 using the 3681 genes related to neuroblastoma biased. Results for both of these scenarios have been added obtained from the Gene Cards Suite. Have authors made to the manuscript. some analysis from the remaining genes? Could be genes Reviewer comment: Finally, as minor comments: - The related with the problem and not related with the disease? Bibliography Section must be revised, there are some It could be interesting to do a cluster analysis to see if the incomplete reference as for example number 14. - In grouped genes using prior knowledge are also clustered Table 1,one ofthemodelsisnamedsimpleforthe baseline together in this analysis. model. It should be names null model as authors explained Author’s response: These are interesting suggestions that before. deserve a separate analysis to be fully addressed. The Author’s response: The bibliography section has been main purpose in using the Gene Cards database was to corrected, and the tables and figures have been relabeled provide an initial filtering to speed up computation. We to be consistent with the text. also re-ran the analysis without this step and found lit- tle difference in predictive performance. We are careful Abbreviations not to place too much emphasis on the interpretation of AFT: Accelerated failure time; CI: Confidence interval; EFS: Event-free survival; elnet: Elastic net; HR: High-risk; INSS: International neuroblastoma staging the gene sets obtained in this analysis. As you’ve pointed system; lasso: Least absolute shrinkage and selection operator; LPS: Low out, there are many new questions that have been uncov- predicted survival; OS: Overall survival; PLS: Partial least squares; RMSE: Root ered and deserve careful consideration. We’ve added some mean squared error; SPLS: Sparse partial least squares comments regarding this in the discussion section of the Acknowledgements manuscript. We acknowledge the CAMDA 2017 committee reviewers for their valuable Reviewer comment: Furthermore, a reference about comments. the Cox proportional hazards model or the R package used Funding should be added. Publication of this article was funded in part by the University of Florida Open Author’s response: We thank the author for pointing Access Publishing Fund. out this omission. The revised manuscript now contains Availability of data and materials additional references. The datasets used in this article be accessed from the GEO repository with Reviewer comment: Section Results, classification of series accession number GSE49711. All R code is available on GitHub at https://github.com/tgrimes/CAMDA-2017-Neuroblastoma. high-risk patients should be rewritten. The second and third paragraph is confused and difficult to see which plot Authors’ contributions corresponds with each sentence. All authors took part in research discussion and the final manuscript preparation. TG performed the data analysis and wrote the first draft of the Author’s response: This section has been reworded to paper. SoD planned the study and the selection of statistical methods. All clarify which table or figure each sentence is referring to. authors read and approved the final manuscript. The titles for each plot have been changed in concordance Ethics approval and consent to participate to the labels used to identify each model within the Not applicable. manuscript. Reviewer comment: In section Pathway analysis, Competing interests The authors declare that they have no competing interests. authors claim that several genes are involved in several pathways.Thatmeans,dogenes appear inthepathways Publisher’s Note or are the pathways enriched by the set of genes? If it is Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. the second case, authors should add a table with the list of pathways, the number of entities in the pathways and Received: 12 October 2017 Accepted: 1 May 2018 the number of genes from the set which appear in the pathway. Author’s response: We thank the reviewer for prompt- References 1. Bosse KR, Maris JM. Advances in the translational genomics of ing this clarification. Previously, the interpretation was neuroblastoma: From improving risk stratification and revealing novel that genes appear in the pathways. But this initial biology to identifying actionable genomic alterations. Cancer. approach seems uninformative, particularly since we use 2016;122(1):20–33. Grimes et al. Biology Direct (2018) 13:11 Page 15 of 15 2. Brodeur GM, Seeger RC, Schwab M, Varmus HE, Bishop JM. Amplification of n-myc in untreated human neuroblastomas correlates with advanced disease stage. Science. 1984;224(4653):1121–4. 3. Formicola D, Petrosino G, Lasorsa VA, Pignataro P, Cimmino F, Vetrella S, Longo L, Tonini GP, Oberthuer A, Iolascon A, et al. An 18 gene expression-based score classifier predicts the clinical outcome in stage 4 neuroblastoma. J Transl Med. 2016;14(1):142. 4. Tan Q, Thomassen M, Jochumsen KM, Mogensen O, Christensen K, Kruse TA. Gene selection for predicting survival outcomes of cancer patients in microarray studies. Adv Comput Inf Sci Eng. 2008;1(1):405–9. 5. Boulesteix A-L, Strimmer K. Partial least squares: a versatile tool for the analysis of high-dimensional genomic data. Brief Bioinform. 2006;8(1): 32–44. 6. Chung D, Chun H, Keles S. Spls: Sparse Partial Least Squares (SPLS) Regression and Classification. 2018. R package version 2.2-2. https:// CRAN.R-project.org/package=spls. Accessed 28 Apr 2018. 7. Chun H, Keles¸ S. J R Stat Soc Series B (Stat Methodol). 2010;72(1):3–25. 8. Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Series B (Stat Methodol). 1996;58(1):267–88. 9. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1. 10. Zou H, Hastie T. J R Stat Soc Series B (Stat Methodol). 2005;67(2):301–20. 11. Mostajabi F, Datta S, Datta S. Predicting patient survival from proteomic profile using mass spectrometry data: an empirical study. Commun Stat Simul Comput. 2013;42(3):485–98. 12. Shah J, Datta S, Datta S. A multi-loss super regression learner (msrl) with application to survival prediction using proteomics. Comput Stat. 2014;29(6):1749–67. 13. Datta S. Estimating the mean life time using right censored data. Stat Methodol. 2005;2(1):65–9. 14. Kleinbaum DG, Klein M. Kaplan-meier survival curves and the log-rank test. In: Survival Analysis, 3rd edn. New York: Springer; 2012. p. 55–96. Chap. 2. 15. Therneau TM. A Package for Survival Analysis in S. 2015. version 2.38. https://CRAN.R-project.org/package=survival. Accessed 28 Apr 2018. 16. Su Z, Fang H, Hong H, Shi L, Zhang W, Zhang W, Zhang Y, Dong Z, Lancashire LJ, Bessarabova M, et al. An investigation of biomarkers derived from legacy microarray data for their utility in the rna-seq era. Genome Biol. 2014;15(12):523. 17. Zhang W, Yu Y, Hertwig F, Thierry-Mieg J, et al. Comparison of rna-seq and microarray-based models for clinical endpoint prediction. Genome Biol. 2015;16(1):133. 18. Safran M, Dalah I, Alexander J, Rosen N, Iny Stein T, Shmoish M, Nativ N, Bahir I, Doniger T, Krug H, et al. Genecards version 3: the human gene integrator. Database. 2010;2010:020. 19. Cox DR. J R Stat Soc Series B (Stat Methodol). 1972;34(2):187–220. 20. Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer; 2000. 21. Harrell Jr FE, Califf RM, Pryor DB, Lee KL, Rosati RA, et al. Evaluating the yield of medical tests. J Am Med Assoc. 1982;247(18):2543–6. 22. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53(282):457–81. 23. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using david bioinformatics resources. Nat Protocol. 2009;4(1):44. 24. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B (Methodol). 1995;57(1):289–300. 25. Schmid M, Hothorn T. Flexible boosting of accelerated failure time models. BMC Bioinformatics. 2008;9(1):269. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Biology Direct Springer Journals

Predicting survival times for neuroblastoma patients using RNA-seq expression profiles

Free
15 pages
Loading next page...
 
/lp/springer_journal/predicting-survival-times-for-neuroblastoma-patients-using-rna-seq-MU12VRGtH4
Publisher
BioMed Central
Copyright
Copyright © 2018 by The Author(s)
Subject
Life Sciences; Life Sciences, general
eISSN
1745-6150
D.O.I.
10.1186/s13062-018-0213-x
Publisher site
See Article on Publisher Site

Abstract

Background: Neuroblastoma is the most common tumor of early childhood and is notorious for its high variability in clinical presentation. Accurate prognosis has remained a challenge for many patients. In this study, expression profiles from RNA-sequencing are used to predict survival times directly. Several models are investigated using various annotation levels of expression profiles (genes, transcripts, and introns), and an ensemble predictor is proposed as a heuristic for combining these different profiles. Results: The use of RNA-seq data is shown to improve accuracy in comparison to using clinical data alone for predicting overall survival times. Furthermore, clinically high-risk patients can be subclassified based on their predicted overall survival times. In this effort, the best performing model was the elastic net using both transcripts and introns together. This model separated patients into two groups with 2-year overall survival rates of 0.40 ± 0.11 (n = 22) versus 0.80 ± 0.05 (n = 68). The ensemble approach gave similar results, with groups 0.42 ± 0.10 (n = 25) versus 0.82 ± 0.05 (n = 65). This suggests that the ensemble is able to effectively combine the individual RNA-seq datasets. Conclusions: Using predicted survival times based on RNA-seq data can provide improved prognosis by subclassifying clinically high-risk neuroblastoma patients. Reviewers: This article was reviewed by Subharup Guha and Isabel Nepomuceno. Keywords: Accelerated failure time, Sparse PLS, Lasso, Elastic net, Data imputation Background been one of the most important markers for stratify- Neuroblastoma is the most frequently diagnosed cancer ing patients. Genome-wide association studies have found in the first year of life and the most common extracra- many other SNPs associated with an increased risk of nial solid tumor in children. It accounts for 5% of all neuroblastoma. However, while aberrations of these genes pediatric cancer diagnoses and 10% of all pediatric oncol- indicate an increased susceptibility to the disease, these ogy deaths [1]. These numbers have improved over the markers are less useful for stratifying patients into risk past decade, but accurate prognosis for the disease has groups after diagnosis. remained a challenge [1]. The difficulty is due to the highly The Children’s Oncology Group stratifies patients into heterogeneous nature of neuroblastoma; cases can range three risk groups using the International Neuroblastoma from tumors that spontaneously regress on their own, to Staging System (INSS) and various prognostic markers aggressive tumors that spread unabated by treatment. including age at diagnosis, tumor histology, MYCN ampli- In 1984, the MYCN oncogene was identified as a fication, and DNA ploidy. According to the American biomarker for clinically aggressive tumors [2]. It has since Cancer Society, the 5-year survival rate for these low-risk, intermediate-risk, and high-risk groups are > 95%, 90% - 95%, and < 50%, respectively. The high-risk group typi- cally consists of patients older than 18 months with INSS *Correspondence: somnath.datta@ufl.edu stage 4 or patients of any age with MYCN amplification. Department of Biostatistics, University of Florida, 2004 Mowry Rd, 32611, Gainesville, USA © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Grimes et al. Biology Direct (2018) 13:11 Page 2 of 15 Predicting survival outcomes using gene expression data PLS is implemented using the “spls” R package [6]. The has been explored with promising results [3, 4]. These number of latent factors v is a tuning parameter, which studies use gene expression profiles with classification is determined from 10-fold cross validation. The optimal methods to stratify patients into risk groups. However, value of v is searched over v = 1, ... , 10. patients that are clinically labeled as high-risk pose a SPLS particular challenge, and classifiers tend to struggle sep- Like PLS, the sparse partial least squares (SPLS) also con- arating those patients into subgroups. In this paper, we structs latent factors, but it incorporates L regularization take the approach of modeling survival time directly using in the process [7]. This induces sparsity in each linear RNA-seq data. This leads to two objectives: the first is to combination of the original covariates that make up the evaluate the accuracy of the model in predicting exact sur- latent factors. There are two tuning parameters, the num- vival times. The second is to determine whether the pre- ber of latent factors v < n and the shrinkage parameter dicted times can be used to subclassify high-risk patients η ∈ (0, 1) for the regularization. Both of these are deter- into distinct groups. mined from 10-fold cross validation using the “spls” R package [6]. The optimal values of v and η are searched Methods over the grid of points with v =1, ..,10 and η = 0.1, ...,0.9. Accelerated failure time (AFT) model Note, to implement PLS the shrinkage parameter, η,is Theacceleratedfailuretime(AFT)modelrelates thelog set to zero. survival times to a linear combination of the predictors. Lasso log(y) = Xβ + ,(1) The least absolute shrinkage and selection operator (lasso) fits the model using least squares subject to an L con- where y ∈ R denotes the vector of n observed sur- straint on the parameters |β |≤ λ,where λ> 0is j=1 vival times, X the n × p matrix with columns containing a tuning parameter that affects the amount of shrinkage the predictor variables for each observation, β ∈ R the [8]. This constraint induces sparsity in the estimated coef- vector of regression coefficients, and  ∈ R avector ficients, setting many coefficients to zero and shrinking of independent random errors with an unspecified dis- others. tribution that is assumed to be independent of X.The The model is fit using the “glmnet” R package [9], which predictors X are centered and scaled so that each col- performs 10-fold cross validation to select λ. umn X ,for i = 1, ... , p, has zero mean and unit variance There are two challenges to fitting this model: the high Elastic net dimensionality of X and the right censoring of y.Since The elastic net (elnet) uses a similar approach as the lasso. It combines both L and L penalties; the estimator p > n, ordinary least squares (OLS) should not be used as 1 2 minimizes the convex function it will simply overfit on the data. Instead, four approaches for dimension reduction are considered, which include 1 1 2 2 ||Y − Xβ|| + λ (1 − α)||β|| + α||β|| ,(2) both latent factor and regularization techniques. To han- 2 2 2 2 dle right censoring, a nonparametric, iterative imputation where λ> 0and α ∈ [0, 1] are two tuning parameters procedure is proposed, which allows the model to be fit as [10]. When α = 1, this reduces to the lasso. By includ- though complete data were available. ing some component of the L penalty, groups of strongly Each of the dimension reduction techniques require the correlated variables tend to be included or excluded in the selection of one or more tuning parameters. These param- model together. The “glmnet” R package [9]isusedtofit eters are determined by 10-fold cross validation, which the model and determine both tuning parameters. is implemented in R using two packages discussed in the following sections. Imputation for right censoring PLS Let (y , δ , X )|i = 1, ... , n denote the set of observed i i i With partialleast squares(PLS),acollectionof v < n survival times, indicators for death from disease, and the orthogonal latent factors are computed as linear combi- p-dimensional vector of covariates for the n patients in the nations of the original covariates. The construction of the dataset. Let T denote thetruesurvivaltimesforpatient latent factors takes into account both X and y;thisisin i = 1, ... , n.Ifthe ith patient’s survival time is censored contrast to principal component analysis (PCA), which (i.e. δ = 0) then we only observe y < T .Thatis, T is i i i i only considers X. A review of PLS and its application to unobserved. genomic data can be found in [5]. Once the v latent factors To deal with this right censoring, the dataset imputa- are computed, the AFT model is fit using OLS with these tion procedure from [11] is used. This procedure is briefly (0) new variables. summarized here. To begin, an initial estimate β is Grimes et al. Biology Direct (2018) 13:11 Page 3 of 15 obtained by fitting the AFT model using only the uncen- censored data is imputed once at the start of the ensemble sored data. Then, in each of k = 1, ... , n iterations, do training; the censored survival times are replaced with the the following. predicted times from the single best model (TI-4). (k) 1 Calculate the Kaplan-Meier estimate S (e) of the Classification: LPS vs. non-LPS { } distribution of model error using (e , δ )|i = 1, ... , n i i The second goal is to subclassify clinically high-risk T (k−1) where e = log(y ) − X β . i i patients. A new dichotomous variable is created to classify 2Impute n new datasets by replacing each censored patients: If the predicted survival time is less than t > 0 T (k−1) ∗ ∗ log(y ) with X β + e ,where e is a sampled i i i years, we say the patient has low predicted survival (LPS). model residual from the conditional distribution Otherwise, the patient is non-LPS. For patient i = 1, ..., n (k) S (e|e > e ). This condition ensures that the with predicted survival time y ˆ ,let imputed observation will be larger than the observed 1if y ˆ ≤ t right-censored time. i LPS = .(3) i,t 0otherwise 3 Use the new datasets to compute n new estimates (k) β for j = 1, ... , n . Two cutoffs were considered with t = 2and t = 5years. 4 Average the n estimates to obtain a final estimate For clinically high-risk patients, the t = 2cutoffisuse- n (k) 1 D (k) ˆ ˜ β = β . ful for identifying those with a significantly lower survival j=1 j rate. In the general population of neuroblastoma patients, The process is repeated for n iterations, and the final the t = 5 cutoff is useful as an alternative way to iden- (n ) ˆ K estimate β is returned. tify high-risk patients, but it cannot tease out the more To balance between computation time and simulation extreme cases. variability, we chose to run n = 5 iterations, imputing n = 5datasetsineach. Performance measures Performance is evaluated on the testing dataset by four Ensemble method different measures. The ensemble method incorporates bagging with rank The first involves the prediction error of survival times. aggregation over each performance measure. The 12 This is measured by the root mean squared error, adjusted models using genes, transcripts, and introns each with to account for the censoring by reweighting each error by PLS, SPLS, lasso, and elnet are considered, along with the the inverse probability of censoring [13]. This is given by, clinical data only model. These 13 models are combined 1/2 n 2 using the ensemble method presented in [12], which is δ y −ˆy i i i RMSE = ,(4) briefly summarized here. n ˆ S T − i=1 i For i = 1, ... , B iterations, do the following where n isthesamplesizeofthe testingdataset, δ is 1 1 From the original training dataset, resample n if the ith patient is uncensored and 0 otherwise, y is the observations with replacement. This set is referred to observed survival time for patient i, y ˆ is the predicted as the bag and will be used to train the ensemble. The survival time, and S is the survival function of censor- out of bag (OOB) samples consist of those not ing. Note that S can be estimated by the Kaplan-Meier chosen for the bag and are used to test the ensemble. estimator with δ replaced by 1 − δ. 2Eachofthe M = 13 models are fit on the bag samples. A reviewer suggested Harrell’s c-index as an alternative 3 Compute K performance measures for each model measure to RMSE. The c-index measures the concordance using the OOB samples. of predicted survival times with true survival times. It is 4 The models are ordered R ,for j = 1, ... , M,by (j) computed as rank aggregation of the K measures. The best model δ I y ˆ < y ˆ I y < y R is collected. i i j i j i=j (1) = .(5) δ I(y < y ) i i j i=j This process results in a collection of B models. The ensemble method uses the average of the predicted sur- In contrast to RMSE, the c-index only considers the rel- vival times from each of these B models. ative ordering of the predicted times. The c-index ranges In this study, we consider K = 3 different measures: from 0 to 1, with values close to 1 indicating strong the RMSE and two logrank test statistics described below. performance. A total of B = 20 iterations are conducted, which keeps The final two measures are based on the LPS classifica- the computational burden at a minimum while maintain- tion of patients using cutoffs t = 2and t = 5. A model is ing desirable results. In addition, to avoid repeating the considered to peform well if it is able to separate patients imputation procedure for each model at each iteration, the into two groups having distinctly different survival curves. Grimes et al. Biology Direct (2018) 13:11 Page 4 of 15 To measure this property, the logrank test [14]isused, transcripts, and 340,414 introns, respectively. A hierarchi- which compares the estimated survival curves for each cal version of the transcript annotation was also available group (LPS versus non-LPS). The test statistic is given by butwas notused. Normalization of the RNA-seq data was performed by O − E g g [16]. The gene counts were normalized as the log of the ,(6) Var O − E g g number of bases aligned in the gene, divided by the num- ber of terabases aligned in known genes and by the length where O − E = d − d (n /n ) is the sum g g g,f f g,f f f ∈F of the gene, with several corrections. The same normal- of observed minus expected deaths in group g = 1, 2, ization is used for the transcript counts. The expressions where F is the set of all observed survival times, d is g,f for the introns are computed as the number of deaths in group g at time f, n is the g,f number of patients at risk in group g at time f,and n (1 + number of supporting reads) ∗ 10 log . is the total number at risk at time f. The survdiff func- number of reads supporting an intron in this data tion in the “survival” R package [15]isusedtocompute The RNA-seq data are filtered prior to model fitting. this statistic. Under the null hypothesis of no difference Genes and transcripts without an NCBI ID are removed. between survival curves, the logrank test statistic has an Any variables with over 80% zero counts in the training asymptotically χ distribution with 1 degree of freedom. dataset are also omitted. A database of 3681 genes related The performance measures for each model are shown in to neuroblastoma was obtained from the GeneCards Suite Figs. 1 and 2. For RMSE and the logrank tests, smaller val- [18]. This dataset is used to subset the remaining genes ues correspond to better performance. For c-index, values and transcripts, resulting in 3389 genes and 47276 tran- close to 1 are better. The error bars are 95% confidence scripts. For the introns, their predictive ability for survival intervals obtained by bootstraping on the testing dataset; time is ranked by fitting each intron in a Cox proportional observations are resampled with replacement and each hazards model [19, 20]. This is repeated for both OS and measure is recomputed. The process is repeated B = 1000 EFS times of patients in the training set. The Cox model times. The 2.5th and 97.5th percentiles are used for the is fit using the “survival” R package [15]. The top 2000 lower and upper confidence limits, respectively. introns with the smallest p-values (testing that the coeffi- cient is zero) are used. This ranking is also performed on Datasets the remaining genes and transcripts; the top 2,000 of each The datasets can be accessed from the GEO database are retained. with accession number GSE49711 [16, 17].Thedataare comprised of tumor samples from 498 neuroblastoma Results patients from seven countries: Belgium (n = 1),Germany Eighteen models are considered in total. Each model is (n = 420),Israel (n = 11),Italy (n = 5),Spain (n = 14), used to estimate overall survival (OS) and event-free sur- United Kingdom (n = 5),andUnited States (n = 42). vival (EFS). For a baseline of comparison, a “null” model Several clinical variables are available for each patient, is fit using clinical covariates alone. Models are then con- along with the RNA-sequencing information from their structed by first selecting a set of predictors: genes, tran- tumor sample. In [16], thedatawererandomlysepa- scripts, introns, or both transcripts and introns (labeled rated into a training set and testing set; this partition was G, T, I, and TI, respectively); and then choosing one of the recorded with the clinical data and is used here. four dimension reduction techniques: PLS, SPLS, lasso, or elastic net (labeled 1-4, respectively). This gives 16 Clinical data possible combinations. Finally we consider an ensemble The clinical data consist of 11 variables. In this study, three model, which pools together the null model and individual of these variables are used as clinical covariates: sex, age, models containing genes, transcripts, or introns. and MYCN status. There are two outcomes of interest: overall survival and Predicting survival times directly event-free survival. Overall survival is calculated as the The models using RNA-seq data tend to perform bet- time from diagnosis to the time of death from disease or ter than the null model in predicting survival times. A the last follow-up date, if the patient survived. Event-free 95% confidence interval (CI) for the adjusted root mean survival is calculated as the time from diagnosis to the squared error (RMSE) of each model is estimated via boot- time of tumor progression, relapse, or death from disease, strapping on the testing set; these are shown in Figs. 1 or to the last follow-up date if no event occurred. and 2. RNA-seq data ForOS,theestimated95% CIforRMSEofthe null The RNA-seq data provide annotations at three feature model is (2.66, 7.61). Every other model besides G-1, G-3, levels, giving datasets comprised of 60,776 genes, 263,544 and G-4 (genes using PLS, lasso, and elnet, respectively) Grimes et al. Biology Direct (2018) 13:11 Page 5 of 15 Fig. 1 Performance measures for overall survival. Each of the 18 models are assessed using the testing dataset. Four measures of performance are considered: the adjusted root mean squared prediction error (RMSE); the logrank test statistic from using the predicted survival time as a classifieron high-risk patients, thresholded at 2 years (LPS2) and 5 years (LPS5); and Harrell’s c-index. 95% confidence intervals are obtained by bootstraping on the testing dataset. This is done by resampling observations with replacement and recomputing each measure. The process is repeated for B = 1000 times, and the middle 95% of measures are used for the confidence interval have smaller RMSE estimates than the null model. How- which overlaps slightly with the null model’s. The I-1 and ever, only the TI-2 model (transcripts and introns using I-2 models (introns using PLS and SPLS) have confidence SPLS) has a confidence interval bounded below the null intervals bounded below the null model’s (Fig. 7). model’s, with an estimated 95% CI of (1.23, 2.60) (Fig. 6). Overall, the performance of predicting exact survival For EFS, the improvements of the RNA-seq models over times is not completely satisfactory. For a patient with the null model appear to be less substantial. The estimated high predicted survival, say 20 years or more, an RMSE 95% CI for RMSE of the null model is (4.37, 5.52). Only five of 1-2 years is acceptable; we can reliably conclude that of the 16 RNA-seq models have lower RMSE estimates this is a low-risk patient who won’t require intensive treat- than the null model. The TI-2 model still performed well ment. However, a clinically high-risk patient may have a in comparison with a 95% CI for RMSE of (2.02, 4.49), predicted survival time of 5 years or less, in which case an Grimes et al. Biology Direct (2018) 13:11 Page 6 of 15 Fig. 2 Performance measures for event-free survival. Each of the 18 models are assessed using the testing dataset. Four measures of performance are considered: the adjusted root mean squared prediction error (RMSE); the logrank test statistic from using the predicted survival time as a classifier on high-risk patients, thresholded at 2 years (LPS2) and 5 years (LPS5); and Harrell’s c-index. 95% confidence intervals are obtained by bootstraping on the testing dataset. This is done by resampling observations with replacement and recomputing each measure. The process is repeated for B = 1000 times, and the middle 95% of measures are used for the confidence interval. Note, the upper limit of RMSE for T-2 is not visible in the plot RMSE of 1-2 years is troublesome; it is unclear whether or Classification of high-risk patients not an agressive course of treatment should be used. These models can be used as a classifier by compar- A reviewer suggested the use of Harrell’s c-index as ing the predicted survival times to a chosen threshold. an alternative measure to RMSE. This measure consid- Since the clinically high-risk group is notorious for hav- ers the relative ordering of predicted survival times with ing poor prognosis, our goal is focused on subclassifying the observed times [21]. We find that models provide pre- these patients. A threshold of 2 years is used. If a patient dicted times that are strongly concordant with observed has a predicted survival time less than 2 years, they are times (Figs. 1 and 2), which indicates an accurate rela- labeled as LPS (low predicted survival). Otherwise, they tive ordering of patients. These results suggests that the are non-LPS. A classifier is considered successful if the models may be useful as a classifier. two resulting groups (LPS versus non-LPS) have distinct Grimes et al. Biology Direct (2018) 13:11 Page 7 of 15 survival curves. The Kaplan-Meier estimates [22]ofthese between LPS and non-LPS groups than is found with OS curves for each RNA-seq model are shown in Figs. 3, 4, (Figs. 3, 4, 5 and 6). The T-1 model provides the largest 5 and 6, and the null model and ensemble are shown in distinction in 2-year EFS rates: 0.29 ± 0.06 versus 0.56 ± Fig. 7. 0.10 (Table 1). Using OS as the outcome, almost every RNA-seq model In general, subclassification is more successful with OS is able to partition high-risk patients into two distinct than with EFS. The ensemble approach (Fig. 7)reflectsthe groups, providing a substantial improvement over the null overall performance in both cases: the LPS and non-LPS model. The TI-4 model produces groups with the largest groups are well separated by the ensemble in OS (0.42 ± difference in 2-year OS rates: 0.40±0.11 versus 0.80±0.05 0.10 versus 0.82 ±0.05) but not for EFS (0.36 ±0.06 versus (Table 1). With EFS as the outcome, there is less separation 0.39 ± 0.09) (Table 1). Fig. 3 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the gene annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test Grimes et al. Biology Direct (2018) 13:11 Page 8 of 15 Fig. 4 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the transcripts annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test Pathway analysis fit using elnet are considered, as these contain an amount Pathway enrichment analysis provides a biological sum- of sparsity that is appropriate for pathway analysis. Two mary of the genes selected by the AFT model. Gene sets gene sets are constructed, one associated with OS and the are constructed by collecting the predictors with nonzero other with EFS. Pathway enrichment analysis (on KEGG coefficients in the fitted G-4, T-4 and TI-4 models. The I- pathways) is performed using DAVID 6.8 [23] and sum- 4 model with introns only is not considered, since introns marized in Tables 2 and 3. cannot easily be interpreted in the pathway analysis. The When predicting OS, a total of 354 unique genes are PLS and SPLS methods gave each predictor some weight given nonzero coefficients by one of the three models. in the AFT model, while the predictors selected by lasso Of these genes, 186 are annotated in KEGG pathways. are a subset of those selected by elnet. Hence, only models DAVID uses a modified fisher exact test to compute Grimes et al. Biology Direct (2018) 13:11 Page 9 of 15 Fig. 5 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the introns annotation from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test p-values for enrichment, and the Benjamini-Hochberg genes annoted in KEGG pathways. However, the RNA-seq correction is applied to account for multiple testing [24]. data used in this study are filtered based on the GeneCards Two pathways are found to be significantly enriched: Path- database. Hence, the pathway enrichment may be more ways in Cancer and ErbB signaling pathway (Table 2). appropriately conducted using those GeneCard genes For EFS, 246 unique genes have nonzero coefficients, of as the background. The GeneCards database contained which 135 are annoted in KEGG pathways. However, no 3512 genes related to neuroblastoma, of which 2044 pathways are enriched for EFS at the 0.05 significance are annoted in KEGG pathways. Relative to this back- level. ground, three pathways are enriched for OS: ErbB sig- The preceeding enrichment analysis uses the entire naling pathway, Salivary secretion, and Inflammatory human genome as a background, which contains 6910 mediator regulation of TRP channels (Table 3). Five Grimes et al. Biology Direct (2018) 13:11 Page 10 of 15 Fig. 6 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using both the transcript and intron annotations from the RNA-seq data. Rows 1-4 correspond to PLS, SPLS, lasso, and elnet fitting procedures. The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test pathways are enriched for EFS: Terpenoid backbone procedure, to predict overall survival (OS) and event-free biosynthesis; Metabolic pathways; Valine, leucine and survival (EFS) times of neuroblastoma patients. Three fea- isoleucine degradation; Biosynthesis of antibiotics; and ture levels of an RNA-seq dataset were considered, includ- Fatty acid metabolism (Table 3). These pathways have p- ing genes, transcripts, and introns. Models were fit using valuesbelowthe0.05significancelevel,butarenonsignif- the three features independently and with transcripts and icant after applying the Benjamini-Hochberg correction. introns together. In terms of RMSE, the predictive performance of OS Discussion is greatly improved in the RNA-seq models over the null In this study we used the AFT model, fit using various model, but this improvement is curtailed when predict- dimension reduction techniques and a dataset imputation ing EFS. The high rate of censoring that is found in this Grimes et al. Biology Direct (2018) 13:11 Page 11 of 15 Fig. 7 Kaplan-Meier estimates for HR and LPS2. Kaplan-Meier estimates for overall survival (left column) and event-free survival (right column) of clinically high risk patients using the null model (first row) and the ensemble approach (second row). The orange line corresponds to patients labeled as LPS2 (predicted survival time less than 2 years), and blue lines are non-LPS2. The p-values are for the logrank test data will be a hinderance for any nonparametric model. OS, two cancer-related pathways are enriched. This anal- Alternative approaches can be considered: One possibil- ysis may be biased, however, since the RNA-seq data are ity is to switch to semiparametric estimation, but this initially filtered using the GeneCards database. If the back- ground is altered to reflect this filtering, we find that approach will be computationally intensive in this high- dimensional setting. A more practical solution may be to one of the two cancer-related pathways remains relatively employ a boosting algorithm (see [25] for example). These enriched. This alteration also reveals additional enriched alternatives were not explored in detail in this paper. pathways for the OS and EFS gene sets, but their relevance The second goal is to subclassify clinically high-risk to neuroblastoma is questionable. Since the prediction (HR) patients. In this venture the AFT model produces of EFS had limited success, it is no surprise that the very promising results. High-risk patients with low sur- genes selected for EFS appear to have limited biological vival times are more sensitive to the amount of error relevance. remaining in predicted times, but the estimates tend to The predictive accuracy and pathway enrichment for OS be in the right direction. That is, the relative ordering suggests that the AFT model with elastic net is able to pick of the patients by their predicted survival times is accu- out biologically meaningful genes. A future study pursu- rate. A reviewer suggested the use of Harrell’s c-index ing this kind of interpretation will need to consider the [21] to measure this effect. The c-index is above 0.8 for stochastic nature of the fitting procedure and determine a each model when predicting OS, indicating strong con- stable set of genes selected by the model. As suggested by a cordance between predicted OS time and true OS times reviewer, we can also explore relationships between these (Fig. 1). The concordance is less strong when predicting genes and those excluded by the initial filtering process. EFS (Fig. 2). Such an investigation may produce biological insights into Using a cutoff of 2 years, each model is converted to a the subgroups of high-risk patients. classifier. The TI-4 model provides the best results for OS. An ensemble of models was considered, which incorpo- For EFS, the I-4 model appears to be the best. A classi- rates bagging with rank aggregation of three performance fier using 5 years as a cutoff is also considered, but the measures. The performance of the ensemble method is performance is not as good; setting the threshold to a comparable to that of the best individual model. This suggests that the ensemble method is able to effectively value below 5 years seems to be necessary in order to iden- tify those patients who are at the highest risk in the HR combine models fit on separate datasets. If additional group. datasets are incorporated, such as copy number varia- A pathway analysis of the gene sets selected by the elas- tion or other -omics data, the AFT model can be fit tic net when predicting OS and EFS is performed. With by simply concatenating the datasets together, but the Grimes et al. Biology Direct (2018) 13:11 Page 12 of 15 Table 1 Summary of Kaplan-Meier estimates for 2-year OS and 2-year EFS for clinically high-risk patients using each of the 18 proposed models LPS non-LPS Outcome Prob SE N Prob SE N Data Model P-value OS 0.57 0.09 31 0.79 0.05 59 Null model 0.574 OS 0.60 0.07 57 0.91 0.05 33 G 1 0.081 OS 0.60 0.06 60 0.93 0.05 30 G 2 0.082 OS 0.55 0.09 29 0.79 0.05 61 G 3 0.166 OS 0.50 0.10 24 0.79 0.05 66 G 4 0.164 OS 0.57 0.07 55 0.94 0.04 35 T 1 0.128 OS 0.56 0.07 47 0.88 0.05 43 T 2 0.002 OS 0.41 0.12 17 0.78 0.05 73 T 3 0.048 OS 0.57 0.09 33 0.79 0.05 57 T 4 0.081 OS 0.61 0.07 56 0.88 0.06 34 I 1 0.124 OS 0.51 0.08 36 0.85 0.05 54 I 2 0.027 OS 0.47 0.10 25 0.81 0.05 65 I 3 0.137 OS 0.43 0.10 25 0.82 0.05 65 I 4 0.017 OS 0.57 0.07 56 0.94 0.04 34 TI 1 0.092 OS 0.54 0.07 49 0.90 0.05 41 TI 2 0.001 OS 0.40 0.11 21 0.81 0.05 69 TI 3 0.036 OS 0.40 0.11 22 0.80 0.05 68 TI 4 0.009 OS 0.42 0.10 25 0.82 0.05 65 Ensemble 0.016 EFS 0.39 0.08 37 0.36 0.07 53 Null model 0.616 EFS 0.35 0.07 53 0.39 0.08 37 G 1 0.299 EFS 0.34 0.07 52 0.41 0.08 38 G 2 0.266 EFS 0.30 0.06 58 0.50 0.09 32 G 3 0.125 EFS 0.31 0.06 59 0.48 0.09 31 G 4 0.125 EFS 0.29 0.06 62 0.56 0.10 28 T 1 0.036 EFS 0.35 0.07 51 0.40 0.08 39 T 2 0.946 EFS 0.29 0.07 38 0.43 0.07 52 T 3 0.065 EFS 0.31 0.07 47 0.43 0.08 43 T 4 0.150 EFS 0.35 0.06 66 0.43 0.10 24 I 1 0.256 EFS 0.35 0.06 66 0.43 0.10 24 I 2 0.256 EFS 0.29 0.07 43 0.45 0.07 47 I 3 0.052 EFS 0.29 0.07 42 0.44 0.07 48 I 4 0.082 EFS 0.31 0.06 57 0.47 0.09 33 TI 1 0.084 EFS 0.33 0.06 68 0.50 0.11 22 TI 2 0.183 EFS 0.29 0.07 42 0.45 0.07 48 TI 3 0.062 EFS 0.30 0.07 44 0.44 0.08 46 TI 4 0.085 EFS 0.36 0.06 58 0.39 0.09 32 Ensemble 0.599 Patients with predicted survival of less than 2 years are labeled as Low Predicted Survival (LPS), and otherwise are non-LPS. Columns labeled “Prob.”, “SE”, and “N” correspond to the estimated probability of 2-year survival, the standard error of the estimate, and the number of patients in the given cohort. The P-values are for the logrank test comparing LPS to non-LPS survival. The “Data” column refers to the type of RNA-seq data used, and the “Model” column refers to the dimension reduction technique used computational requirement quickly becomes too bur- that this heuristic works well in combining different anno- densome. The ensemble approach may provide a useful tations of RNA-seq data, but further investigation is heuristic for combining several datasets. We have shown needed to verify the performance with disparate datasets. Grimes et al. Biology Direct (2018) 13:11 Page 13 of 15 Table 2 Pathway enrichment analysis of genes selected by the is a very well-written paper. Overall, the analysis is scien- G-4, T-4, and TI-4 models when predicting OS (no pathways were tifically compelling and relies on creative applications of significantly enriched for EFS) sound statistical techniques. The classifier comparing the Outcome Pathway Count Size P-value BH predicted survival times to a 2-year threshold is success- OS Pathways in cancer 26 393 < 0.001 0.010 ful when it is based on transcript and intron annotations. The ensemble method and its potential application to fit- OS ErbB signaling pathway 11 87 < 0.001 0.012 ting disparate datasets holds much promise for future The columns include KEGG pathway name, the number of genes in the gene set work. that are in the pathway, the total number of genes annotated for the pathway, the p-value from a modified fisher’s exact test, and the Benjamini-Hochberg corrected Reviewer comment: As a suggestion for future p-value research, but entirely unrelated to the current paper which is more than satisfactory, I have the following suggestion. Conclusion From the second paragraph of the Discussion, it appears In this study, we explored the performance of the AFT that it may be helpful to explore Harrell’s C-index as an model in predicting survival times for neuroblastoma alternative measure of accuracy. This may be a better mea- patients. A classifier was constructed by comparing pre- sure than RMSE for the parametric models, especially dicted survival times to a 2-year threshold. Using both because they appear to get the relative ordering of the transcript and intron annotations in the model gave the survival times right rather than the actual magnitudes. best performance. We are able to subclassify clinically Author’s response: We thank Dr. Guha for this sug- high-risk patients into two distinct groups, one with a 40% gestion. The performance of each model using Harrell’s 2-year overall survival rate and the other at 80%. This c-index has been added to the revised manuscript. suggests that the AFT model is useful in subclassifying Reviewer comment: On Line 7 of page 2, should the high-risk patients, which can help clinicians in choosing comma following INSS be deleted? 2. On Line 7 of page 6, effective treatment plans. Only RNA-seq data was con- what is K? sidered in this study, but other types of data can be used Author’s response: Grammatical corrections have been as well. The ensemble method is a useful heuristic for made to the manuscript. For the latter point, there are K = combining several high-dimensional datasets under this 3 performance measures in this study. This is now clarified framework, and it has been shown capable of maintaining in the text. optimal performance. Reviewers’ comments Reviewer’s report 2: Isabel Nepomuceno, Universidad de Reviewer’s report 1: Subharup Guha, University of Florida, Sevilla, Seville, Spain Gainesville, USA In this paper, authors used the accelerated failure time (AFT) model with four dimension reduction techniques The authors explore the performance of the AFT model in and a dataset imputation scheme to predict overall sur- predicting survival times for neuroblastoma patients. This vival and event-free survival times of neuroblastoma Table 3 Pathway enrichment analysis of genes selected by the patients. Three feature levels of and RNA-Seq dataset G-4, T-4, and TI-4 models were considered. Authors shown that the use of RNA-Seq Outcome Pathway Count Size P-value BH data improves accuracy in comparison to using clinical data alone. In general the paper is appropriate to the OS ErbB signaling pathway 11 60 0.029 0.999 journal. The analysis presented in this paper is very inter- OS Salivary secretion 6 23 0.042 0.995 esting. I have several suggestions and comments to be OS Inflammatory mediator 9 48 0.049 0.983 revised: regulation of TRP channels Reviewer comment: The Method section is written in a clear manner but is difficult to reproduce. Authors EFS Terpenoid backbone 4 8 0.010 0.906 biosynthesis mentioned the R package used but they don’t provide the EFS Metabolic pathways 29 304 0.016 0.847 Rcodeofthe study. Author’s response: We thank Dr. Nepomuceno for her EFS Valine, leucine and 5 20 0.032 0.911 isoleucine degradation comments and suggestions. All R code and output is available from GitHub at https://github.com/tgrimes/ EFS Biosynthesis of 12 98 0.037 0.882 antibiotics CAMDA-2017-Neuroblastoma. The session info is also EFS Fatty acid metabolism 5 21 0.037 0.820 reported, which includes the R version, computer specifica- tions, and a list of the packages used during the analysis. In this analysis, the GeneCards genes are used at the background. The columns include survival outcome (OS or EFS), KEGG pathway name, the number of genes in Reviewer comment: The Ensemble Method subsec- the gene set that are in the pathway, the total number of genes annotated for the tion, authors use bagging with rank aggregation over each pathway, the p-value from a modified fisher’s exact test, and the performance measure and set B to 20. Why this parameter Benjamini-Hochberg corrected p-value Grimes et al. Biology Direct (2018) 13:11 Page 14 of 15 is fixed to 20 should be explained. And authors should the GeneCards database to subset on genes, which would explain why the use bagging instead of cross validation. bias our selection to genes in cancer-related pathways. In Author’s response: The choice of 20 iterations for bag- response, we have modified this section and now conduct ging is a compromise between computation time and model a pathway enrichment analysis. However, a question is performance. We also considered B = 50 but did not find a raised regarding the choice of background: should our gene sets be compared to all genes in the genome (as is usu- substantial change in performance. Reviewer comment: The description of the RNA-Seq ally done) or to the GeneCards genes that we subset on? Data, authors reduce the "raw data" with 60776 genes With the former, there is a concern that the analysis may be into 3401 using the 3681 genes related to neuroblastoma biased. Results for both of these scenarios have been added obtained from the Gene Cards Suite. Have authors made to the manuscript. some analysis from the remaining genes? Could be genes Reviewer comment: Finally, as minor comments: - The related with the problem and not related with the disease? Bibliography Section must be revised, there are some It could be interesting to do a cluster analysis to see if the incomplete reference as for example number 14. - In grouped genes using prior knowledge are also clustered Table 1,one ofthemodelsisnamedsimpleforthe baseline together in this analysis. model. It should be names null model as authors explained Author’s response: These are interesting suggestions that before. deserve a separate analysis to be fully addressed. The Author’s response: The bibliography section has been main purpose in using the Gene Cards database was to corrected, and the tables and figures have been relabeled provide an initial filtering to speed up computation. We to be consistent with the text. also re-ran the analysis without this step and found lit- tle difference in predictive performance. We are careful Abbreviations not to place too much emphasis on the interpretation of AFT: Accelerated failure time; CI: Confidence interval; EFS: Event-free survival; elnet: Elastic net; HR: High-risk; INSS: International neuroblastoma staging the gene sets obtained in this analysis. As you’ve pointed system; lasso: Least absolute shrinkage and selection operator; LPS: Low out, there are many new questions that have been uncov- predicted survival; OS: Overall survival; PLS: Partial least squares; RMSE: Root ered and deserve careful consideration. We’ve added some mean squared error; SPLS: Sparse partial least squares comments regarding this in the discussion section of the Acknowledgements manuscript. We acknowledge the CAMDA 2017 committee reviewers for their valuable Reviewer comment: Furthermore, a reference about comments. the Cox proportional hazards model or the R package used Funding should be added. Publication of this article was funded in part by the University of Florida Open Author’s response: We thank the author for pointing Access Publishing Fund. out this omission. The revised manuscript now contains Availability of data and materials additional references. The datasets used in this article be accessed from the GEO repository with Reviewer comment: Section Results, classification of series accession number GSE49711. All R code is available on GitHub at https://github.com/tgrimes/CAMDA-2017-Neuroblastoma. high-risk patients should be rewritten. The second and third paragraph is confused and difficult to see which plot Authors’ contributions corresponds with each sentence. All authors took part in research discussion and the final manuscript preparation. TG performed the data analysis and wrote the first draft of the Author’s response: This section has been reworded to paper. SoD planned the study and the selection of statistical methods. All clarify which table or figure each sentence is referring to. authors read and approved the final manuscript. The titles for each plot have been changed in concordance Ethics approval and consent to participate to the labels used to identify each model within the Not applicable. manuscript. Reviewer comment: In section Pathway analysis, Competing interests The authors declare that they have no competing interests. authors claim that several genes are involved in several pathways.Thatmeans,dogenes appear inthepathways Publisher’s Note or are the pathways enriched by the set of genes? If it is Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. the second case, authors should add a table with the list of pathways, the number of entities in the pathways and Received: 12 October 2017 Accepted: 1 May 2018 the number of genes from the set which appear in the pathway. Author’s response: We thank the reviewer for prompt- References 1. Bosse KR, Maris JM. Advances in the translational genomics of ing this clarification. Previously, the interpretation was neuroblastoma: From improving risk stratification and revealing novel that genes appear in the pathways. But this initial biology to identifying actionable genomic alterations. Cancer. approach seems uninformative, particularly since we use 2016;122(1):20–33. Grimes et al. Biology Direct (2018) 13:11 Page 15 of 15 2. Brodeur GM, Seeger RC, Schwab M, Varmus HE, Bishop JM. Amplification of n-myc in untreated human neuroblastomas correlates with advanced disease stage. Science. 1984;224(4653):1121–4. 3. Formicola D, Petrosino G, Lasorsa VA, Pignataro P, Cimmino F, Vetrella S, Longo L, Tonini GP, Oberthuer A, Iolascon A, et al. An 18 gene expression-based score classifier predicts the clinical outcome in stage 4 neuroblastoma. J Transl Med. 2016;14(1):142. 4. Tan Q, Thomassen M, Jochumsen KM, Mogensen O, Christensen K, Kruse TA. Gene selection for predicting survival outcomes of cancer patients in microarray studies. Adv Comput Inf Sci Eng. 2008;1(1):405–9. 5. Boulesteix A-L, Strimmer K. Partial least squares: a versatile tool for the analysis of high-dimensional genomic data. Brief Bioinform. 2006;8(1): 32–44. 6. Chung D, Chun H, Keles S. Spls: Sparse Partial Least Squares (SPLS) Regression and Classification. 2018. R package version 2.2-2. https:// CRAN.R-project.org/package=spls. Accessed 28 Apr 2018. 7. Chun H, Keles¸ S. J R Stat Soc Series B (Stat Methodol). 2010;72(1):3–25. 8. Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Series B (Stat Methodol). 1996;58(1):267–88. 9. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1. 10. Zou H, Hastie T. J R Stat Soc Series B (Stat Methodol). 2005;67(2):301–20. 11. Mostajabi F, Datta S, Datta S. Predicting patient survival from proteomic profile using mass spectrometry data: an empirical study. Commun Stat Simul Comput. 2013;42(3):485–98. 12. Shah J, Datta S, Datta S. A multi-loss super regression learner (msrl) with application to survival prediction using proteomics. Comput Stat. 2014;29(6):1749–67. 13. Datta S. Estimating the mean life time using right censored data. Stat Methodol. 2005;2(1):65–9. 14. Kleinbaum DG, Klein M. Kaplan-meier survival curves and the log-rank test. In: Survival Analysis, 3rd edn. New York: Springer; 2012. p. 55–96. Chap. 2. 15. Therneau TM. A Package for Survival Analysis in S. 2015. version 2.38. https://CRAN.R-project.org/package=survival. Accessed 28 Apr 2018. 16. Su Z, Fang H, Hong H, Shi L, Zhang W, Zhang W, Zhang Y, Dong Z, Lancashire LJ, Bessarabova M, et al. An investigation of biomarkers derived from legacy microarray data for their utility in the rna-seq era. Genome Biol. 2014;15(12):523. 17. Zhang W, Yu Y, Hertwig F, Thierry-Mieg J, et al. Comparison of rna-seq and microarray-based models for clinical endpoint prediction. Genome Biol. 2015;16(1):133. 18. Safran M, Dalah I, Alexander J, Rosen N, Iny Stein T, Shmoish M, Nativ N, Bahir I, Doniger T, Krug H, et al. Genecards version 3: the human gene integrator. Database. 2010;2010:020. 19. Cox DR. J R Stat Soc Series B (Stat Methodol). 1972;34(2):187–220. 20. Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer; 2000. 21. Harrell Jr FE, Califf RM, Pryor DB, Lee KL, Rosati RA, et al. Evaluating the yield of medical tests. J Am Med Assoc. 1982;247(18):2543–6. 22. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53(282):457–81. 23. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using david bioinformatics resources. Nat Protocol. 2009;4(1):44. 24. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B (Methodol). 1995;57(1):289–300. 25. Schmid M, Hothorn T. Flexible boosting of accelerated failure time models. BMC Bioinformatics. 2008;9(1):269.

Journal

Biology DirectSpringer Journals

Published: May 30, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off