Introduction

The fifth edition of the World Health Organization (WHO) Classification of Tumors of the Central Nervous System (CNS) (WHO 2021), published in 2021, classifies gliomas into adult-type diffuse gliomas, pediatric-type diffuse low-grade gliomas, pediatric-type diffuse high-grade gliomas, circumscribed astrocytic gliomas, glioneuronal and neuronal tumors, and ependymal tumors based on histological and molecular pathology characteristics1. Among these, adult-type diffuse gliomas are the most common malignant tumors of the CNS2. The 2021 updated classification of diffuse gliomas depends largely on isocitrate dehydrogenase (IDH1/2) mutation status and 1p/19q co-deletion status, resulting in three primary disease groups: IDH-mutant, 1p/19q co-deleted oligodendroglioma (O); IDH-mutant, non-co-deleted astrocytoma (A); and IDH-wildtype glioblastoma (GBM). This new edition combines established histological diagnoses and molecular markers for the integrated classification of adult-type diffuse gliomas1,3. Survival rates vary widely depending on the glioma subtype, with low- and high-grade gliomas having 5-year survival rates as high as 80%, and below 5%, respectively4. Moreover, non-invasive methods are needed for patients who cannot safely undergo resection and for whom biopsy is performed only to obtain tissues for diagnosis. In addition, preoperative noninvasive prediction for the classification of adult diffuse gliomas is valuable for developing appropriate treatment plans.

Routine preoperative imaging for gliomas includes computed tomography (CT) and magnetic resonance imaging (MRI). Imaging is widely used in clinical practice to aid decision-making5, and advances in computational techniques have led to further developments in radiomics6. Texture analysis, which uses texture features to reveal the discriminant factors for tissue classification and image segmentation first appeared in MRI medical applications at the beginning of the eights7 and was successfully used for brain tumor classification and grading in the previous study8,9. In radiomics, medical images containing biological information are transformed into high-dimensional data (i.e., pre-established quantitative features) that can be further analyzed. This method has been utilized in clinical decision-making systems to improve diagnosis, treatment, and prognosis10,11. In general, deep learning features outperform radiomics features5,12,13. Segmentation of brain tumors using multimodal MRI plays a crucial role in biomedical and clinical analyses. Manual segmentation is the gold standard for tumor segmentation; however, it is costly, time-consuming, is limited by interobserver variability14.

Therefore, we used a deep learning method incorporating a U-net-based convolutional neural network (CNN)15 for automatic brain tumor segmentation directly from conventional MRI-acquired images, including the core and peritumor edema regions of the tumor and constructed a prediction model for the subtypes of glioma with machine learning techniques based on radiological features from the tumor core and edema, and explore its molecular pathological characteristics.

Materials and methods

Data source

1,251 preoperative multimodal MRI scans of gliomas for tumor segmentation task were obtained from organizers of the 2021 Brain Tumor Segmentation Challenge (BraTS2021)16. The Brain Tumor Segmentation Challenge (BraTS) aims to evaluate state-of-the-art methods for the semantic segmentation of brain tumors by providing a 3D MRI dataset with voxel-wise ground truth labels that are annotated by physicians17,18. In 2021, the organizers of the 2021 Brain Tumor Segmentation Challenge (BraTS2021) provided the largest training dataset. The training dataset includes 1251 subjects, each with four 3D MRI modalities: (a) native (T1) and (b) post-contrast T1-weighted (T1Gd, T1CE), (c) T2-weighted (T2), and (d) T2 Fluid-attenuated Inversion Recovery (T2-FLAIR), acquired with different protocols and various scanners from multiple institutions. A series of 184 patients including multimodal imaging and clinicopathological characterization data from the Zhongnan Hospital of Wuhan University was employed in this study to explore the correlation between glioma morphological features and pathological characteristics. Further, the same 184 cases of glioma obtained from the Zhongnan Hospital of Wuhan University and 240 cases of glioma obtained from the University of California San Francisco Preoperative Diffuse Glioma MRI dataset (UCSF-PDGM), known as The Cancer Imaging Archive (TCIA) were used for glioma subtype prediction task. In addition, 29 MRI scans of gliomas were included for testing the segmentation model and 184 cases of glioma were employed for exploring molecular pathological characteristics from Zhongnan Hospital of Wuhan University.

