Radiogenomic-Based Survival Risk Stratification of Tumor Habitat on Gd-T1w MRI Is Associated with Biological Processes in Glioblastoma

Glioblastoma (GBM) is a highly aggressive, rapidly fatal, and the most common primary malignant tumor in the brain. Current standard-of-care treatment for GBM is maximal safe surgical resection followed by Stupp protocol—cranial radiotherapy and concomitant chemotherapy, with temozolomide (TMZ; ref. 1). Despite this aggressive and multimodal treatment, median survival has only slightly improved to approximately 15 months, with less than 25% of patients surviving up to 2 years, and less than 10% surviving for over 5 years (2). This poor prognosis can be attributed to genetic instability and intratumor and intertumor heterogeneity of GBM that leads to treatment resistance, progression, and tumor recurrence (3, 4).

Recent investigative studies have identified several driver mutations [i.e., isocitrate dehydrogenase (IDH)], chromosomal anomalies (i.e., 1p19q chromosome deletions; ref. 5), epigenetic alterations [such as O-6-methylguanine-DNA dethyltransferase (MGMT); ref. 6] as specific targets for personalizing anticancer therapeutics, enhancing treatment response, and improving survival. Although these genomic biomarkers have shown some prognostic value in GBM, biopsies are tissue destructive, rely on a small sample of tissue that do not account for the inherent intratumoral heterogeneity extant in GBMs, and are usually expensive and time-consuming to perform during routine screening. With almost 40% of GBM patients diagnosed with cancer recurrence within 6 months even after aggressive treatment management, alternate treatment options such as antiangiogenesis therapies, immunotherapy vaccinations, and clinical trials could be more suited in patients who may be at high risk of tumor recurrence, to achieve prolonged progression-free survival (PFS; refs. 7–9). There is, hence, a need to develop imaging-based biomarkers that are prognostic of patient outcome to potentially assist in building comprehensive and patient-centric treatment plans.

Multiple studies suggest that tumor-associated cellular alterations are not just limited to the visible tumor enhancement of GBM, but infiltrate beyond the tumor margins into the peritumoral brain zone (10). It is also known that 90% of GBM tumors recur in the peritumoral region and also the inner necrotic compartments (11). Radiomics has provided a surrogate mechanism to noninvasively characterize the tumor by capturing subvisual cues of morphologic diversity (e.g., roughness, image homogeneity, regularity, and edges) on routine MRI scans for disease diagnosis, patient prognosis, and evaluate treatment response (12). Recently, a few research groups have demonstrated the added prognostic value in investigating radiomic features, not just from within these visible tumor margins, but also from the surrounding peritumoral regions (i.e., edema) as well as the features from necrotic core (13–15). This concept of “tumor habitat” involves interrogating the radiomic features from the enhancing tumor as well as the necrotic core, and peritumoral edema subcompartments across different MRI protocols.

Although a few recent studies have attempted to develop prognostic radiomic-based models from the tumor habitat to predict progression-free survival (PFS) and overall survival (OS) in GBM, there has been limited work on creating a reliable survival risk assessment model using radiomic features from the tumor habitat, using a large multi-institutional cohort. Further, a more critical missing link in previous work has been on identifying the molecular associations of habitat-specific radiomic features with the underlying signaling pathways that drive different biological processes via radiogenomic analysis. Improved understanding of how the changes in biological processes at the molecular level affect changes at a radiologic scale could help drive adoption of the prognostic radiomic features in a clinical setting. Further, such cross-scale associations of radiomic-based prognostic markers with molecular processes in GBM could allow for designing patient-centric treatments based on their risk profile for poor survival.

In this work, we have two objectives. First, we seek to create a prognostic radiomic risk score (RRS) using radiomic features from different subcompartments of the tumor habitat on routine multiparametric MRI protocols (i.e., Gd-T1w, T2w, T2w-FLAIR) to risk-stratify GBM patients based on their PFS. The RRS model will be trained using a multi-institutional training cohort of 130 patients via least absolute shrinkage and selection operator (LASSO) Cox regression model. We will then evaluate the ability of this RRS to stratify an independent test cohort (n = 73) into two groups of “low risk” and “high risk” of poor survival, based on the fixed median cutoff value of the RRS. In our second objective, in an attempt to provide a biological basis for the radiomic features driving the RRS, we will seek to identify radiogenomic correlations between the most prognostic radiomic features from different subcompartments of the tumor habitat with molecular signaling pathway networks in GBM using Gene Ontology (GO) and single-sample gene set enrichment analysis (ssGSEA; refs. 16–18). Figure 1B illustrates an overview of our framework.

