A total of 129 patients were retrospectively included in the current study. All patients had clinical information available from the Cancer Genome Atlas GBM Collection (TCGA-GBM¹) and corresponding multimodal MRI data available from the Cancer Imaging Archive (TCIA²). The complete radiological data of the TCIA-GBM consists of 262 multimodal MRI scans obtained from eight institutions (Clark et al., 2013). The data used in the present study included the pre-operative baseline scans acquired with T1WI, T1CE, T2WI, and FLAIR sequences. The patients were assigned to two cohorts: a training cohort comprising 99 patients from institutions 2, 6, 8, and 12, and an independent test cohort comprising 30 patients from institutions 14, 19, 27, and 76. No institutional review board approval was required since TCGA is a publicly available dataset without patient identifiers. Except for the multi-parametric MRI data, the clinicopathologic information, including gender, age, KPS, prognostic treatment, TP53, PTEN, EGFR, and IDH1, were obtained from TCGA GBM Project for all patients.

The pre-processed and labeled MRI data are available through TCIA (Bakas et al., 2017). The pre-operative MRI volumes were initially co-registered to the T1WI and then skull-stripped. Consistent with the BRAin Tumor Segmentation (BRATS) challenge, the segmentation labels delineate three parts of each tumor, including enhanced (ET) and non-enhanced tumor (NET) regions, as well as the peritumoral edema (ED). Tumor segmentation was achieved by integrating an automated pre-segmentation process with manual corrections from a board-certified neuroradiologist. The automated hybrid generative-discriminative method, which won first place at the International Multimodal BRATS 2015, was applied to produce a set of labels. Then, these automatically segmented labels were modified and the misclassified labels were manually corrected.

A variety of imaging features were extracted from the MRI datasets, including spatial information, intensity, volumetric, morphologic, histogram, and textural parameters. A total of 537 radiomics features, including 179 features from ET, NET or ED, were extracted for each subject. A detailed description of the extracted spatial features can be found in Supplementary Table 1.

The number of subjects (independent samples) should exceed the number of selected features by a factor of at least 10 according to the guideline (Harrell, 2015). Thus, the combination of the least absolute shrinkage and selection operator (LASSO) algorithm and Cox survival model was used to select the most useful prognostic features in the training dataset. The LASSO method utilizes a regularization parameter λ to shrink the coefficients of all irrelevant features to zero. Here, λ was optimized to maximize the area under the receiver operating curve (AUROC) in a 10-fold cross-validation procedure. Then, the Cox regression model was constructed from the features with non-zero coefficients. After that, radiomics scores (Radscores) were computed for each patient through a linear combination of the selected LASSO features, weighted by their respective coefficients. These scores can be useful by further stratifying the GBM patients.

Specifically, when considering the influence of different tumor regions, two kinds of feature selection strategies were applied and compared in this study:

By using Kaplan-Meier survival analysis, the potential correlation between radiomics signatures and OS was assessed in the training cohort and validated in the test cohort. The optimal cutoff values of Radscores were selected based on their association with the patients’ OS in the training cohort using the maximally selected rank statistics from the “maxstat” R package (Delgado et al., 2014). Then, the same cutoff values were applied to the test cohorts. Patients were classified into high- and low-Radscore groups accordingly for further analysis. Evaluation of the radiomics signature built from different features was performed by adding the Radscore as an independent factor in the multivariable Cox proportional hazard model (backward step-down selection; the Akaike information criterion, AIC). This model also integrated the general clinical risk factors as well.

Radiomics and clinical nomograms were both applied to the training cohort based on the multivariate Cox analysis to indicate the increased predictive value of the radiomics signatures to the clinical risk factors for individualize OS assessment. Specifically, the radiomics nomograms were generated by integrating the radiomics signatures and clinical risk factors, while the clinical nomogram contained only the clinical risk factors. Then, the increased predictive value of the radiomics signatures relative to clinical risk factors was assessed by calculating the discrimination, reclassification, calibration, and clinical usefulness. The effect of the radiomics signatures was compared with that of the clinical nomogram. Before this, two radiomics nomograms corresponding to two feature selection strategies were also constructed and compared.

For model fitting assessment, AIC was calculated to assess the risk of model overfitting. And the performance of different models was evaluated using the integrated Brier score (IBS) by calculating the prediction error over time and the compared, which represents the weighted average squared distance between the predicted probability of the established model and the observed survival status. IBS values can range from 0, i.e., a perfect model, to 0.25, i.e., a non-informative model with a 50% survival/non-survival prediction.