The study was conducted in accordance with the Declaration of Helsinki and approved by Medical Ethics Committee, Zhongnan Hospital of Wuhan University (No.2019048). Informed consent was obtained from all individual participants included in the study.

Glioma automatic segmentation

The network used in this study has an encoder-decoder structure, in which the encoder is down-sampled five times using convolution; similarly, the decoder follows the same structure using transpose convolution to up-sample five times, between the encoder-decoder. The encoder portion of the feature maps is then spliced to the corresponding decoder feature maps of the same level through jump links to provide richer feature information. The feature maps were normalized using instance normalization and LeakyReLU, an activation function, was used to compute the output. The T1-CE and T2W-FLAIR images were synthesized into a two-channel image as the input. Because the boundary between the tumor and edema is not obvious on T2W-FLAIR images, unlike conventional segmentation, the tumor and edema were outputted as two mutually exclusive labels, where the labels are preprocessed into two overlapping regions: namely, the tumor core region and the entire tumor region (containing both the tumor core and edema regions). Accordingly, the last layer of the network, the softmax layer, is replaced by the sigmoid activation layer, converting a multi-classification task of tumor region and edema region into two binary classification subtasks. A deep supervision mechanism was introduced to improve the gradient propagation computation process in the early layers. Initially, 32 convolutional filters were applied, up to a maximum of 320 convolutional filters in the top layer. The modified network structure is illustrated in Fig. 1.

Fig. 1
figure 1

Structure of the larger U-Net with asymmetric scaling on the encoder. The number at the top and bottom indicates the number of channels and the dimension of the features vectors, respectively.

Full size image

