Preoperative MRI-radiomics features improve prediction of survival in glioblastoma patients over MGMT methylation status alone.

BACKGROUND: Glioblastoma (GBM) is the most common malignant central nervous system tumor, and MGMT promoter hypermethylation in this tumor has been shown to be associated with better prognosis. We evaluated the capacity of radiomics features to add complementary information to MGMT status, to improve the ability to predict prognosis.
METHODS: 159 patients with untreated GBM were included in this study and divided into training and independent test sets. 286 radiomics features were extracted from the magnetic resonance images acquired prior to any treatments. A least absolute shrinkage selection operator (LASSO) selection followed by Kaplan-Meier analysis was used to determine the prognostic value of radiomics features to predict overall survival (OS). The combination of MGMT status with radiomics was also investigated and all results were validated on the independent test set.
RESULTS: LASSO analysis identified 8 out of the 286 radiomic features to be relevant which were then used for determining association to OS. One feature (edge descriptor) remained significant on the external validation cohort after multiple testing (p=0.04) and the combination with MGMT identified a group of patients with the best prognosis with a survival probability of 0.61 after 43 months (p=0.0005).
CONCLUSION: Our results suggest that combining radiomics with MGMT is more accurate in stratifying patients into groups of different survival risks when compared to with using these predictors in isolation. We identified two subgroups within patients who have methylated MGMT: one with a similar survival to unmethylated MGMT patients and the other with a significantly longer OS.

INTRODUCTION

Glioblastoma (GBM) is the most common malignant central nervous system tumor (46.6%) and has a poor prognosis with only 5.5% (95% CI 5.2%-5.8%) survival at 5 years [1]. Imaging plays an important role in determining patient diagnosis, prognosis, risk assessment, response evaluation and follow up [2]. The standard treatment consists of surgical resection followed by radiotherapy with concomitant and adjuvant chemotherapy (temozolomide) [3]. Recently, several molecular biomarkers have been identified to predict prognosis and outcomes to therapy as well as to guide personalized therapies [4, 5]. The hypermethylation of the MGMT promoter (O6-methylguanine–DNA methyltransferase), which results in the silencing of this DNA repair enzyme, provides key prognostic data and has been shown to be associated with better prognosis and better response to temozolomide [6, 7]. We studied the utility of combining pre-treatment, diagnostic magnetic resonance imaging (MRI)-based radiomics features with MGMT methylation to further stratify patient prognosis.

Previous studies have applied qualitative visual feature descriptors assigned by radiologists from MRI scans, such as those based on the VASARI (Visually AcesSAble Rembrandt Images) criteria, showing that these may predict patient outcomes [8, 9]. Recently, radiomics features of MR images have been used to extract quantitative measures of whole tumor appearance heterogeneity to predict outcomes in multiple cancers including GBM [10]. Radiomics involves the extraction of quantitative feature descriptors of tumor appearance through edges, textures, shape and histogram features through automated high-throughput computer analysis of medical images and has been applied to predict outcomes in several cancers [10-12]. In brain cancer, radiomic features have shown feasibility to predict outcomes including survival [13-18].

Radiomics when applied to predict underlying genotypic features is called radiogenomics [19-22] and has been used to non-invasively detect association of images with underlying molecular markers such as MGMT methylation, IDH1 mutation and 1p/19q co-deletion [13, 16, 18, 23, 24]. In our study, we evaluated whether the combination of the molecular biomarker MGMT and image-based tumor heterogeneity computed using radiomics analysis would produce more accurate overall survival (OS) classification than MGMT alone. We focused on the capacity of radiomics to improve the effectiveness of MGMT methylation as a biomarker to identify patient groups with distinct survivals among the patients with or without MGMT methylation. Our aim was to build an easy-to-compute and interpretable model incorporating the most relevant radiomics features to augment the capacity of MGMT methylation as biomarker of outcome.

Analysis consisted of several steps including relevant features selection followed by a univariate analysis to construct models that could potentially be incorporated into routine clinical workflows to inform patient prognosis and guide treatment decision making. All the analyses were performed on a training cohort and then validated on an independent validation cohort following the protocol described in Figure 1.

