2.1

Patient cohorts

The current study utilizes data from the multimodal brain tumor segmentation (BraTS) challenge[17, 18] and the other part of the dataset is from Sun Yat-sen University Cancer Center (SYSUCC). These sources provide comprehensive imaging data, essential for advancing the accuracy and reliability of glioma segmentation methodologies. We randomly split the BraTS into training and internal test data, and the SYSUCC dataset employed as an external testing set. The MR data of glioma patients admitted from November 2018 to March 2022 to SYSUCC. Ethical approval of the SYSUCC dataset was obtained for this retrospective study and informed consent from patients was waived. Figure 1 shows the specific details of patient enrolment pathway.

Figure 1.

Enrolment pathway for patients in two centers

2.2

Preprocessing and segmentation

To minimize heterogeneity, all MR images were standardized by reorienting to the left-posterior superior coordinate system, co-registering to the T1 anatomical template from the MNI152 atlas, resampling to a 1mm isotropic resolution, and applying skull stripping using the brain extraction tool [18, 19]. Each imaging dataset was manually segmented by one to four raters following a consistent annotation protocol, with all annotations validated by experienced neuroradiologists. The delineation process for each ROI in the SYSUCC dataset involved collaboration between two readers with 3 and 4 years of experience in brain tumor diagnosis, respectively, and was reviewed by a senior reader with 11 years of experience. Based on the BraTS challenge guidelines, three types of ROIs were considered: non-enhancing tumor, enhancing tumor, and edema. This study focuses on the combined analysis of these glioma regions.

2.3

Feature fusion

This study considered three categories of radiomics features, encompassing shape features, first-order statistics features, and texture features. The texture features included features derived from the grey level co-occurrence matrix (GLCM), grey level run length matrix (GLRLM), grey level size zone matrix (GLSZM), grey level dependence matrix (GLDM), and neighbor grey tone difference matrix (NGTDM). Prior to feature extraction, seven types of image filters were applied, including the origin filter (no filter applied), Laplacian of Gaussian (LoG) filter, Square filter, square root filter, exponential filter, gradient filter, and wavelet filter. In the case of LoG filtering, a parameter sigma, which emphasizes either fine or coarse textures, was set to 3, 4, or 5. Wavelet filtering involved 8 decompositions per level, considering all possible combinations of applying either a high or a low pass filter in each of the three dimensions. The grey level within the region of interest was quantized separately with a specific bin count of 32. A total of 1470 features were extracted for each MR sequence, resulting in 5880 features per patient.

Before utilizing the extracted features, batch effect correction was applied to alleviate batch effects from different center datasets[20]. Shapley values, derived from coalitional game theory, address how groups of individuals can collaboratively work towards a common goal [21]. This explanatory method integrates optimal credit allocation with local explanations by calculating Shapley values. In this context, each feature value of the subject is treated as a player in a cooperative game, where the prediction represents the total payout. The significance of a feature, denoted as z, is quantified using the Shapley value, defined as follows:

where F is the complete set of all features, A represents the subsets obtained from F without the feature z. f(A) is the output of a specific model to be explained using A. ψ_z is the Shapley value or the contribution of feature z. The proposed feature Shapley analysis (FSA) aims to estimate modality weights to fuse features from multiple classifiers based on training prediction performance and Shapley values. Specifically, let us define N classifiers C_i. The predictive performance of each classifier during the training phase can be quantified using evaluation metrics, including accuracy, area under the receiver operating characteristic curve (AUC), F1 score, specificity, and sensitivity. These metrics collectively form an evaluation matrix E_i,j(i = 1, 2, …, N, j = 1, 2, …, K) that provides a comprehensive assessment of the classifier’s effectiveness. It is organized as a matrix with rows corresponding to N classifiers and columns representing K evaluation criteria. The modality importance weight of each classifier can be described as S_i,k(i = 1, 2, …, N, k = 1, 2, …, M), where M represents the number of modalities. The FSA estimates the fusion weights ω_k of each modality by using

where the c_i and S_i,k are firstly normalized as

In (3), E_max,j and E_min,j are defined as

Here, E_max,j and E_min,j are the Euclidean distance between the vectors , for each classifier i, where e_max and e_min represent the maximum and minimum values for each criterion, respectively, and are defined as follows:

In (4), S_i,k is a matrix of N rows of ψ_k vectors with rows and columns representing N classifiers and M modality features with Shapley value. With the normalized weights ω_k, the final fusion feature can be calculated as:

There were totally 1470 features for each patient after fusion.

2.4

Model construction

We implemented a coarse-to-fine feature selection strategy to reduce the dimensionality of the feature space to mitigate the risk of overfitting. Initially, a univariate logistic regression (LR) analysis was conducted to assess the correlation between individual features and glioma grades. Features with statistically significant correlations (p < 0.05) were retained for further analysis. Next, we employed a minimum redundancy and maximum relevance feature selection framework to diminish mutual redundancy among features while preserving those with maximum relevance [22]. In the final stage of feature selection, we utilized an embedding method, specifically the multivariate LR with L1-regularization, to identify the most informative features.

To identify the most effective prediction model, we employed four machine-learning algorithms for constructing combined models: LR, random forest (RF), support vector machine (SVM), and gradient boosting tree (GBT). It’s worth noting that our dataset exhibited a degree of class imbalance, with the size of the minority class approximately one-fourth that of the majority class. To address the class imbalance, we applied the synthetic minority oversampling technique [23] for data augmentation, ensuring an equal ratio of HGG and LGG. During the training phase, a 10-fold cross-validation approach was implemented to determine hyperparameters for the models.

