Survival prediction in gliomas based on MRI radiomics combined with clinical factors and molecular biomarkers

Patients

The patient data used for this study was approved by the Institutional Review Board of First Hospital of Shanxi Medical University, with the approval number: 2021 K-K073. For another TCGA/TCIA database, since all information on the patients was accessible for public download, no institutional review board permission was required.

Public datasets can be downloaded from The Cancer Imaging Archive (TCIA) official website (https://www.cancerimagingarchive.net/collections/). A total of 709 participants from the TCGA/TCIA project were enrolled, including TCGA-LGG, TCGA-GBM, and UCSF-PDGM datasets as the training cohort and internal validation cohort. Additionally, 326 patients with pathologically diagnosed grade 2–4 glioma were retrospectively collected between October 2011 and April 2019 from the First Hospital of Shanxi Medical University (FHSXMU) and Shanxi Provincial People’s Hospital (SPPH) as an external validation cohort. All patients determined in this study met the following criteria: (a) pathological confirmation of Grade 2–4 glioma post-surgery; (b) availability of complete preoperative MRI data, including CE-T1WI and T2FLAIR images; (c) determination of tumor molecular markers IDH and MGMT; (d) a follow-up time exceeded 2 years or endpoint events occurred. The data for the 414 patients from TCGA/TCIA and the 165 patients from our hospital included in the study are shown in Table S1. Overall survival was calculated from the moment of postoperative pathological diagnosis to death or the last follow-up if patients were still alive. The detailed patient inclusion process is shown in Fig. 1.

MRI protocol

Preoperative MRI images were acquired on a 3.0T scanner (Signa HDxt, GE Healthcare, Waukesha, WI, USA) equipped with an eight-channel matrix coil. The sequence parameters included CE-T1WI (TR/TE, 195/4.76 ms; FOV, 240 × 240 mm²; slice thickness/spacing, 5.0/1.5 mm; matrix, 256 × 256) and T2FLAIR (TR/TE, 8,000/95 ms; TI, 2,000 ms; FOV, 240 × 240 mm²; slice thickness/spacing, 5.0/1.5 mm; matrix, 256 × 256) images. The CE-T1WI was carried out after the injection of 0.1 mmol/kg gadolinium contrast agent.

IDH genotyping and MGMT promoter methylation detection

For patients in the TCGA/TCIA queue, IDH mutation status and MGMT methylation data were downloaded from TCGA and cBioPortal (http://www.cbioportal.org/). Institutional specimens underwent deparaffinization followed by microdissection to ensure ≥80% tumor cellularity. DNA was extracted using the Simlex OUP® FFPE DNA kit. IDH1/2 mutations were analyzed by Sanger sequencing after PCR amplification (ABI9700) using specific primers: IDH1-F: 5′-CGGTCTTCAGAGAAGCCATT-3′ and R: 5′-GCAAAATCACATTATTGCCAA-3′, and for IDH2-F: 5′-CAAGAGGATGGCTAGGCGAG-3′ and R: 5′-CAAGCTGAAGAAGATGTGGAAAAG-3′.

Pyrosequencing was used to assess MGMT methylationb (Quillien et al., 2014). The extracted DNA was modified with bisulfite using the BisulFlashTM DNA modification kit (Epigentek, Farmingdale, NY, USA). PCR amplification was performed using the DRR006 kit (Takara, Shiga, Japan) with a reaction volume of 40 μl. PCR amplification (DRR006 kit, 40 μl) conditions: 94 °C/2 min; 50 cycles of 94 °C/20 s, 55 °C/20 s, 72 °C/20 s; final extension 72 °C/5 min. Pyrosequencing of the PCR products was performed on the PyroMark Q96 system (Qiagen, Hilden, Germany) according to the provided instructions. Ten CpG sites were quantitatively analyzed within the MGMT promoter, with methylation thresholds defined as ≥8% (methylated) and <8% (unmethylated) based on averaged values.

Image preprocessing and tumor ROI segmentation

UCSF-PDGM MRI data were preprocessed using an open-source deep-learning tool to remove skull structures. Images were aligned and resampled to 1 mm isotropic resolution in T2FLAIR space through nonlinear registration. Tumor regions (enhancing tumor, necrosis, and edema) were automatically segmented using Brain Tumor Segmentation (BraTS) Challenge-based algorithms. Detail the segmentation methodologies applied to the UCSF-PDGM dataset in the PDF referenced by the DOI: https://doi.org/10.48550/arXiv.1811.02629. The code for brain tumor segmentation algorithms based on the BraTS Challenge can be found in the GitHub project brats-unet: UNet for brain tumor segmentation (https://github.com/icerain-alt/brats-unet). We selected the contrast-enhanced tumor regions for feature extraction and subsequent analysis. The preprocessing pipeline for the CE-T1WI and T2FLAIR images from TCGA-LGG, TCGA-GBM, and our institutional data involved spatial resolution normalization through resampling to isotropic 1 mm³ voxel dimensions. We performed z-score normalization on the image grayscale values, resulting in a mean of 0 and a standard deviation of 1 (Carré et al., 2020). The T2FLAIR images were co-registered with the corresponding CE-T1WI images by affine transformation using the Oxford Centre for Functional MRI of the Brain (FMRIB) Linear Image Registration Tool (FLIRT) from the FMRIB Software Library (FSL). All CE-T1WI and T2FLAIR MRI images underwent resampling and normalization using Matlab. Tumors were manually segmented on the axial thin-section 2D slices by two radiologists with 10 and 15 years of experience using ITK-SNAP software (http://www.itksnap.org). For enhanced tumors, the region of interest (ROI) was described based on the enhanced edge on CE-T1WI images. For unenhanced tumors, the tumor signal on CE-T1WI images was lower compared to the peritumoral edema area, and the ROI was determined based on the scope of unenhanced lesions on CE-T1WI images. The final ROI was verified by a senior radiologist with 20 years of experience. The intraclass and interclass correlation coefficients (ICC) were utilized to assess the stability and repeatability of radiomic features. To evaluate intergroup consistency, the first two radiologists independently depicted 30 arbitrary gliomas. Two weeks later, the first physician re-delineated the same 30 cases to estimate intragroup consistency. The flowchart of this study is illustrated in Fig. 2.

Radiomics feature extraction and selection

A total of 1,781 radiomics features were independently extracted from the tumor regions of CE-T1WI and T2FLAIR MRI images using the open-source software FAE (https://github.com/salan668/FAE) based on PyRadiomics. These features included 14 shape features, 18 first-order features, 75 texture features (consisting of 24 GLCM features, 14 GLDM features, 16 GLRLM features, 16 GLSZM features, and five NGTDM features) of the original images, 144 first-order features and 600 texture features of wavelet transform, 18 first-order features and 75 texture features drawn on square, square root, gradient, laplace, exponent, logarithm, 2D transformation, and 279 first-order and texture features based on 3D transformation. Features with an ICC higher than 0.75 were used for subsequent analysis. Z-score was used to normalize the CE-T1WI and T2FLAIR feature matrices before feature screening in the training cohort, internal and external validation cohorts. Initially, univariate Cox proportional hazards regression was conducted to filter out variables demonstrating no significant correlation with OS in the training cohort (P > 0.05). Subsequently, the least absolute shrinkage and selection operator (LASSO) method was performed to select significant features that were strongly associated with patient overall survival.

Construction of clinical-molecular prediction model

Univariate Cox regression analysis was performed to analyze the correlation between clinical-molecular characteristics and OS, and those with P < 0.05 were included in the multivariate Cox regression analysis. Furthermore, independent risk predictors of OS were selected for the construction and validation of clinical-molecular models.

Radiomics signature construction and validation

These specific features were integrated into the radiomics signature, and the radiomics risk score was calculated based on the linear combination of weight coefficients assigned to each selected feature. Three sets of radiomics models (CE-T1WI, T2FLAIR, and combined-sequences) were established. Firstly, the correlation between the radiomics riskscore and OS was profiled by univariate Cox regression. The C-index of radiomics models was calculated, and the area under the curve was used to evaluate the diagnostic efficacy of the radiomics signatures in predicting 1-, 2-, and 3-year survival. The optimal radiomics signature was attained. The patients were divided into high-risk and low-risk groups according to a fixed cutoff value for Kaplan-Meier survival analysis. Finally, to demonstrate the incremental value of radiomics features, risk stratification was performed with the radiomics signature in various clinical and molecular subgroups. The detailed computational workflow is as follows:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=(\mathrm{c}\mathrm{o}\mathrm{e}{\mathrm{f}}_{\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_1}\times \mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_1)+\left(\mathrm{c}\mathrm{o}\mathrm{e}{\mathrm{f}}_{\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_2}\times \mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_2\right)+\cdots \left(\mathrm{c}\mathrm{o}\mathrm{e}{\mathrm{f}}_{\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_{\mathrm{n}}}\times \mathrm{f}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}{\mathrm{e}}_{\mathrm{n}}\right).$$

The coefficients were obtained from the multivariate Cox regression model of selected features.

Construction and evaluation of nomogram

A nomogram integrating the radiomics signature and clinical molecular risk factors independently predicted OS in glioma patients. The net reclassification index (NRI) was used to evaluate the enhancement of the radiomics signature to the clinical-molecular model. The calibration curve was used to assess the consistency between the predicted probability of survival and the actual survival. Decision curve analysis (DCA) compared the net benefits of the radiomics model, clinical-molecular model, and combined model at different threshold probabilities, which confirmed the clinical utility of the nomogram.

Statistical analysis

All statistical analyses were performed using R software, version 4.2.3 (http://www.R-project.org). The R packages used in this analysis were as follows: the “glmnet” package was used for LASSO Cox regression; Kaplan-Meier survival analysis and multivariate Cox regression were done using the “survival” package and the “survminer” package. The “timeROC” package was used to plot ROC curves; Nomogram and calibration curve were performed using the “rms” package; C-index calculation was implemented using the “Hmisc” package. A two-sided p-value of 0.05 indicated statistical significance.

Patient clinical characteristics

The clinical characteristics of all patients are shown in Table 1. For the TCGA/TCIA dataset, 414 patients were split into a training set (n = 331) and an internal validation set (n = 83) at an 8:2 ratio. The median OS of the training set and the internal validation set was 498 days (95% CI [443.49–552.51]) and 561 days (95% CI [402.29–719.71]), respectively. In the external validation set (n = 165), 103 patients died and 62 patients survived by the end of follow-up. The median follow-up duration was 1,245 days (95% CI [1,118.21–1,371.79]), with a median OS of 777 days (95% CI [628.85–925.15]). There were significant differences in the distribution of age, gender, pathological grade, IDH, and OS among the three groups of cohorts (P < 0.05), while only MGMT showed no statistically significant difference (P > 0.05).

Construction and validation of a clinical-molecular prediction model

According to univariate Cox regression analysis, age, pathological grade, IDH, and MGMT were significantly correlated with OS (P < 0.05), while gender was not statistically significant (P > 0.05). Multivariate Cox regression analysis further demonstrated that age, pathological grade, IDH, and MGMT remained as independent risk factors for glioma prognosis (P < 0.05). The specific P-values and HR are shown in Table 2. The C-index of the predictive model constructed based on clinical and molecular factors in the training set, the internal validation set, and the external validation set was 0.719, 0.713, and 0.689, respectively.

Radiomics signature construction and validation

As a whole, 19, 16, and 31 features were screened using LASSO Cox regression from CE-T1WI, T2FLAIR, and combined sequence images, respectively. The process of combined-sequences filtering features is shown in Fig. 3. In the training set and internal validation set, the combined sequences radiomics model demonstrated a higher area under the curve (AUC) (0.84 & 0.84, 0.84 & 0.91, 0.86 & 0.91) in predicting 1-, 2-, and 3-year survival compared to the CE-T1WI radiomics model (0.84 & 0.77, 0.85 & 0.80, 0.84 & 0.82) and the T2FLAIR radiomics model (0.77 & 0.65, 0.77 & 0.70, 0.76 & 0.73). Additionally, the multiparametric radiomics signature showed a better C-index compared to the single sequence model, making it the optimal radiomics model (Table 3). Hazard ratio (HR) values with 95% confidence interval (CI) and P values of each selected feature of the combined sequences radiomics model are shown in Fig. S1. The Spearman correlation heatmap illustrates the relationships between the selected significant radiomics features and clinical parameters in the training cohort (Fig. 4). Among the five clinical parameters, gender and MGMT status demonstrated stronger positive correlations with most radiomics features, while age, tumor grade, and IDH status exhibited weaker associations. The process for constructing the multiparametric radiomics signature was summarized in the supplementary documentation. The patients were split into high-risk and low-risk groups based on the median risk score (1.143) of the training set in three cohorts (Fig. 5).

Sub-stratified survival analysis of the radiomics signature

To prove the valuable predictive effectiveness of the radiomics model, subgroup survival analysis was performed to evaluate differences in OS between different clinical and molecular subgroups (Figs. 6 and 7). In the training set of clinical (age, gender) and molecular (MGMT) subgroups, there were statistically significant differences between high-risk and low-risk groups. The P-values of the log-rank test were all under 0.0001. In the internal validation set, significant differences were exhibited in each subgroup, with P-values less than 0.05. In the age, gender, IDH, and MGMT subgroups of the external validation set, the radiomics signature significantly stratified high-risk and low-risk groups according to the risk score cut-off value. While in the IDH mutation group of the training set and GBM group of the external validation set, risk stratification was not achieved.

Nomogram establishment and evaluation

The radiomics signature and selected risk factors generated a nomogram for OS prediction (Fig. 8a). Compared to the clinical-molecular model (C-index = 0.719 and 0.689) and radiomics signature (0.771 and 0.711), the nomogram showed superior ability to predict prognosis in the training (0.774) cohort and external validation cohort (0.776). The calibration curves showed delightful agreement between the predicted values for 1-, 2-, and 3-year survival probabilities and the actual values for the combined model in the training cohort, internal validation, and external validation cohort (Figs. 8B–8D). Upon further analysis, after integrating the radiomics signature into the clinical-molecular model, the integrated model produced an NRI of 0.545 (95% CI [0.115–1.058]), improving the proportion of correct classification. The comprehensive model provided more net benefit across most threshold ranges in contrast to other models in the decision curves analysis (Fig. 8E).