Standardized pre-processing has been applied to all these multimodal MRI scans, including alignment of patients’ modality images using rigid alignment, co-registration to the same anatomical template (SRI24)19, removal of background, skull, and scalp, isotropic resampling resolution interpolation to 1 mm. The pre-processing pipeline is publicly available through the Cancer Imaging Phenomics Toolkit (CaPTk) (https://cbica.github.io/CaPTk/) and Federated Tumor Segmentation (FeTS) tool (https://github.com/FETS-AI/Front-End/). After preprocessing, the image sizes of all the MRI scans and associated labels were 240 × 240 × 155 pixels. In addition, owing to the principle of MRI, which is influenced by the scanning equipment and individual patient characteristics, we normalized the images for each case by subtracting the mean grayscale of the brain region and then dividing by the standard deviation. We randomly divided the cases into training (n = 833 annotated images) and validation (n = 209 cases) sets. After that, the images were first randomly cropped to a patch size of 192 × 192 × 192 pixels to facilitate the up/down sampling operation in UNet and maintain a cubic like shape feature map, which the first two dimensions require cropping and the third dimension is padding to 192 with background value. A stochastic gradient descent (SGD) optimizer was used with the Nesterov momentum set to 0.99. The initial learning rate was set to 1e− 3 and decayed according to the following definition:

$$\:Ir=0.001\times\:{(1-\frac{epoch}{200})}^{0.9}$$

The loss function was defined as the sum of the binary cross-entropy and Dice losses. The T1-CE and T2W-FLAIR images were synthesized as two-channel inputs to train the defined U-Net network, and the pre-training model was obtained after 200 epochs. The training process was then conducted continually on the pre-trained model using clinical data. After 500 training epochs or when an early stop condition was satisfied (the value of the training loss did not decrease for at least 20 successive epochs), the final segmentation model was acquired. The training process was performed using PyTorch 1.12 on an NVIDIA RTX 3090 GPU (24 GB).

Glioma morphological characteristics analysis

In this study, we employed clinical data from Zhongnan Hospital of Wuhan University to analyze the correlation between pathological characteristics and morphological characteristics of glioma, including tumor volume, tumor surface, tumor longest diameter, tumor sphericity, edema volume and edema index which was calculated by the volume of the tumor core and edema by that of tumor core. Data of morphological characteristics were obtained using PyRadiomics toolkit based on segmentation mask. After tumor segmentation, the PyRadiomics toolkit can provide the calculation of parameters such as volume, surface area, sphericity, longest diameter, etc. Data of pathological characteristics were shown in Supplementary Material Table S1.

Prediction of subtypes of adult-type diffuse gliomas

424 cases of glioma, stratified sampling with three tumor subtypes, were then randomly divided into training (358 cases) and test (66 cases) data. The detailed distribution after data division is shown in Supplementary Material Table S2. Using the open source “PyRadiomics” package20, 96 features were extracted from the tumor core + edema 3D volumes of interest (VOI) obtained by segmentation from T2W-FLAIR and T1CE sequences, including First Order Statistics, Shape-based (3D), Gray Level Co-occurrence Matrix, Gray Level Run Length Matrix, Gray Level Size Zone Matrix, Neighbouring Gray Tone Difference Matrix, and Gray Level Dependence Matrix features in seven categories. Detailed information on the extracted radiomic features is provided in the Supplementary Material Table S3. Z-scores were used to normalize each feature. The median was used to fill for data with missing feature value. There are errors in the segmentation results, which may not correctly identify the tumor or edema area with poor MRI image quality. In these cases, some radiomic features are not computable, such as entropy, which requires logarithmic operations. We still keep the special cases to improve the robustness of the model and consider data of non-computable features were missing feature values. We applied five classifiers: namely, support vector machine (SVM), logistic regression, decision tree, random forest, and XGBoost. The hyperparameters of each classifier model were optimized using a random grid search and five-fold cross-validation. L2 Regularization penalty was used to avoid overfitting in SVM and Logistic regression by lightly smoothing the model weight. In other tree-based methods (decision trees, random forests, and XGBoost classifiers), the regularization parameters were used, such as “max_depth”, “min_samples_split”, and “min_samples_leaf”. For instance, by setting max_depth = 10, min_samples_split = 2, min_samples_leaf, = 1, the largest depth and the branches of a tree was controlled. These regularization parameters were also determined using a random grid search and five-fold cross-validation as stated previously.

Model evaluation

To evaluate the tumor detection results of the segmentation model, we used image data from 29 cases of glioma from the Zhongnan Hospital of Wuhan University for testing. An experienced neurooncology radiologist and a neurosurgeon used the 3D Slicer platform21 to manually segment the image in the test set for comparison. We used the following metrics: Dice score, Hausdorff distance, and volumetric similarity coefficient. The Dice score measures the overlap between the ground truth and the predicted segmentation (where 1 indicates perfect overlap). The Hausdorff distance measures the agreement between the boundaries of the ground truth and the predicted segmentation, in which a lower value indicates better agreement. Finally, the volumetric similarity coefficient measures the agreement in volume (where 1 indicates perfect agreement). The area under the curve (AUC), accuracy, sensitivity, and specificity were calculated to comprehensively evaluate each tumor subtype classifiers. The model with the highest average AUC was selected as the final classification model.

Statistical analyses

Mann–Whitney Test (95% confidence interval) was used to analyze differences in groups of two independent samples (IDH, 1p/19q, O6-methylguanine-DNA methyltransferase [MGMT], p53, epidermal growth factor receptor [EGFR]). The Kruskal–Wallis test was used for comparisons of groups with three or more independent samples (subtype, grade). Finally, the Pearson correlation coefficient was used for correlation analysis of continuous-discrete data (e.g., Ki67) Two-sided p < 0.05 indicated statistically significant differences. All data were analyzed using IBM SPSS Statistics for Windows, version 26.0 (IBM Corp., Armonk, New York, USA). Continuous variables are expressed as medians and interquartile ranges.

Results

Tumor segmentation algorithm performance

For automated segmentation, we achieved a mean Dice score of 0.884 for the tumor core and 0.889 for the whole tumor (tumor core + edema), a mean Hausdorff distance of 16.179 mm for tumor core and 4.34 mm for whole tumor, and a mean volumetric similarity coefficient of 0.918 for the tumor core and 0.9 for the intact tumor. Figure 2 shows the boxplots of these metrics for different patients in the test set. Figure 3 shows the automated segmentation performed using our method.

Fig. 2
figure 2

Dice scores, Hausdorff distances, and volumetric similarity for all patients in test set. (a) Dice scores for test set. (b) Hausdorff distances for test set. (c) Volumetric similarity for test set.

Full size image
Fig. 3
figure 3

Automatically segmented regions from segmentation model.

Full size image

Correlation between glioma morphological characteristics and subtypes

We observed morphological characteristics of glioma differed among clinical data, molecular levels, and pathological levels, respectively. Glioma tumor volume differed significantly across subtypes, with GBM having a significantly larger tumor volume than subtypes A (p = 0.022) and O (p = 0.002); however, subtypes A and O did not differ significantly (p = 1.000). The edema volume differed significantly across subtypes, with GBM having a significantly larger edema volume than subtypes A (p = 0.000) and O (p = 0.001); however, subtypes A and O did not differ significantly (p = 0.684). The edema index also differed significantly across subtypes, with GBM having a significantly larger edema index than A (p = 0.002) but no significant difference between A and O (p = 0.061) or between GBM and O (p = 1.000) (Fig. 4).

Fig. 4
figure 4

Relationships between glioma morphological characteristics and glioma subtypes. (a) Differences in tumor volume between glioma subtypes. (b) Differences in tumor surface between glioma subtypes. (c) Differences in tumor longest diameter between glioma subtypes. (d) Differences in tumor sphericity between glioma subtypes. (e) Differences in edema volume between glioma subtypes. (f) Differences in edema index between glioma subtypes.

Full size image

Glioma tumor volume, edema volume and edema index differed significantly among IDH mutation statuses (p = 0.000, 0.000, 0.009, Table S4, Fig. S1-1), with a significantly larger volume for IDH-wildtype than that for IDH-mutant, respectively. The tumor volume and edema volume both differed significantly among the 1p/19q states (p = 0.012, p = 0.009, Fig. S1-2), with a significantly larger tumor volume in 1p/19q non-co-deletion than the 1p/19q co-deletion. Tumor volume differed significantly across grades, with WHO grade IV having a significantly larger tumor volume than grade II (p = 0.000), but no significant differences between grades III and IV (p = 0.073) or grades II and III (p = 1.000). Edema volume also differed significantly between grades, with a significantly greater edema volume in grade IV than in grades II (p = 0.000) and III (p = 0.028) but no significant difference between grades III and II (p = 0.142). The edema index differed significantly across grades, with a significantly greater index in grade IV than in grade II (p = 0.005), but no significant difference between grades III and II (p = 0.253) or between grades III and IV (p = 1.000) (Fig. S1-3). The edema index differed significantly among p53 statuses, with a significantly greater index for the p53 wild type than for p53 mutant (p = 0.029, Fig. S1-4) while glioma showed no significant differences in any of the morphological characteristics according to EGFR amplification status (Fig. S1-5). The tumor surface area differed significantly across MGMT promoter methylation statuses, with a significantly larger area for the highest level of methylation compared with the lowest level (p = 0.023, Fig. S1-6). We also observed a significant association between Ki67 levels and morphological characteristics (tumor volume, tumor surface area, tumor diameter, and edema volume) in 177 patients (p = 0.012, 0.013, 0.016, 0.017, Table S5; Fig. 5). However, the correlations between Ki67 levels and these morphological characteristics are weak which we need to collect more pathological information data of patients for further analysis.

Fig. 5
figure 5

Relationships between glioma morphological characteristics and Ki67. (a) Pearson correlation coefficient between tumor volume and Ki67. (b) Pearson correlation coefficient between tumor surface and Ki67. (c) Pearson correlation coefficient between tumor longest diameter and Ki67. (d) Pearson correlation coefficient between tumor sphericity and Ki67. (e) Pearson correlation coefficient between edema volume and Ki67. (f) Pearson correlation coefficient between edema index and Ki67.

Full size image

Algorithm performance for the prediction of adult-type diffuse glioma subtypes

The AUC, accuracy, sensitivity and specificity were used as performance metrics for evaluating different models. Since the prediction of adult-type glioma subtypes was a multi-class prediction task, the macro-average AUC was finally used to compare performance between different classifiers. In the test set, random forest classifier achieved a highest AUC of 0.945 for the molecular subtypes of adult-type diffuse glioma prediction, which denotes it outperforms other classifiers on macro-average performance. The prediction of the “Glioblastoma, IDH-wildtype”, “Astrocytoma, IDH-mutant”, and “Oligodendroglioma, IDH-mutant, 1p/19q-coded” subtypes showed AUCs of 0.960, 0.914, and 0.961, respectively. The full results and corresponding receiver operating characteristic (ROC) curves are shown in Table 1; Fig. 6, respectively.

Fig. 6
figure 6

Receiver operating characteristic curves of the subtype prediction model based on random forest classifier in the test set.

Full size image
Table 1 Evaluation results of prediction model based on random forest classifier in the test set.
Full size table

Discussion

We developed a method for automatic brain tumor segmentation using 3D U-Nets integration. We further analyzed relevant influencing factors and clinical significance of the brain glioma tumor core area and peritumor edema and reported that edema volume of glioma differs among GBM and astrocytoma, GBM and oligodendroglioma, but not among astrocytoma and oligodendroglioma. These findings indicated that image features of edema play a role in glioma subtypes classification. Using the feature extraction parameters of the tumor and edema images obtained from this segmentation, we then developed a method to predict the molecular subtypes of adult-type diffuse gliomas based on the WHO 2021 classification. We constructed five classifiers that enabled the classification of glioma among extracted radiomics features and random forest showed the highest performance (AUC = 0.945 in the test set).

In the first phase of this study, we used BraTS21 public datasets, based on the nnUNet segmentation framework, and selected T1-CE and T2W-FLAIR sequences according to the actual imaging protocols of the clinical dataset for segmentation of the tumor core and whole tumor (tumor core and peri-tumor edema). Finally, utilizing pre-training of a public dataset and micro-tuning of a clinical dataset using only two imaging sequences, we achieved a performance close to those of two programs22,23 that performed well with BraTS21.The consensus recommendations for a standardized Brain Tumor Imaging Protocol24,25 suggest acquiring T1W, FLAIR, DWI, T2, and T1CE sequences. However, in clinical radiology practice, it is not always possible to acquire all imaging sequences of sufficiently complete and high quality for various reasons such as contrast intolerance and motion artifacts caused by unavoidable patient movement. Therefore, the method proposed in this study, which is based on fewer sequences while guaranteeing high accuracy in tumor segmentation, has a high value in clinical practice. Polina et al.26 evaluated the contributions of individual sequences to multimodal tumor segmentation, and reported that just two imaging sequences, T1CE and T2W-FLAIR, achieved comparable performance as the four full imaging sequences under the same experimental conditions; however, their study used only conventional 3D U-net. Without adjusting the network structure and training parameters, the performance was not outstanding. Raphael et al.27 investigated the optimal order of tumor imaging sequences incrementally, showing that the T2W-FLAIR sequence improved the segmentation of peritumoral edema. Based on the above studies, the results of the present study demonstrated that the segmentation model based on the nnUNet framework can be used for good segmentation performance of the tumor core and peritumor edema using a limited number of sequences (T1-CE + T2W-FLAIR).

WHO 2021 expands upon the trend that started in 2016, using key molecular biomarkers to define neoplastic entities and greatly reducing the dependency on morphological features for tumor classification2. Due to the close relationship between IDH mutations, 1p/19q co-deletion, and patient prognosis28, previous studies have predicted IDH status and 1p/19q co-deletion using multimodal MRI radiomics features29, residual convolutional neural networks30,31, and 3D-Dense-UNet32. Furthermore, a previous study showed that the molecular subtypes of diffuse gliomas can be predicted comprehensively using radiomics analysis33. If two separate methods are used, one predicting IDH status and one predicting 1p/19q co-deletion status, an IDH wild-type glioma can be predicted to have a 1p/19q co-deletion; however, as IDH wild-type and 1p/19q co-deletion are mutually exclusive34, this is irrelevant to clinical reality.

In this study, we predicted the subtypes of adult-type diffuse gliomas using a trained model based on the latest classification of adult-type diffuse gliomas by the World Health Organization in 2021, using a well-trained segmentation network that extracts features from the tumor core and edema 3D VOI obtained by automated segmentation from T2W-FLAIR and T1CE sequences as inputs. The prediction performance achieved an AUC of 0.945. The prediction of each of the three subtypes was also highly accurate (GBM: 90.9%, A: 86.3%, and O: 92.4%; Table 1). The network identified “GBM” (100%) with higher sensitivity than “A” (57.1%) and “O” (33.3%). This result was not surprising because of the remarkable intra-tumor heterogeneity of IDH wild-type GBM with areas of necrosis, blood-brain barrier breakdown, and extensive perifocal edema2, which likely aided the network in analyzing and identifying relevant MRI features to correctly identify IDH wild-type GBM with high sensitivity. Cluceru et al.31 and Golestan Karami et al.35 also performed a three-group classification analysis and observed higher sensitivity in the IDH-wildtype group compared with IDH-mutant astrocytoma and IDH-mutant 1p/19q-co-deleted oligodendroglioma.

Wang et al.3 proposed a CNN-based comprehensive diagnostic model for the integrated classification of adult-type diffuse gliomas using deep learning from pathology images. However, pathology images can only be acquired after surgery or biopsy and cannot provide preoperative guidance, especially in cases that are difficult to diagnose preoperatively or when an operation is not performed. In contrast, our method provides a non-invasive alternative. Karami et al.35 combined multi-shell diffusion with conventional MRI to apply deep learning for the molecular diagnosis of diffuse gliomas; however, their accuracy for the three molecular subtypes according to WHO 2021 was low (60 ± 5%). In addition, their study used an independent dataset from only one hospital. In contrast, our study utilized data from 424 cases in two datasets, the Zhongnan Hospital of Wuhan University and the UCSF-PDGM (TCIA) public database, to prevent an overestimation of the performance of our method. We used Z-scores to normalize each radiomic feature and we used the median for data with missing feature values considering the multi-center data source. Z-score is a commonly used normalization method which is suitable for outlier data beyond the value range, thus greatly improving the stability of the model. This also allowed our method to be robust to the heterogeneity naturally present in clinical imaging data, enabling its broad application in clinical practice. Tang et al.36 used a multitask network to predict multiple tumor genotype-related features and overall survival in patients with glioblastoma. Because their method only works for patients with glioblastoma, the tumor grade must be known in advance, which complicates its use in preoperative situations when tumor grading is not yet known. Moreover, studies on the molecular subtypes of gliomas5,36,37 generally require tumor segmentation as input, which is a time-consuming task. The automatic segmentation method for glioma tumors used in this study can be applied to prediction models.

Our tumor segmentation method has room for expansion. Tumor segmentation was used for tumor subtype prediction in this study. We did not study the effects of different sequences or different combinations of sequences on the segmentation tasks and results. Moreover, we applied the basic nnUNet framework and did not investigate the mechanism of attention, group normalization, or other potential methods for optimization. In addition, a larger number of samples and more clinical data can also enable further improvement in the performance and applicability of the methods developed in the present study.

The correlation analysis results between pathological features and morphological characteristics are unstable due to the factors of outliers. The reason for this situation may be the imbalance of data or we probably need more detailed glioma subtype classification to analysis the differences of each other. Moreover, although our method showed good overall performance, the performance differed between tumor categories. The sensitivity for predicting oligodendrogliomas was low, which could be attributed to the lack of a central pathology review. This difference in the predictive performance between the subgroups can also be attributed to an imbalance in the data which possible reason is that cases of the other two subtypes are rare. Therefore, although our method may be relevant for certain subgroups, further improvement is necessary to ensure its relevance to the entire patient population and improve the performance and robustness of model, especially for these complicated cases.

In future studies, we aim to optimize our method by including perfusion-weighted imaging (PWI) and diffusion-weighted imaging (DWI) as these imaging modalities may contain additional information that can be linked to tumor genetic traits and aggressiveness38.

In conclusion, we developed a prediction model for molecular glioma subtype in brain glioma using parameters obtained from the automated segmentation method we also developed in this study. We further identified influencing factors and confirmed the clinical significance of the tumor core area and peritumoral edema in brain gliomas. Additional studies are needed to further train and enhance our model and demonstrate its potential in the clinical setting.