For discrimination performance assessment, calibration curves were derived to compare the consistency between the OS predicted by the radiomics and clinical nomograms and the actual OS. The Harrell concordance index (C-index) was calculated to quantify the discrimination performance. Meanwhile, decision curve analysis (DCA) enabled the determination of clinical usefulness of the radiomics nomograms by measuring the net benefit at several probability thresholds.

Statistical analysis performed using the R software. R packages that were utilized are summarized as follows: the “surv_cutpoint” function of “survminer” package was applied for cutoff point calculation; the “glmnet” and “survival” packages were applied for LASSO Cox regression; the “rms” package was applied for multivariate Cox regression, calibration plots and nomograms; the “Hmisc” package was applied for C-indexes comparisons; the “rmda” package was applied for DCA; and the “ipred” package was applied for IBS calculation.

Clinical characteristics of all patients are shown in Supplementary Table 2. We used the LASSO Cox regression model to select features and build prognostic classifiers from these features. Ten potential predictors were identified from the 179 features for each subregion, i.e., ET, NET, and ED, in the training cohort Eight potential predictors were identified when all the features were concatenated together. The radiomics signatures were then constructed with weights applied to the coefficients in each model and the Radscores were obtained. The relevant equations for this process can be found in Supplementary Table 3.

The distributions of the Radscores of different sub-regions in the training cohort are shown in Figure 1. The optimal cutoff value was −0.233, −0.295, and 0.055 for the Radscore-ET, Radscore-NET and Radscore-ED, respectively, and −0.321 for the Radscore-Con. Accordingly, patients were stratified into a low-Radscore group and high-Radscore group. The OS rate in the high-Radscore group and low-Radscore group in the training and test cohorts are provided in Supplementary Table 4. Combined with Figure 2, we found that patients with lower Radscores generally exhibited a better OS. This finding was assessed in the training cohort and then confirmed in the test cohort.

To validate the radiomics signatures and to provide a quantitative and clinical method to predict the probability of the 1-, 3-, and 5-year OS of GBM patients, a clinical nomogram and two radiomics nomograms were constructed based on the training cohort data (Figure 3). Four clinical risk factors, including days to birth, KPS, prognostic treatment and spatial_frontal, were significantly associated with survival and were included in the models (Table 1).

As for the radiomics nomograms, we found that the discrimination performance of the one constructed from all the independent RadScore risk factors, including Radscore-ET, Radscore-NET, and Radscore-ED (C-index/AIC: training cohort: 0.717/452; test cohort: 0.655/143) was better than the nomogram constructed from the Radscore-Con as a single risk factor (C-index/AIC: training cohort: 0.656/468; test cohort: 0.535/146). Thus, the nomogram constructed from the combination of the Radscore-ET, Radscore-NET, and Radscore-ED signatures was used for further analysis.

Meanwhile, the correlation between the subregion RadScores was examined to test whether overfitting exist. The correlation coefficients between RadScore-ET and Radscore-NET, between Radscore-ET and Radscore-ED, and between Radscore-NET and Radscore-ED were 0.36, 0.28, and 0.16, respectively (Figure 4A). Figure 4B indicates the scatter plot based on these three RadScores and 1-year survival difference between two groups can be observed obviously.

We found good agreement between three actual survival observations and the survival estimates determined by the radiomics nomograms at 1, 3, and 5 years in both training and test cohorts, as is depicted in the calibration curves of the nomograms (Figure 5). The C-index, IBS and AIC estimations for the clinical and radiomics nomogram models are also summarized in Table 2. As assessed with the IBS (training cohort/test cohort: 0.077/0.066, lower values indication better model performance), the C-index (training cohort/test cohort: 0.717/0.655, higher values indication better discriminative ability) and AIC (training cohort/test cohort: 452/143, lower values indicating better model performance), the radiomics nomogram performed better than the clinical nomogram in both the training and test cohorts.

Meanwhile, multivariable Cox regression analysis after adjustment for clinical factors showed that the Radscores were an independent and powerful prognostic factor for OS in the training and test cohorts (Supplementary Figure 1). The p-values for the Radscore-ET, Radscore-NET, and Radscore-ED were 0.96, 0.002, and 0.003, respectively. Among these, the Radscore-ED had the highest hazard ratio and played an important role in prognosis.

The DCA (Figure 6) showed that in most of the reasonable threshold probabilities, the net benefit of the radiomics nomogram was slightly higher than that of the clinical nomogram.