Ethics approval for this study was obtained from the Human Scientific Ethics Committee of the First Affiliated Hospital of Zhengzhou University (No. 2019-KY-176). Among the 604 patients receiving craniotomy for tumor resection and pathologically diagnosed with LGG in the Department of Neurosurgery, the First Affiliated Hospital of Zhengzhou University, from July 2009 to July 2019, 335 were further selected according to the following criteria: (1) adult patients (age ≥18 years); (2) histopathological diagnosis of primary grade II/III glioma; (3) availability of IDH and 1p/19q status; (4) availability of preoperative multiparametric MRI, including axial T1WI, CE-T1WI, T2WI, FLAIR, and ADC imaging; and (5) availability of sufficient image quality without significant artifacts, as determined by neuroradiologists and neurosurgeons. The selection procedure is shown in Supplementary Figure S1 . Clinical factors (gender and age) were obtained from the medical record system.

All patients were examined on either 1.5- or 3.0-T clinical MR scanners from Siemens, Philips, or GE Healthcare. The brain MRI protocol included the following sequences: (a) axial and sagittal T1WI, (b) axial T2WI, (c) axial FLAIR imaging, and (d) DWI and the corresponding ADC maps generated with the software incorporated into the MRI unit; and (e) axial, sagittal, and coronal CE-T1WI obtained immediately after intravenous administration of a gadolinium-based contrast agent. Details of the MR machines and sequence parameters are provided in the Supplementary Material .

The workflow of this study is shown in Figure 1 . Image preprocessing was performed to standardize the images. First, the N4 bias field correction was applied to remove any low frequency intensity non-uniformity (20). After all voxels were isotropically resampled into 1 × 1 × 1 mm³ using trilinear interpolation, multiparametric MRI samples were co-registered to the corresponding CE-T1WI using a rigid transformation. Histogram matching was used to normalize signal intensity. A batch-effect correction tool ComBat (21) was used to remove scanner and site effects. The three-dimensional volume of interest (VOI) of tumor contours was manually delineated using the open-source software ITK-SNAP (www.itk-snap.org) by a neuroradiologist (JY with 11 years of experience), primarily on axial FLAIR images where T2WI and CE-T1WI were used to cross-check the extension of the tumor and fine-tune the tumor contour. According to BraTS subvolume convention (22), the VOIs were delineated as the whole tumor region, including the edema, enhancing core, the non-enhancing core, and the necrotic/cystic core. The VOIs were re-delineated in 60 cases (IDHwt, n = 20; IDHmut-noncodel, n = 20; IDHmut-codel, n = 20), which were randomly selected among the included patients by a senior neuroradiologist (JC with 20 years of experience). The segmented VOI was then overlaid with the co-registered resampled T1WI, CE-T1WI, T2WI, FLAIR, and ADC images. Neuroradiologists were blinded to clinical, pathological, and molecular data.

All radiomics features were extracted using Pyradiomics extractor. Three groups of features were extracted: shape features, first-order intensity features, and higher-order texture features. Texture features were extracted using five different methods, including the gray-level co-occurrence matrix (GLCM), gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), gray-level dependence matrix (GLDM), and neighborhood gray-tone difference matrix (NGTDM). Two filters [wavelet transform and Laplacian of Gaussian (LoG) with four sigma levels (2.0, 3.0, 4.0, and 5.0)] were enabled in the extracted intensity and texture features. The extracted features are summarized in Supplementary Table S1 . The extracted features were consistent with the Imaging Biomarker Standardization Initiative (IBSI) (23).