Figure 1

Workflow of the study divided on four steps

Training cohort: 1) Radiomics feature selection using the Least Absolute Shrinkage and Selection Operator (LASSO) regression model and prognosis value of these features for the OS prediction 2) Prognosis value of MGMT and combination of MGMT and radiomics features for the OS prediction. Validation cohort: 3) Validation of the prognosis value of the selected radiomics features for the OS prediction. 4) Prognosis value of MGMT and the combination of MGMT and radiomics features for the OS prediction on the validation cohort.

RESULTS

Patient characteristics

Median patient age at diagnosis was 61.5 years (range, 18-87) and 40% were female. Seventy percent of the patients had a KPS score >80%. Among the 121 patients with known MGMT status, 40% were methylated. A large majority of the patients were IDH wild-type (n=137). All patient characteristics including details specific to training and external validation cohorts are summarized in Table 1.

Table 1

Patient characteristics and outcome

		Training Cohort (n=98)	Validation Cohort (n=61)	All (n=159)
Age (median [range])		64y [20–87]	61y [18–80]	61.5y [18–87]
Gender
	F	38 (39%)	26 (43%)	64 (40%)
	M	60 (61%)	35 (57%)	95 (60%)
KPS (%)
	50%	1 (1%)	1 (2%)	2 (1%)
	60%	5 (5%)	7 (11%)	12 (7%)
	70%	23 (23%)	2 (3%)	25 (16%)
	80%	26 (27%)	25 (41%)	51 (32%)
	90%	31 (32%)	4 (7%)	35 (22%)
	100%	12 (12%)	13 (21%)	25 (16%)
	Unknown	0 (0%)	9 (15%)	9 (6%)
MGMT status
	Unmethylated	57 (58%)	15 (25%)	73 (46%)
	Methylated	29 (30%)	20 (33%)	49 (31%)
	Unknown	12 (12%)	26 (42%)	37 (23%)
IDH status
	Wild-type	90 (92%)	47 (77%)	137 (86%)
	Mutant	5 (5%)	3 (5%)	8 (5%)
	Unknown	3 (3%)	11 (18%)	14 (9%)
Extent of resection
	Biopsy	0 (0%)	10 (16%)	10 (6%)
	Subtotal resection	52 (53%)
	Near total resection	16 (16%)	51 (84%)	149 (94%)
	Gross total resection	30 (31%)
OS
	median [range]	587d [79 - 1324]	478d [6 - 2368]	567d [6 - 2368]
	1yr survival probability	71%	63%	68%
	2yr survival probability	43%	25%	36%
	3yr survival probability	28%	17%	26%
Vital Status
	Dead	56 (57%)	45 (74%)	101 (64%)

Feature selection

286 MRI-radiomics features (Figure 2D) were used in a least absolute shrinkage selection operator (LASSO) regression model to reduce the number of features in the training cohort to identify the features most relevant to outcome. Different bin sizes of MRI intensity were found to have a limited impact on feature selection: the model using a 128-bin setting selected the same set of features as the model using a 32 or 64-bin setting, and there was also a large overlap in the set of selected features with that of the model using a 16-bin setting (Supplementary Table 1).

Figure 2

From the image acquisition to the radiomics features extraction

(A). Acquisition of T1WI, T1WICE and FLAIR sequences. (B). Image segmentation using simultaneously the information from the 3 sequences (in red) and reporting this segmentation on all sequences. (C). Image post treatments: Gabor filtering and binning. (D). Radiomics feature extraction.

The model based on a 128-bin setting selected 8 MRI-radiomic features as relevant to outcome: 1 feature from fluid attenuation inversion recovery (FLAIR) images (intensity histogram feature category: interquartile range), 2 features from T1-weighted images (T1WI) images (Gray-Level Co-occurrence Matrix (GLCM) feature category: cluster prominence and difference variance), 2 features from the contrast-enhanced T1-weighted (T1WICE)(Gabor feature category: Gθ=0, f=2 Skewness and Gθ=0, f=2 Mean) and 1 feature from the T1WICE images (intensity histogram feature category: 10th percentile), and 2 shape features (sphericity and surface-to-volume ratio).