Our retrospective cohort consisted of treatment-naïve multiparametric MRI scans from three cohorts. Two of the three data sets were curated from publicly available data sets: The Cancer Imaging Archive (TCIA; refs. 19, 20) and The Ivy Glioblastoma Atlas Project (Ivy GAP; refs. 21, 22). TCIA is an open archive of cancer-specific medical images and associated clinical metadata established by the collaboration between the NCI and participating institutions in the United States. Similarly, Ivy GAP is freely accessible online data resource that aims to identify important morphologic hallmarks of GBM.

The third data set was curated from a Health Insurance Portability and Accountability Act compliant and institution review board (IRB)–approved participating institution—Cleveland Clinic Foundation (CCF), where the need for an informed consent from all patients was waived. Between December 1, 2011, and May 1, 2018, radiology image archives of this participating institution were searched to identify 80 GBM patients as confirmed via histopathology. Further information on the images acquired can be found in Supplementary Section S1.

The retrospectively collected 369 GBM cases from three sites were further triaged using inclusion criteria that involved the availability of (i) routine MRI sequences (Gd-T1w, T2w, and T2w-FLAIR) for treatment-naïve patients with diagnostic image quality and (ii) PFS information for all individuals. A total of 166 patients with the following criteria were excluded: (i) absence of baseline/pretreatment surgery scans (n = 47), (ii) missing either Gd-T1w, T2w, or T2w-FLAIR MRI protocol scans (n = 76), (iii) presence of MRI artifacts (n = 27), and (iv) missing PFS information (n = 16). Following the patient exclusion criteria, 203 patients were enrolled for the radiomic analysis. The patient enrolment flowchart is shown in Fig. 1A.

Of the 203 GBM patients who followed the inclusion/exclusion criteria for radiomic analysis, 130 patients from TCIA cohort comprised our training cohort, and the remaining 73 cases were used as independent test cohort. Further, from the 203 GBM cases, 157 tumors were selected for radiogenomic analysis after applying an additional inclusion criterion that required RNA-sequencing data obtained from cellular tumor. A total of 32 GBM tumors were further excluded from the study due to the absence of corresponding gene-expression data acquired on Affymetrix U133A platform. Finally, a total of 125 GBM tumors were investigated for the radiogenomic analysis.