The predictive model performance was evaluated using different metrics such as the receiver operating characteristic (ROC) curve, AUC accuracy, F1 score, sensitivity, and specificity. Both internal and external test cohorts were used to validate the generalizability of the predictive model. All above were conducted using the Python programming language (version 3.6.13). The specific flowchart is shown in Figure 2.

Figure 2.

Overview of the modality fusion analysis framework for the grade of gliomas.

3.1

Patient characteristics

We utilized the available section of the BraTS 2020 dataset, reserving 70% for the training set. In the development of our glioma grading model, cross-validation was implemented within a training set comprising 257 patients with gliomas (204 HGGs and 53 LGGs). The remaining 30% of the dataset was designated for internal model testing. For the independent external testing dataset, we retrospectively gathered 136 patients who underwent MRI examinations and received pathological confirmation of gliomas. It’s essential to note that each patient contributed four MRI modalities, including T1, T2, FLAIR, and T1C. Consequently, a total of 2016 MR images from 504 glioma patients were included in this study. Table 1 provides a summary of the datasets employed.

Table 1.

Summary of the dataset used in this study.

3.2

Model evaluation

In Figure 3, ROC curves illustrate the performance of all models in the training cohort, internal test cohort, and external test cohort. LR exhibits a slightly inferior overall performance. RF demonstrates strong performance across all datasets, achieving an AUC of 0.864 on the external validation set. SVM attains the highest AUC on the internal validation set (0.951). In addition, we carried out integrated learning of above models and get the relatively best model. The integration method we used is the voting strategy and the voting ensemble model (VEM) achieved the best AUCs of 0.954 on the internal validation set and 0.856 on the external validation set.

Figure 3.

Receiver operating characteristic (ROC) curves of prediction performance of modality SHAP model.

3.3

Fusion interpretability

To visualize the features fusion process, we generated a cumulative histogram illustrating the fusion of modal features, as depicted in Figure 4. This plot conveys the cumulative results of SHAP values assigned to features selected from multiple modalities before normalization. The visualization highlights the varying contributions of different modalities to each feature; certain modalities exhibit a significant impact on specific features, while others have a more modest influence.

Figure 4.

Histogram and cumulative histogram of the modal weights of contribution among radiomic features.

To quantitatively elucidate the ensemble model’s capability in distinguishing glioma grades, we employed the SHAP method to assess feature importance and their contributions to the grading process. Figure 5 illustrates the SHAP summary plot of features post-feature selection, highlighting their relative significance and impact on the model. The y-axis represents multimodal features ranked by importance, while the x-axis depicts SHAP values. Each point on the plot corresponds to a sample’s SHAP value, with color coding indicating the magnitude of each feature. The analysis reveals that features derived from GLDM, first-order statistics, GLRLM, GLCM, and GLSZM—extracted from filtered and transformed MRI images—are critical in differentiating glioma grades. Additionally, certain features show pronounced correlations between LGG and HGG patients, with higher feature values often associated with a greater likelihood of HGG. For example, a lower value of wavelet-HHL_gldm_LargeDependenceHighGrayLevelEmphasis, synthesized from four modalities, correlates with a higher SHAP value, suggesting an increased probability of HGG development.

Figure 5.

SHAP summary plots of the top 15 features.

To gain a deeper understanding of how the interaction between features influences the model output, SHAP dependence contribution plots can unveil these interactions’ impact. We chose the most critical feature (square_glrlm_ShortRunEmphasis) to elucidate this impact. As illustrated in Fig. 6 (a), the red points represent patients with higher values of interactive features, while the blue points represent those with lower values. Fig. 6 (a) demonstrates that the SHAP value for higher square_glszm_ZoneEntropy is elevated when square_glrlm_ShortRunEmphasis is less than zero, indicating a lower risk of HGG. Conversely, when square_glrlm_ShortRunEmphasis surpasses one, patients with smaller square_glszm_ZoneEntropy may be developing HGG. In Fig. 6 (b), it is evident that the SHAP value for square_ngtdm_Busyness increases as square_glrlm_ShortRunEmphasis decreases, signifying a decreased likelihood of developing HGG.

Figure 6.

Scatter plots of fusion features SHAP dependence analysis. (a) is impact of fused square_glszm_ZoneEntropy and square_glrlm_ShortRunEmphasis on final model output, and (b) is impact of fused square_ngtdm_Busyness and square_glrlm_ShortRunEmphasis on final model output.

3.4

Comparison with different fusion methods

We compare the proposed FSA with the following fusion methods. First, the current mainstream methods of features concatenation are compared. Then, simple and effective fusion methods, maximum and average value of features which always used in pooling of deep learning [24, 25] were also compared. Additionally, discriminant correlation analysis (DCA) [26] associates class relationships with the correlation analysis of feature sets, effectively implementing feature fusion. DCA maximizes the pairwise correlation between two feature sets while eliminating inter-class correlation and restricting correlation within each class. In addition, two deep learning fusion methods were implemented to compare with our proposed fusion method. Cheng et al proposed e a glioma grading model termed as multimodal disentangled variational autoencoder (MDVAE) [27] by integrating the complementary information from multimodal MRI images. Liang et al. proposed an association-based fusion method (AF) [28] for multi-modal classification, aimed at integrating complementary information from various modalities to enhance classification performance. This method leverages the strengths of each modality, thereby improving the overall accuracy and robustness of the classification system. After implementing the four methods mentioned above, following the same feature selection and model training procedures, we obtained the comparative results for each model as shown in Table 2.

Table 2.

Comparison method result.