Overall survival analysis in training cohort

MGMT methylation status

MGMT methylation alone was associated with improved OS with a median survival of 36.0 months for methylated patients and 16.8 months for unmethylated patients (p=0.02) (see Figure 3A).

Figure 3

Overall survival curves for (A) Methylation of the MGMT promoter (n=86 patients); (B) Edge enhancement the tumor in T1WICE images (skewness of histogram after Gabor filtering with a direction of 0° and a frequency of 2 – patients are split according to the optimal cut-off value τ = −0. 49), (n=98 patients); (C) Combination of the methylation status of the MGMT promoter and edge enhancement in T1WICE images (skewness of histogram after Gabor filtering with a direction of 0° and a frequency of 2, the reported p-value is for the difference in prognosis for the patients with methylation and heterogeneity under the median vs. the rest, (n=86 patients).

Radiomic features

Table 2 presents the optimal thresholds for each selected feature obtained with Receiver Operating Characteristic (ROC) curve analysis. Of the 8 selected features, 4 were significantly associated with OS after adjusting for multiple comparisons using the Bonferroni procedure: 2 shape features (sphericity and surface-to-volume ratio) and 2 features from T1WICE images (Gθ=0, f=2 Skewness and Gθ=0, f=2 Mean). Tumors with high sphericity (> 0.73) or surface-to-volume ratio < 1.73 were associated with a longer median survival (36.0 months) vs. tumors with low sphericity (17.5 months) and tumors with large surface-to-volume ratio (15.5 months), respectively. Patients with the highest edge enhancement on T1WICE (i.e., Gθ=0, f=2 Skewness value lower than the optimal threshold [−0.49]) had a longer median survival (36.0 months) vs. their counterpart (16.3 months) (Figure 3B). Median survival times are reported in Table 2.

Table 2

Overall survival analysis in the training and validation sets. The 8 features selected by the LASSO regression model were used in the training set and the 4 significant features in the training set were used in the validation set. All p-values were corrected for multi-testing by the Bonferroni procedure

				Testing set				Validation set
				median survivals (95% LCL; 65% UCL) (months)				median survivals (95% LCL; 65% UCL) (months)
Image sequence	Feature category	Feature name	cut-off value	group 1 (≥ cut-off)	group 2 (<cut-off)	p-values†	group sizes	group 1 (≥ cut-off)	group 2 (< cut-off)	p-values†	group sizes
FLAIR	Intensity histogram	Interquartile range	175.5	21.7 (18.4; NA‡)	16.3 (12.1; 29.0)	0.8	52/46	-	-	-	-
T1WI	GLCM	Cluster prominence	1.77E+07	25.9 (21.0; NA‡)	17.5 (14.6; 28.3)	0.32	35/63	-	-	-	-
		Difference variance	253.7	NA‡ (11.6; NA‡)	18.7 (16.3; 25.9)	0.24	22/76	-	-	-	-
T1WICE	Intensity histogram	10th percentile	29.5	18.2 (15.5; 22.2)	NA‡ (NA‡; NA‡)	0.06	72/26	-	-	-	-
	Gabor features	Gθ=0, f=2 Mean	144.1	29.0 (18.4; NA‡)	16.8 (12.1; 22.2)	0.02	54/44	15.7 (13.3; 22.3)	14.1 (6.8; NA‡)	1	51/10
		Gθ=0, f=2 Skewness	-0.49	16.3 (12.4; 21.7)	36.0 (25.9; NA‡)	0.02	56/42	12.2 (8.5; 19.7)	20.6 (15.7; NA‡)	0.04	31/30
NA	Shape	Sphericity	0.73	36.0 (19.3; NA‡)	17.5 (12.4; 22.1)	0.007	47/51	16.0 (13.9; 22.3)	8.5 (6.8; NA‡)	0.8	50/11
		Surface to volume ratio	1.73	15.5 (11.6; 21.0)	36.0 (19.3; NA‡)	0.002	34/64	20.6 (13.3; NA‡)	14.7 (11.7; 20.8)	1	16/45

†Significance of difference in the OS between the two groups (p-values were corrected for multiple comparisons with Bonferroni procedure)

‡NA values appears when the survivorship function does not reach 0.5

95% LCL and 95% UCL are the 95% lower and upper confidence limits around the median survivals

Combination of MGMT and Radiomics

The combination of MGMT methylation status and the selected Gθ=0, f=2 Skewness radiomics feature identified a subgroup of patients with a better prognosis than their counterparts. Patients with MGMT methylation and high negative skewness of Gabor edge enhancement had a better prognosis (survival probability of 0.61 after 43 months of observation) compared with other patients (survival probability of 0.15 and median survival of 17.5 months, (95% CI: 13.4, 21.0)) (p=0.0005). On the other hand, patients with MGMT methylation and Gθ=0, f=2 Skewness > optimal threshold [−0.49]) did not have a significantly different prognosis than patients without MGMT methylation (p=0.5) (Figure 3C). MGMT methylation status combined with the remaining selected features (sphericity, surface-to-volume ratio and Gθ=0, f=2 Mean) did not identify subgroups of patients with significantly different prognosis.