We selected 25 semantic descriptors of imaging features from VASARI annotations based on preoperative MRI: f1, tumor location; f2, side of lesion center; f3, eloquent brain; f4, enhancement quality; f5, proportion enhancing; f6, proportion non-contrast-enhancing tumor; f7, proportion necrosis; f8, cysts; f9, multifocal or multicentric; f10, T1/FLAIR ratio; f11, thickness of enhancing margin; f12, definition of the enhancing margin; f13, definition of the non-enhancing margin; f14, proportion of edema; f15, edema crosses midline; f16, hemorrhage; f17, diffusion characteristics; f18, pial invasion; f19, ependymal invasion; f20, cortical involvement; f21, deep white matter invasion; f22, non-enhancing tumor crosses midline; f23, enhancing tumor crosses midline; f24, satellites; and f25, calvarial remodeling. The exact description of all these features can be found in The Cancer Imaging Archive VASARI research project webpage (https://wiki.cancerimagingarchive.net/display/Public/VASARI+Research+Project). Each tumor was independently scored by a neuroradiologist (JY) and a neurosurgeon (ZZ with 11 years of experience), according to the VASARI scoring system based on the five MR sequences using ITK-SNAP. Furthermore, the T2-FLAIR mismatch sign, which has been defined as an easily detectable imaging sign on routine clinical MRI studies for the diagnosis of IDHmut-noncodel gliomas (24), was also assessed. Any disagreement between the two raters was resolved through discussion and consensus. A qualitative feature dataset was obtained by combining the VASARI features and T2-FLAIR mismatch signs.

For radiomics features, all features were standardized using z-score normalization. First, the stability of the extracted features was evaluated by interobserver reproducibility of the two image readers. Intraclass correlation coefficient (ICC) values were calculated for each feature of the 60 patients. Features with ICC value ≥0.90 were selected in this study. The correlation coefficient between each pair of features was calculated to eliminate redundancy. For feature pairs with correlation coefficients >0.75, the feature with the worst univariate predictive power (smaller Mann–Whitney U-test p-value) was removed. Based on the qualitative feature dataset and the remaining robust and non-redundant radiomics feature data set, the R package Boruta (25) was used to select the optimal all relevant features. Boruta is a random forest-based all-relevant feature selection wrapper algorithm that iteratively compares the importance of original features with the importance of artificially added random features, progressively removing irrelevant features. The most important features of the qualitative and radiomics feature datasets were obtained.

First, a radiomics model based on selected radiomics features was constructed using the R package randomForest to classify the three molecular subtypes. Then, a combined model with selected radiomics features, selected qualitative features, and clinical factors was constructed using the random forest algorithm. For comparison, a clinical model based on clinical factors (gender and age) and a qualitative model based on VASARI features and T2-FLAIR mismatch signs were constructed using the same algorithm. Besides, a radiomics model without ADC sequence was also built. Gini index was used as importance measure (26).

All statistical analyses were performed using the R software (version 4.0.5, http://www.Rproject.org). Statistical significance was set at p < 0.05. The patients in this study were randomly divided into a training cohort and a testing cohort at a ratio of 4:1, where the distribution of the clinical characteristics was balanced. Differences in gender, age, and molecular subgroups between the training and testing cohorts were assessed using the t-test or χ² test. Differences in patient characteristics across the three molecular subtypes were assessed using the Kruskal–Wallis test. In the training cohort, 10-fold cross-validation was applied to optimize the parameters of random forest classifiers in all four classification models (combined model, radiomics model, qualitative model, and clinical model). The testing cohort was used for the final model evaluation. The classification performance (one specific class versus all other classes) was assessed using receiver operating characteristic (ROC) analysis according to the area under the curve (AUC), balanced accuracy, sensitivity, and specificity. The maximum value of the Youden index (sensitivity + specificity − 1) was chosen as the optimal cutoff for each binary classification. All indices were calculated for both training and testing cohorts. The AUCs were statistically compared between different classifiers using DeLong analysis.

A total of 335 patients were included in the current study according to the selection criteria. The patients were divided into training (n = 269) and testing (n = 66) cohorts. There were no significant differences in clinical factors and molecular subtypes between the training and testing cohorts, as shown in Table 1 . The distribution of patient characteristics across the three molecular subtypes is shown in Supplementary Table S2 .

We extracted 1,197 features from each sequence ([14 shape features, 234 intensity features (18 were from original images, 72 were from LoG images, and 144 were from wavelet images), 949 texture features (73 original texture features, 292 LoG texture features, and 584 wavelet features)]. In total, 5,929 radiomic features were extracted from the five MRI sequences for each patient. After the robustness tests, 3,103 out of 5,929 features remained. After redundancy reduction, 335 features were selected for the subsequent analyses. The heat maps of the correlation coefficients of both the 3,103 features and the selected 335 features are shown in Supplementary Figure S2 . After the Boruta feature selection, the 17 most important features for an optimal model fit were finally selected, including 11 texture features and 6 intensity features, as shown in Table 2 . The results of the Boruta feature selection are shown in Supplementary Figure S3 , where the boxplots of the importance of all features fed to Boruta are shown. The univariate association of each selected feature with the molecular subtype was significant (false discovery rate adjusted, p < 0.001). There were no significant differences in selected 17 radiomics features between different scanners and different field strengths (ANOVA test, p > 0.05), as shown in Supplementary Figures S4 and S5 , respectively.

The Boruta algorithm revealed 10 qualitative features that were significantly associated with the three molecular subtypes. The results of the Boruta feature selection are shown in Supplementary Figure S6 . Finally, nine VASARI features and T2-FLAIR mismatch sign were selected to build a qualitative model: f1, tumor location (V1); f4, enhancement quality (V2); f5, proportion enhancing (V3); f6, proportion non-contrast enhancing tumor (V4); f11, thickness of enhancing margin (V5); f13, definition of the non-enhancing margin (V6); f20, cortical involvement (V7); f22, non-enhancing tumor crosses midline (V8); f23, enhancing tumor crosses midline (V9); and T2-FLAIR mismatch (V10).

The ROC curves of the combined model, radiomics model, clinical model, and qualitative model for the training and testing cohorts are shown in Figure 2 . The AUC of the radiomics model was 0.6557 for IDHwt, 0.6830 for IDHmut-noncodel, and 0.72579 for IDHmut-codel in the testing cohort. When combining the radiomic features with the qualitative features and clinical factors, the AUCs of the combined model were 0.8623 for IDHwt, 0.8056 for IDHmut-noncodel, and 0.8036 for IDHmut-codel, in the testing cohort, with balanced accuracies of 0.8924, 0.8066, and 0.8095, respectively. Significant differences in AUCs between the radiomics model and the combined model were found for IDHmut-noncodel and IDHmut-codel in the training and testing cohorts (DeLong p < 0.05). The AUCs of the clinical model in the testing cohort were 0.6551 for IDHwt, 0.4841 for IDHmut-noncodel, and 0.5873 for IDHmut-codel. In the qualitative model, the AUCs in the testing cohort were 0.7488 for IDHwt, 0.7599 for IDHmut-noncodel, and 0.7892 for IDHmut-codel. The AUCs of the radiomics model without ADC sequence in the testing cohort were 0.7234 for IDHwt, 0.5144 for IDHmut-noncodel, and 0.6533 for IDHmut-codel, which have significant differences compared with the radiomics model for IDHmut-noncodel and IDHmut-codel (DeLong p < 0.05). The classification performance of the radiomics model and the combined model in both the training and testing cohorts are summarized in Tables 3 and 4 , respectively. The performance of the clinical, qualitative, and no ADC sequence radiomics models is shown in Supplementary Tables S3–S5 , respectively.

To describe the univariate contribution of each parameter used for subtype classification, a heat map of the subtype-specific parameter importance in the classification is shown in Figure 3 . The meanings of the 17 radiomics features are detailed in Supplementary Table S6 . The Gini index was calculated as the importance value for building the combined model, indicating the univariate contribution to the classification. A larger value indicates greater importance in classifying a specific subgroup.