Multisite data distribution has been reported in Table 1A. Our data selection for training and testing was performed in a randomized fashion, while ensuring that the distribution of patient demographics such as clinical factors (patient's age at diagnosis, gender) and molecular characteristics (MGMT promoter methylation status, IDH status) was balanced between the training and test sets (n = 130 in the training and n = 73 in the independent test set). The clinical characteristics and treatment follow-up of the patients within the training and test cohort have been summarized in Table 1B. Within the training cohort, the extent of resection [gross total resection (GTR) or near total resection (NTR) or subtotal resection (STR)] was not known in 106 cases. From the test cohort (n = 73), 40 GBM cases had undergone GTR, 12 GBM cases underwent NTR, and 12 GBM tumors underwent STR. Five cases within the test cohort underwent laser ablation.

No statistically significant difference was found between the training and test cohorts with respect to the age (P = 0.75), gender (P = 0.45), and PFS (P = 0.37).

Preprocessed TCIA cases along with annotations on every 2D slice were obtained from Bakas and colleagues (23, 24), Cases obtained from Ivy GAP and CCF were preprocessed using the same pipeline as the TCIA cohort (25). For every study, MR images were reoriented to the LPS (left–posterior–superior) coordinate system. All MRI protocols were then coregistered to the T1w anatomic template of SRI24 Multi-Channel Normal Adult Human Brain Atlas via affine registration and resampled to 1 mm³ voxel resolution (26). Subsequent skull-stripping was implemented on all cases (25). N4-bias correction algorithm (27) was not implemented in our preprocessing pipeline as it is known to obliterate the T2w-FLAIR signal (23). However, a low-level image processing smoothing filter, Smallest Univalue Segment Assimilating Nucleus (SUSAN), was used to reduce high-frequency intensity variations of noise (28). Lastly, we corrected the MRI protocols for intensity nonstandardness, which refers to the issue of MR image “intensity drift” across different imaging acquisitions. The intensity nonstandardness results in MR image intensities lacking tissue-specific numeric meaning within the same MRI protocol, for the same body region, or for images of the same patient obtained on the same scanner. Therefore, intensity standardization was implemented in MATLAB R2014b (Mathworks) using the method presented in Madabhushi and colleagues (29).

Every GBM tumor was defined into 3 subcompartments of the tumor habitat: (i) tumor necrosis, (ii) enhancing region of the tumor, and (iii) peritumoral edema. Gd-T1w MR images were used to delineate hypointense signals to capture the necrotic core of GBM tumor, usually located in the central region of the tumor. Similarly, enhancing region of GBM is represented as hyperintense regions on Gd-T1w MR images when compared with precontrast T1w. When precontrast T1w scans were not available, T2w images and T2w-FLAIR scans were used to evaluate the presence of blood (30). T2w and T2w-FLAIR scans were also used to identify edema and necrotic core.

The independent test cohort, which consisted of Ivy GAP and CCF cases, was manually annotated using 3D Slicer software (31). Every 2D slice of each MRI scan with visible tumor was manually annotated by the expert readers. A total of three experts were asked to perform the manual annotations on a total of 73 test studies. The senior-most expert (expert 1, V. Hill, >10 years of experience in neuroradiology) independently annotated the studies obtained from Ivy GAP, whereas expert 2 (V. Statsevych) with 7 years of experience in neuroradiology supervised expert 3 (K. Bera, with 3 years of radiology experience), to manually annotate the CCF cases from the test set. In rare cases with disagreement across the two readers (experts 2 and 3), the senior-most radiologist (V. Hill, expert 1) was consulted to obtain the final segmentations.

A total of 936 3D texture-based radiomic and 19 shape-based features were extracted individually from every subcompartment (edema, necrosis, and enhancing tumor) on every MRI protocol. The feature set for every study on each MRI protocol included 19 shape-based features, 52 Haralick features (capturing tumor heterogeneity), 501 Laws energy (captures presence of spots, edges, waves, and ripples in an image), 383 Gabor wavelet (capturing structural detail at different orientations and scales)-based features on a per-voxel basis. Statistics of median, standard deviation, skewness, and kurtosis were then calculated from the feature responses of all voxels within the region of interest. Therefore, a total of 2,850 features extracted for every study on every MRI protocol (Gd-T1w, T2w, and T2w-FLAIR) across all tumor subcompartments. A list of the extracted features is summarized in Supplementary Sheet 1. All feature values were normalized (mean of 0 and standard deviation of 1). Feature extraction and statistic calculations were performed using in-house software implemented in MATLAB R2014b platform (MathWorks).

A list of extracted shape-based features has been provided in Supplementary Section S2 and Supplementary Table S1. Detailed description of the set of features used in this work and its possible relationship to the pathophysiology of GBM is provided in Supplementary Section S2 and Supplementary Table S2. Details of implementation and parameters of every 3D texture feature extracted have been provided in Supplementary Section 2 and Supplementary Table S3. The feature extraction pipeline has been made publicly available at https://github.com/ccipd/BrIC_Lab.

The experimental design was set up to evaluate the prognostic ability of radiomic features from every MRI protocol (Gd-T1w, T2w, T2w-FLAIR). Within the training set, univariable Cox regression was used as the feature pruning method, where all features were investigated to find the most statistically significant prognostic subset radiomic features with a P value less than 0.05. After feature pruning, three LASSO Cox regression models were set up in 5-fold cross validation with 1,000 iterations, to further select the most prognostic radiomic features within each of the MRI protocols based on their frequency of occurrence. LASSO L1 regularization technique iteratively shrinks the feature coefficient estimates toward zero, and chooses an optimal tuning parameter lambda (λ) that increases in a cross-validation setup until features with only nonzero coefficients are selected. Therefore, in cases with a very large number of features, a LASSO model can help both shrink and find the sparse model that involves a small subset of the features. Thus, features picked by the LASSO models within the training cohort and then pooled in a linear combination and multiplied with their respective coefficients to construct an RRS:

where n is the number of features selected by LASSO for a given MRI protocol, β is the weighted coefficient of the selected feature, and τ is the selected MRI protocol–based radiomic feature. Finally, the constructed RRS was then applied in the test set for validation.

For every MRI protocol, each patient within the training and test sets was stratified into high-risk and low-risk groups based on the median value of the RRS as a fixed cutoff. The prognostic ability of the RRS was evaluated by implementing the Kaplan–Meier (KM) survival analysis and Cox regression. The KM survival curves were used to compare survival times across high-risk and low-risk groups within the training and testing cohorts for (i) Gd-T1w MRI protocol, (ii) T2w MRI protocol, and (iii) T2w-FLAIR MRI protocol. The horizontal axis on the KM survival curve represents the time, and the vertical axis shows the probability of survival. At any given point on the survival curve, the probability that a patient in each group to remain alive at that time is presented (13).

Further, hazard ratios (HR) were used to quantify the effect of individual feature on survival. Features yielding negative regression coefficients (i.e., low feature values correlated with long-term survival) in our Cox model produce a HR between 0 and 1; features yielding positive regression coefficients (i.e., low feature values correlated with short-term survival) produce a HR between 1 and infinity. We also computed concordance indices (C indices or C statistic) for each of our univariable and multivariable analysis experiments in R. C indices is the fraction of all pairs of subjects whose predicted survival times are correctly ordered (i.e., concordant with actual survival times). C indices = 1 indicates that the model has perfect predictive accuracy, and C indices = 0.5 indicates that the model is no better than random chance.

Within the training cohort of TCIA patients, preprocessed gene-expression data using Robust Multichip Average (RMA) algorithm (32) were available for a total of 125 patients (Table 1A) and hence used for our radiogenomic analysis. This TCGA GBM gene-expression data (level3, Affymetrix HT HG U133A) is publicly available for download from the Broad Institute (http://gdac.broadinstitute.org/runs/stddata__2016_01_28; ref. 33). Further information regarding the acquisition and expression profiling has been published by TCGA network (34). This transcriptomic data, consisting of 12,042 annotated genes, was used to investigate the possible underlying biological processes of the RRS. The most differentially expressing genes (DEG) of the RRS were selected using Wilcoxon rank sum test. The Benjamini and Hochberg method was used to adjust the P values and control for the false discovery rate (FDR) in multiple testing (35). Statistically significant DEGs (P < 0.03) were then used to identify distinct GO-based biological processes (16, 17). GO highlights the most overrepresented genes and finds the systematic linkages between those genes and biological processes.

To gain further insight into the GO-based biological processes and their association with the individual radiomic features from different subcompartments that were used to create the Gd-T1w RRS, we performed ssGSEA (18, 36). For a predefined set of genes, ssGSEA captures the significantly enriched or depleted biological processes and calculates an enrichment score for every patient in the cohort. These predefined set of genes for the GO-based biological processes were acquired from the Molecular Signatures Database (MSigDB v6.2) platform (http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection= BP; ref. 37). Within the training cohort, ssGSEA score was calculated for 125 GBM patients, and then a Spearman rank correlation coefficient (ρ) matrix was constructed to identify the contribution of the 25 individual Gd-T1w radiomic features (which were used to create the RRS) and their association with the GO biological processes.

All the statistical analyses were performed on the MATLAB R2014b platform (MathWorks). P values <0.05 were considered significant. The Cox proportional hazard regression models were built using the “coxphfit” command in MATLAB.

Within the training cohort, 212 prognostic radiomic features were obtained from 2,850 Gd-T1w radiomic features after feature pruning from univariable Cox regression. For the Gd-T1w protocol, the LASSO model selected 25 radiomic features with lambda value of 0.088 (see Fig. 2A). Details of the features selected and their coefficients have been listed in Fig. 2B. The complete RSS formulation can be found in Supplementary Section S3. Of the 25 Gd-T1w radiomic features, 8 were picked from the peritumoral edema region, 9 were selected from the enhancing region of the tumor, and the remaining 8 were picked from the necrotic core of GBM. The Gd-T1w RRS for every patient within the training and testing cohorts has been provided in Supplementary Sheet S2.

In a multivariable Cox analysis, the RRS built using the 25 Gd-T1w radiomic features resulted in statistically significant KM curves (log-rank test, P < 0.00001, HR = 4.5, Fig. 2C, left). The Gd-T1w RRS was also found to be a prognostic indicator of PFS in the independent test cohort (log-rank test, P = 0.0117, HR = 2.0735; Fig. 2C, right). Other clinical parameters (such as age and gender) between the low-risk and high-risk groups based on Gd-T1w RRS median cutoff were also investigated (Supplementary Section S4 and Supplementary Table S3).

Similarly, 321 and 496 prognostic features were obtained from 2,850 T2w and T2w-FLAIR MRI-based radiomic features, respectively, using the univariable Cox regression method. Then, a total of 25 T2w and 12 T2w-FLAIR radiomic features were picked by the LASSO models, respectively. Similar to Gd-T1w RRS, statistically significant KM survival curves were reported with 25 T2w (log-rank test, P < 0.00001, HR = 28.6736) and 12 T2w-FLAIR (log-rank test, P value = 0.0013, HR = 4.8434) MRI radiomic features, respectively. However, the RRS created using radiomic features from the tumor habitat of T2w and T2w-FLAIR studies did not yield statistically significant differences in the “high-risk” and “low-risk” patient populations upon validation (log-rank test, T2w P = 0.49, T2w-FLAIR P = 0.25).

We further assessed the added clinical utility of the RRS to existing clinical parameters (age and gender), molecular features (MGMT and IDH mutation status), and extent of resection (EOR; on test set only) to evaluate improvement in prediction of PFS in GBM patients. In a multivariable setting, we found that the Gd-T1w RRS resulted in a training concordance index (C-index) of 0.76 [95% confidence interval (CI): 0.67–0.85] and predicted PFS on the independent test cohort with a C-index of 0.75 (95% CI, 0.68–0.84, P = 0.012). Interestingly, combining clinical (age, gender), molecular features (MGMT, IDH mutation status) and EOR with the Gd-T1w RRS, improved the prediction of PFS within the training (C-index = 0.81; 95% CI, 0.71–0.90, P value = <0.0001) and test set (C-index = 0.84; 95% CI, 0.67–0.95, P value = 0.03), respectively, as compared with validation on test set using either clinical (age: C-index = 0.55; 95% CI, 0.46–0.64; gender: C-index = 0.52; 95% CI, 0.46–0.64) or molecular features (MGMT status: C-index = 0.65; 95% CI, 0.57–0.72; IDH status: C-index = 0.55; 95% CI, 0.49–0.60) or EOR (C-index = 0.56; 95% CI, 0.44–0.67) or Gd-T1w RRS (C-index = 0.75; 95% CI, 0.68–0.84) alone. Table 2 lists the HR and concordance indices for the clinical parameters (age and gender), molecular features (MGMT status and IDH status), EOR, and combined radiomic features for predicting PFS.

Empirical analysis of the 12,042 annotated genes (obtained from 125 patients of the training cohort) using the Wilcoxon rank sum test resulted in a total of 192 DEGs that had a P value of less than 0.03 (with an FDR of 3%). These 192 DEGs were identified to be differentially expressed between the “high-risk” and “low-risk” groups of the Gd-T1w RRS and were further used for GO analysis. Figure 3A shows the supervised hierarchical clustering of these 192 DEGs using the correlation distance metric. Figure 3B shows the boxplots of two representative DEGs “ID1” (inhibitor of differentiation/DNA binding) and “BMP4” (bone morphogenetic protein 4). A complete list of the 192 DEGs can be found in Supplementary Sheet S3.

Radiogenomic analysis of the training cohort using GO revealed that the Gd-T1w RRS was associated with 57 biological processes that have been listed in Fig. 3C. A complete list of all the biological processes has been provided in Supplementary Sheet S4. It was observed that 24 biological processes were implicated in cell adhesion, cell proliferation, differentiation, and angiogenesis. A directed acyclic graph investigating the interrelationships between these 57 biological processes has been provided in Supplementary Section S5 (38). In Fig. 3D, the fold enrichment change and number of genes involved in the various biological processes have been listed. It can be seen that biological processes of cell adhesion, which involved over 100 genes, had the highest fold-enrichment change. Most of the GO terms involved in cell differentiation, proliferation, and angiogenesis had less than 100 overrepresented genes, but their fold enrichment change was greater than 2.

Predefined gene set annotations for 15 GO-based biological processes (associated with cell differentiation, cell adhesion, and angiogenesis) were available from MSigDB and thus used to calculate the ssGSEA scores. Supplementary Sheet S5 has the complete list of predefined set of genes that were used to describe these 15 GO-based biological processes. ssGSEA scores helped understand the associations of these biological processes with the individual 25 radiomic features that were obtained from different tumor subcompartments and used to create the Gd-T1w RRS. As may be observed in Fig. 4, statistically significant correlations (P < 0.05) were found between the shape features (i.e., sphericity, elongation, and convexity) of the peritumoral edema region and biological processes of cell proliferation, angiogenesis, and cell adhesion. Additionally, Laws energy and Gabor wavelet texture features from within the peritumoral edema region of the GBM tumor habitat were also found to be correlated with cell proliferation (ρ_{Laws Energy} = −0.186; 95% CI, −0.350 to −0.011) and angiogenesis (ρ_Gabor = 0.193; 95% CI, 0.018–0.356), respectively.

Pretreatment estimation of PFS can help predict therapeutic efficacy of conventional first-line treatment. Additionally, recent randomized clinical trials have reported that PFS can be used as a surrogate endpoint that correlates well with OS (39, 40). Currently, there are no clinically validated imaging-based markers that are prognostic of PFS in primary Glioblastoma (GBM; ref. 41).

In this work, we had two objectives. In our first objective, we created and independently evaluated an RRS using radiomic features from the tumor habitat on pretreatment MR imaging to stratify GBM patients based on their PFS on 203 studies from a multi-institutional cohort. In our second objective, we reported significant associations of the prognostic radiomic features that contributed to RRS with specific molecular signaling pathways, via radiogenomic analysis.

In our first objective, the survival risk score created using radiomic features from Gd-T1w MRI was found to be statistically significantly different across the “low-risk” and “high-risk” groups both on the training (P < 0.001, n = 130) and the holdout test sets (P = 0.03, n = 73). In consensus with recent studies (42–44), we further discovered that integrating the risk score from Gd-T1w MRI with clinical parameters (age and gender) and molecular features (MGMT status and IDH status) improved the performance of the prognostic model. Our prognostic RRS consisted of a total of 25 radiomic features including textural features belonging to Laws energy and Gabor wavelet families, and as shape features from the peritumoral edema region, on Gd-T1w MRI. Laws energy radiomic features capture wavy, ripple, and spot-like patterns, which may be reflective of increased peritumoral neovascularization and atypical blood vessels. Our findings also corroborate with our previous work, where we identified Laws energy features from the edematous and enhancing region on T2w-FLAIR and Gd-T1w MRI to be predictive of OS in GBM (13, 43). Similarly, Gabor wavelet features capture changes in structural orientation by quantifying prominent MRI intensity changes across multiple directional gradients in a given region of interest (ROI; ref. 45). These findings are in line with multiple studies that have demonstrated the efficacy of Gabor wavelets in predicting tumor malignancy (46), disease aggressiveness (47), and treatment response (48) on imaging for various cancers.

In our second objective, we identified significant radiogenomic correlations (P < 0.05) between our 25 prognostic Gd-T1w radiomic features from different subcompartments of the tumor habitat with molecular signaling pathways in GBM using GO and ssGSEA. Specifically, using the corresponding gene-expression data available for the training set (n = 125), we identified a total of 192 DEGs across the “low-risk” and “high-risk” groups obtained from the Gd-T1w RRS. Prognostic value of some of these genes has been demonstrated within various cancers. For example, ID1 is known to be overexpressed in breast and GBM cancers (49, 50). In consensus, we found that ID1 had an adjusted P value of 0.0039. Similarly, high expression of BMP4 has been known to be significantly associated with better prognosis in GBM (51). Our analysis revealed that the “low-risk” group was associated with statistically significant BMP4 expression when compared with the “high-risk” group, within the training set (P = 0.0069). TP53 (tumor protein p-53)-inducible proteins such as TP53I11 (tumor protein p53 inducible protein 11) and TP73-AS1 (TP73 antisense RNA 1) were also found to be statistically significant between the “low-risk” and “high-risk” radiomic risk groups with P values of 0.005 and 0.007, respectively. CD34, a transmembrane phosphoglycoprotein, known to be a prognostic marker with potential to be a therapeutic target, was differentially expressed with an adjusted P value of 0.017 (52).

Using the 192 DEGs, GO analysis further revealed that cell adhesion, cell differentiation, proliferation, and angiogenesis were the key biological processes that were significantly different across the high-risk and low-risk groups (as defined via our RRS). Biological processes surrounding cell adhesion, cell differentiation, proliferation, and angiogenesis are key contributors in driving malignant tumor behavior, support resistance to chemoradiotherapy, and contribute to poor PFS in primary GBM (53, 54). Previous work by Liu and colleagues has similarly shown via radiogenomic analysis that biological processes such as immune response, programmed cell death, cell proliferation, and vasculature development dictate PFS in low-grade gliomas (42).

Finally, we demonstrated that multiple biological processes that are associated with poor PFS were also significantly correlated (P < 0.05) with radiomic features extracted from the peritumoral edema region of GBM. For instance, we found that sphericity (a measure of 3D compactness) of the edematous region was negatively correlated with most of the biological processes involved in cell adhesion (ρ_avg = −0.140; 95% CI, −0.308–0.036) and angiogenesis (ρ_avg = −0.178; 95% CI, −0.343 to −0.002). This suggests that a more infiltrative pattern of tumor growth as reflected on peritumoral edema may be associated with poor prognosis in GBM patients. Further, Ismail and colleagues have previously reported that lower values of elongation of the edematous region (ratio between major and minor axes of the tumor) are indicative of tumor progression in conventional posttreatment primary GBM (55). Along similar lines, we found association of lower elongation values of edema with angiogenesis (ρ_avg = 0.118; 95% CI, −0.287 to 0.058) and cell differentiation (ρ_avg = 0.128; 95% CI, −0.048 to 0.299), which suggests that better prognosis may be linked with low RRS in GBM.

This study did have some limitations. The validation was done retrospectively, and radiogenomic analysis was possible only on 1 of the 3 cohorts (due to unavailability of the corresponding gene-expression data in the Cleveland Clinic cohort and substantial heterogeneity in RNA-sequencing processing pipeline of Ivy GAP cohort). Although we compared the clinical parameters (age and gender), molecular (MGMT and IDH status), and EOR with RRS to the extent possible, due to the limited training cohort, the development of RRS in this work could not be controlled for clinical parameters, molecular status, and EOR. This will be part of our future study involving large multisite evaluation of RRS. Further work is also mandated to develop a radiogenomic approach that incorporates transcriptomic data coming from different processing platforms, and thus reducing the influence of batch effects in RNA-sequencing data. In the future, our work will also focus on incorporating complementary imaging parameters (PET, perfusion, and DWI), which may further improve survival prediction using radiogenomic analysis, while also accounting for follow-up treatment. Additionally, we also plan to expand our analysis with corresponding gene-expression data available across all patients and eventually using prospectively collected scans.

In summary, we developed and independently evaluated a prognostic RRS using features obtained from the tumor habitat on Gd-T1w MRI, which could potentially be used as a noninvasive imaging-based surrogate marker of PFS in primary GBM. We further demonstrated that the Gd-T1w RRS, with other clinical and molecular features, can improve patient stratification when compared with using Gd-T1w RRS alone. Further, we identified the biological processes of cell differentiation, cell adhesion, and angiogenesis to be associated with Laws energy features from the enhancing tumor and shape features of infiltrative peritumoral edema, respectively, thus providing a biological basis for the radiomic features driving the RRS with possible underlying processes that dictate tumor behavior in GBM patients. Following additional large-scale prospective validation, our Gd-T1w RRS could potentially be used as a surrogate endpoint in clinical trials (39), predict onset of tumor recurrence (40, 41), and for evaluating the efficacy of chemoradiotherapy in primary GBM patients (56, 57).

V. Varadan reports a paid consultant/advisory board relationship with Curis Inc. and commercial research grants from Curis Inc. and Philips Healthcare. A. Madabhushi reports commercial research grants from Philips and Inspirata Inc, ownership interest (including patents) in Euclid Bioimaging and Inspirata Inc., unpaid consultant/advisory board relationships with Inspirata Inc., AstraZeneca, Bristol Meyers-Squibb, and Merck, and involvement in an NIH U24 grant with PathCore Inc. and in 3 different R01 grants with Inspirata Inc. N. Braman is an employee/paid consultant for IBM Research. V.B. Hill is an employee/paid consultant for Google. M. Ahluwalia is an employee/paid consultant for AbbVie, AstraZeneca, BMS, Monteris, Bayer, Karyopharm, Forma Therapeutics, Varian Medical Systems, Flatiron, Tocagen, CBT Pharmaceuticals and VBI Vaccines, reports receiving commercial research grants from Abbvie, AstraZeneca, Bayer, Bristol-Myers Squibb, Merck, Incyte, Novocure, Novartis, Pharmacyclics, Inspire, and holds ownership interest (including patents) in Mimivax and Doctible. No potential conflicts of interest were disclosed by the other authors.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH, the U.S. Department of Veterans Affairs, the Department of Defense, or the United States Government.

Conception and design: N. Beig, K. Bera, P. Prasanna, M.S. Ahluwalia, A. Madabhushi, P. Tiwari

Development of methodology: N. Beig, K. Bera, P. Prasanna, J. Antunes, N. Braman, R. Verma, M.S. Ahluwalia, A. Madabhushi, P. Tiwari

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): N. Beig, K. Bera, A.S. Bamashmos, R. Verma, V.B. Hill, V. Statsevych, M.S. Ahluwalia

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N. Beig, K. Bera, P. Prasanna, R. Correa, S. Singh, A.S. Bamashmos, M.S. Ahluwalia, V. Varadan, A. Madabhushi, P. Tiwari

Writing, review, and/or revision of the manuscript: N. Beig, K. Bera, P. Prasanna, J. Antunes, S. Singh, N. Braman, R. Verma, V.B. Hill, M.S. Ahluwalia, V. Varadan, A. Madabhushi, P. Tiwari

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N. Beig, R. Correa, M. Ismail, R. Verma, V.B. Hill, V. Statsevych

Study supervision: P. Prasanna, A. Madabhushi, P. Tiwari

Research reported in this publication was supported by the NCI of the NIH under award numbers 1U24CA199374-01, R01CA202752-01A1, R01CA208236-01A1, R01 CA216579-01A1, R01 CA220581-01A1, 1U01 CA239055-01, and 1P20 CA233216-01; National Center for Research Resources under award number 1 C06 RR12463-01; VA Merit Review Award IBX004121A from the United States Department of Veterans Affairs Biomedical Laboratory Research and Development Service; the DOD Prostate Cancer Idea Development Award (W81XWH-15-1-0558); the DOD Lung Cancer Investigator-Initiated Translational Research Award (W81XWH-18-1-0440); the DOD Peer Reviewed Cancer Research Program (W81XWH-16-1-0329); National Institute of Diabetes and Digestive and Kidney Diseases (1K25 DK115904-01A1); the Ohio Third Frontier Technology Validation Fund; the Wallace H. Coulter Foundation Program in the Department of Biomedical Engineering and the Clinical and Translational Science Award Program (CTSA) at Case Western Reserve University; Department of Defense Peer Reviewed Cancer Research Program (PRCRP) Career Development Award; Dana Foundation David Mahoney Neuroimaging Program; and The V Foundation Translational Grant.

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