Overall survival in the validation cohort

In the validation cohort, MGMT methylation alone remained a significant prognostic factor with a median survival of 19.7 months for methylated patients vs. 13.3 months for unmethylated patients, p=0.01 (Figure 4A).

Figure 4

Overall survival curves for a validation cohort (A) Methylation of the MGMT promoter (n=35 patients); (B) Edge enhancement in T1WICE images (skewness of histogram after Gabor filtering with a direction of 0° and a frequency of 2). Patients are split using the optimal cut-off value of −0.49 obtained from the training cohort (n=61 patients); (C) Combination of the methylation status of the MGMT promoter and edge enhancement in T1WICE images (skewness of histogram after Gabor filtering with a direction of 0° and a frequency of 2, the reported p-value is for the difference in prognosis for the patients with methylation and heterogeneity under the median vs. the rest, (n=35 patients).

Using the optimal thresholds identified in the training dataset, only the Gθ=0, f=2 Skewness feature from T1WICE images remained significantly associated with OS after adjusting for multiple comparisons with the Bonferroni procedure (Table 2). Similar to the training set, patients with highest edge enhancement on T1WICE (Gθ=0, f=2 skewness< optimal threshold [−0.49], p=0.04) had a longer median survival of 20.6 months compared with 12.2 months for the group with low edge enhancement (Figure 4B).

Combination of MGMT and Gθ=0, f=2 skewness

The validation set confirmed the results of the training set: those patients who had both MGMT methylation and high negative skewness of Gabor edge enhancement had a longer median OS of 22.7 months vs. patients who did not and who had a median OS of 12.2 months (p=0.004) (Figure 4C).

DISCUSSION

Our study differed from the typical radiogenomic studies which aim to find correlations between radiomics features with genotypic expressions [16, 18, 23, 24]. Instead, we examined the synergistic value of adding radiomics to MGMT methylation status to predict prognosis. Increasingly, studies have employed next generation sequencing to show the predictive utility of underlying tumor genotypic and transcriptomic profiles for predicting survival [25, 26]. Next generation sequencing studies are not readily available at all centers, however, and may be limited by the available amount of tumor tissue. In contrast, MGMT status is widely available via PCR based methods. We therefore tested our hypothesis that textural intra-tumor heterogeneity, tumor shape and intensity histogram features computed from radiomics analysis could be used to augment efficacy of predictive biomarkers such as MGMT status.

Results from our study showed that Gabor features extracted from T1WICE MRI complements the findings using qualitative visual measures using the VASARI criteria [27]. For instance, the length of the lesion's major axis on FLAIR, the proportion of contrast-enhanced on T1WICE and ring enhancement assessed from T1WICE images have shown to be associated with survival [28] and/or MGMT promoter methylation in GBM [29].

We developed a model combining radiomics and MGMT methylation status to predict OS and showed that the combined model predicted survival in an independent validation cohort using the median feature values extracted from the training cohort. Our study identified the skewness of the Gabor filtered edge image (Gθ=0, f=2 Skewness) computed from the T1WICE as a significant predictor of survival. High negative skewness of Gabor edges was associated with longer OS. Gabor edges are directional sensitive edge detectors and large negative skewness indicates that patients with longer OS exhibited a more homogeneous distribution of edges in their tumors compared with those with shorter OS (see Supplementary Figure 1).

In contrast to prior works using radiomics analysis which showed the prognostic value of measures including GLCM [30] and GLSZM (Gray-Level Size-Zone Matrix) features for OS prediction [31], none of these measures were found to be relevant in our study, even if two were selected by the LASSO regression model. We believe that these features were not selected or not found to have a significant prognostic value due to large variations in MR images, such as the voxel size or the magnetic strength, within our cohort and within the multi-institutional cohort. Indeed, it is possible that such effects influenced the value of these features more than patients’ prognosis.

Second, tumor segmentation was performed by simultaneously including all three sequences in order to extract the entire tumor extent including any visible edema on the FLAIR sequence. The Gabor edges were better than the whole tumor summary measures such as the GLCM features to quantify the appearance difference between tumor bulk and edema.

The orientation (0 degree) of the Gabor filter that was selected by our model may be explained with the shape of the tumors that have their longest diameter in the back to front orientation, probably because of the shape of the brain. In addition, the results for this feature were confirmed on the validation cohort using the median value found on the training cohort as the threshold.

The shape features (sphericity and surface-to-volume ratio) were not found to be significantly associated with OS in the validation cohort. This result might be explained by the imbalance in cohort size within the group after splitting the training set by the optimal threshold. Further investigations are needed to evaluate the prognostic value of shape descriptors from MRI images in GBMs.

However, the selection of eight features among the 286 investigated features does not mean that these selected features are the only relevant features. Indeed, LASSO selects and ranks independent features in the order of their relevance for the predicting the outcome variable. In other words, features that are correlated with the top-ranking features are removed, which reduces the false discovery due to multi-testing.

Our analysis identified OS-relevant radiomics features. However, unlike other studies that have reported a significant correlation between MGMT status and radiomics features [23, 32], we found no such correlation. Because MGMT status and OS are indirectly connected, especially for the survival during the first months [7], it is possible that features correlated with MGMT status were not selected in our model.

The combination of the Gθ=0, f=2 Skewness radiomic feature and MGMT status identified a group of patients with a better prognosis in both the training and validation cohorts. These patients had methylated MGMT and a high negative value of Gabor edge skewness. At the opposite, methylated MGMT patients without a high value of Gabor edge skewness had OS indistinguishable from unmethylated MGMT patients. We believe that the radiomic feature quantifies the sharpness of the tumor boundary on the T1WICE and might reflect the degree of tumor infiltration. An interesting extension of this study will be to focus on the correlation between altered genes or pathways and the measurements of tumor edges.

Our radiomics approach differs from others in a few aspects. First, we used automated tumor segmentations with machine learning that was then edited or approved with minimal user input. Second, the tumor segmentations and ultimately the radiomics features were calculated from a single segmented volume combining T1WI, T1WICE and FLAIR. Such a segmented volume in fact includes the edema visible on FLAIR images that will be excluded when using only T1WI, and T1WICE scans. Inclusion of edema components helps the features to quantify the strength of tumor edges.

Our study suffers from some limitations. First, features were extracted without histogram standardization techniques and voxel resampling techniques which could have resulted in the decrease or the loss of the prognostic value of several features [33, 34]; however, this was done because our goal was to build models that could be easily used in a clinical practice across institutions and which would be applicable for new patients without constraints in the image acquisition protocol. Second, we employed one segmentation for all sequences. Extracting features from different segmentations derived from different sequences could have potentially resulted in capturing different aspects of the tumor microenvironment. We plan to pursue such approach especially for correlating image-based features to the genomic pathways in the future. Third, we investigated only one feature selection technique (LASSO) with the dead or alive status of the patients as a binomial outcome. This technique allowed us to identify features that are of potential interest, but it cannot guarantee the best feature selection. However, this method was chosen based on the work of Deist et al., who reported that random forest and elastic net logistic regression yielded higher discriminative performance compared with several other machine learning methods [35]. Fourth, although our cohort consisted mostly of IDH wildtype GBMs, we included a few IDH mutant GBMs that represented 5% of our patients and is concordant with the reported proportion of this rare mutation in primary GBMs [36]. Nevertheless, the model that we developed for untreated GBMs was able to successfully predict survival regardless of IDH status. Supplementary Table 2 provides the repartition of IDH mutant patients into the considered groups and no bias due to these patients was found.

Our study showed the combined value of radiomics with MGMT in predicting OS in patients diagnosed with GBM from pre-treatment MRI using cross-validation training and independent test sets. Our results highlight that radiomics features possibly contain information that is not supplied by molecular markers such as MGMT methylation alone. Particularly, we identified that the methylated MGMT patients can be split into two groups with significantly different survival times. Future studies will focus on methylated MGMT glioblastoma patients and the ability for the radiologist to visually determine analogs to the radiomic characteristics that we have identified.

MATERIALS AND METHODS

Patients

A total of 159 patients with untreated GBM and pre-treatment MRI scans were retrospectively analyzed in our study, from Memorial Sloan Kettering Cancer Center (MSK=112) and from The Cancer Imaging Archive (TCIA=47). All patients were scanned with MRI consisting of T1WI, T1WICE, and FLAIR images. Ninety-eight patients from MSK were included in the training set for our study; 78 underwent surgery followed by concomitant and adjuvant radiochemotherapy (30 fractions of 6000cGy and temozolomide) according to the standard Stupp protocol [3] and the remaining 20 others patients were treated with hypofractionated radiation (15 fractions of 4150±610cGy), Optune tumor treating fields (Novocure, St. Helier, Jersey) and/or bevacizumab (Genetech, San Francisco, CA). After determining the training cohort for our study, 61 patients (14 MSK + 47 TCIA) with brain MRIs and unknown treatments (for the TCIA) were used to build an independent validation cohort. The analyzed outcome was overall survival (OS) defined as defined as the time between the date of diagnosis and the last follow-up date or the date of event (e.g., death). Radiomics analysis was performed in both the training and validation cohorts. To determine if the combined performance of radiomics and MGMT status would better predict prognosis, we used a reduced cohort of 121 patients (86 for training and 35 for validation) with known MGMT status (see Table 1).

Memorial Sloan Kettering Cancer Center Patients

This retrospective study was approved by the MSK institutional review board with a waiver of written informed consent. MR images were acquired using Discovery MR750w and MR450w scanners (GE Medical Systems, Milwaukee, WI) with a standard head or neurovascular coil. The image acquisition details for this cohort are provided in Supplementary Table 3. All patients underwent maximal safe resection as per standard of care. Tumor DNA was tested using a polymerase chain reaction (PCR) based assay for MGMT promoter hypermethylation in region +28 to +47 from the translation start site in exon 1.

TCIA Patients

From the GBM patients in the TCIA database (n=262), we excluded those patients who (a) were missing at least one of the three sequences, (b) had images with severe motion artifacts, and (c) were either recurrent GBMs or determined to be post-operative MRIs. This resulted in a total of 47 patients between 1998 to 2011 with pre-operative MRIs consisting of T1WI, T1WICE and FLAIR images (included patient IDs are in the supplementary material). The image acquisition details for this cohort are provided in Supplementary Table 1. MGMT methylation status of these patients was obtained from the TCGA database using the corresponding patient IDs. Further details of the analyzed patients from TCIA data can be downloaded on the TCIA website: https://wiki.cancerimagingarchive.net.

Radiomics features extraction

All tumors were segmented using grow-cut segmentation combined with active machine learning [37]. The segmentations were validated by an experienced neuroradiologist with >18 years of experience with MRI and brain tumor imaging (R.J.Y) and manually adjusted where required prior to radiomics feature extraction (Figure 2B).

Tumor segmentation

The algorithm from Veeraraghavan and Miller 2011 [37] was modified to simultaneously process multiple MR image sequences consisting of T1WI, T1WICE and FLAIR images. Prior to segmentation, the image sequences were co-registered with a rigid registration using our in-house implementation using a C++ wrapper function implemented around the Insight ToolKit [38].

Next, starting from user-drawn brush strokes to identify the tumor and background regions, the algorithm generated segmentation of the tumor by automatically combining image intensities from three sequences consisting of FLAIR, T1w pre and T1w post-contrast images. To generate the segmentation, the algorithm learned a model of the tumor through brush stroke inputs using a support vector machine (SVM) algorithm, following which any regions that were difficult to segment were automatically identified by the algorithm and presented as queries to the user to indicate whether those regions corresponded to tumor or background. Upon user input in some of those regions, the algorithm automatically refined the tumor segmentation.

Radiomic features were extracted using the CERR software [39], which was demonstrated to be compliant with the Image Biomarkers Standardization Initiative (IBSI) nomenclature [40]. The MR intensities in all images were scaled to Q=128 values using:

With Vmin and Vmax the minimum and maximum intensity value in the tumor segmentation (Figure 2C). Additionally, Gabor filtering was applied to extract edge maps from the images using four orientations (θ=0°, 30°, 45° and 90°) and with a frequency f of 2 and 2√2 [41]. Gabor filtering consists of applying bandpass filtering operation to extract edges in the image. The Gabor edges are directional sensitive filters and are composed of a Gaussian envelope function superimposed on a sinusoidal wave. In two dimensions, Gabor filters can be used to determine the edge response along various orientation (or angles) and scales (or bandwidths) of the Gaussian function. Using four orientations and two bandwidths on the three different MRI sequences resulted in twenty-four Gabor edge images.

A total of 286 MRI-radiomics texture, shape, edge and intensity histogram features were extracted from the unscaled, scaled and Gabor edge filtered MRIs from all three sequences (Figure 2D). The full list of extracted features is provided in the supplementary file. First-order features consisted of intensity histogram features (n=20) computed from the unscaled MRIs. The scaled MRI images were used to compute GLCM [42] and GLSZM [43]. An offset of one voxel in all directions and 26-connexity was used to extract the co-occurrence matrices. One matrix was built to describe the entire tumor volume using 3D neighbors and one unique matrix for all directions (using the same merging method as described in the IBSI [40]), and then used to extract GLCM features [42]. Similarly, one GLSZM was produced and then used for extracting tumor heterogeneity features [43]. In addition, four main histogram metrics (mean, standard deviation, skewness and kurtosis) were extracted from the Gabor edge maps.

Statistical analysis

Our study was designed in accordance with the guidelines proposed by Vallières et al. for responsible radiomics research for faster clinical translation [44]: we focused particularly 1) to describe precisely the workflow from the image acquisition to the feature extraction and 2) to design a study with an independent validation and with corrections for multiple testing comparison.

All statistical analyses were performed using R software (v3.4.3) and the glmnet, roc and survival packages. Feature selection was performed implicitly using the Least Absolute Shrinkage and Selection Operator (LASSO: elastic net logistic regression with the mixing parameter equal to 1) [45] regression model with the patient's dead or alive (56% vs 42% of patients) status as a binomial response variable and a 10-fold cross-validation was used for the tuning parameter (λ). Prior to fitting the model, the features were transformed to be in the same scale using the standardization method integrated in the glmnet package. The features found to be relevant by the model were those that were assigned non-zero weights by the model. Selected features were then used to compute their association with OS. Kaplan-Meier analysis was used to determine the association between the radiomics features and OS by dichotomizing the patients using the best threshold (Youden index) obtained from ROC curves analysis of the individual features. The p-values were corrected for multiple comparisons using the Bonferroni method for controlling the false discovery rate. Only p-values <0.05 after adjustment for multiple comparisons were considered significant. We also combined the relevant radiomics features identified using the LASSO method with the MGMT molecular biomarker to determine whether the combined model produced better discrimination between patients by overall survival.

An independent validation set was used to validate the results on radiomics and on the combination of radiomics with MGMT for OS prediction. In this cohort, the previously selected radiomics features were dichotomized using the median value reported in the discovery set.