Multimodal data integration using deep learning predicts overall survival of patients with glioma

Abstract

Gliomas are highly heterogenous diseases with poor prognosis. Precise survival prediction could benefit further clinical decision‐making, clinical trial incursion, and health economics. Recent research has emphasized the prognostic value of magnetic resonance imaging, pathological specimens, and circulating biomarkers. However, the integrative potential and efficacy of these modalities require to be further validated. After incorporating 218 patients of The Cancer Genome Atlas glioma datasets of and 54 patients of the Huashan cohort with complementary prognostic information, we established a squeeze‐and‐excitation deep learning feature extractor for T1‐contrast enhanced images and histological slides and explored to screen significant circulating 5‐hydroxymethylcytosines (5hmC) profiles for glioma survival by least absolute shrinkage and selection operator‐Cox regression. Therefore, a prognostication predictive model with high efficiency was obtained through survival support vector machine multimodal integration of radiologic imaging, histopathologic imaging features, genome‐wide 5hmC in circulating cell‐free DNA and clinical variables, suggesting a valid strategy (concordance‐index = 0.897; Brier score = 0.118) for improved survival risk stratification of glioma patients.

INTRODUCTION

According to the Global Cancer Statistics 2020, an estimated 308,102 new cases and 251,329 deaths of brain tumors were reported worldwide. Gliomas are reported as the most prevalent primary central nervous system cancers, accounting for approximately 80% of cases.^1,2 The survival for patients with glioma has improved minimally over the past 20 years. The median overall survival (OS) was 78.1, 37.6, and 14.4 months for low-grade gliomas (LGGs), anaplastic gliomas, and glioblastomas (GBMs), respectively.³ The prognosis of gliomas with multiple grades and subtypes varies immensely. An accurate assessment of individual OS is vital not only for enhancing the survival of patients but also for guiding subsequent comprehensive treatment protocols. The early, objective, and precise estimation empowered clinicians to opt for more aggressive surgical interventions and radiotherapy combinations, prolonging the survival duration of patients. Recent studies have concentrated on the prognostic prediction of gliomas using diverse modalities and algorithms, including radiologic imaging, histopathologic results, and proteogenomic data.⁴ However, current predictive models have several nonnegligible technical bottlenecks, impeding their wide applications.

First, these models specifically applied a single modality, neglected the broader clinical context, and restrained validation, leading to efficiency overestimation due to overfitting. Generally, the patient state is characterized by a whole spectrum of modalities, ranging from radiology, histology, and genomics to electronic health records, providing a potential opportunity of multimodal integration for diverse tumors, such as ovarian cancer,⁵ colon adenocarcinoma,⁶ and medulloblastoma.⁷ Besides clinicogenomic features, multiscale clinical imaging is routinely performed during the course of care, including contrast-enhanced magnetic resonance (MR) and hematoxylin and eosin (H&E)-stained slides. This multiscale medical imaging modality and slides provide cancer features and spatial context information of tumors; however, establishing image feature extractors and algorithms is indispensable for further prediction.

Second, most existing models or signatures have limited model robustness and transportability to accommodate independent populations.⁸ Efficient biomarkers for patient stratification, definition of risk groups, and monitoring of response to therapy have become an integral component of clinical trials in oncology. Novel and more specific biomarkers need to be identified. In the human genome, 5-hydroxymethylcytosines (5hmC) are abundant epigenetic features generated by the oxidation of 5-methylcytosines through the ten–eleven translocation enzymes. Recent studies have suggested that 5hmC modifications are related to cancer pathobiology, including the observed global reduction of 5hmC levels in various cancer types.^9–11 The 5hmC modifications in promoters, gene bodies, and gene regulatory elements (e.g., enhancers) faithfully reflect gene expression activation in mammalian genomes, and thus can serve as ideal markers for specific gene/locus activation in chromatin. Molecular signatures in circulating cell-free DNA (cfDNA) based on cytosine 5-hydroxymethylation have been demonstrated previously to potentially define the tissue of tumor origin in various disease types.^9–13 However, only a few studies focused on their prognostic role in cancers,^14,15 which is required for comprehensive validation in further studies.

The present study aimed to evaluate the prognostic potential of multimodal integration for radiologic imaging, histopathologic imaging, circulating 5hmC profiles, and clinical characteristics. We obtained a Cox proportional hazards (CoxPH)–multimodal predictor (MMP) and concise nomograms to stratify patients with glioma by survival for further clinical decision-making.

METHODS AND MATERIALS

Study design, data acquisition, and patient selection

The glioma datasets from The Cancer Genome Atlas (TCGA) program were included and applied in the present study. The TCGA glioma cohort consisted of TCGA-GBM and TCGA-LGG cohorts. The inclusion criteria for patient selection were as follows: (1) patients with pathologically diagnosed gliomas; (2) those with preoperative MR imaging (MRI) data with T1-weighted gadolinium contrast-enhanced (T1-CE) imaging; and (3) those with available survival data and clinical characteristics, including sex, age, World Health Organization (WHO) grade, tumor location (motor and language functional area, MLFA), and isocitrate dehydrogenase 1 (IDH1) mutation status. The MRI scans of individuals were obtained from the TCGA-GBM and TCGA-LGG cohorts and downloaded from The Cancer Imaging Archive (). The pathologic images of TCGA-GBM and TCGA-LGG were downloaded from the TCGA genomic data commons (GDC) Data Portal (). Finally, 218 patients with TCGA gliomas (100 patients with GBM and 118 patients with LGG) were selected for further analysis.

We also enrolled 54 patients with glioma meeting the inclusion criteria from our medical group in Huashan Hospital, Fudan University, during 2020–2021. In addition, we screened patients with glioma (n = 57) in our previous 5hmC cohort (GSE132118) meeting the inclusion criteria. The present study with related protocols was approved by the ethics committee of the Huashan Hospital, Fudan University (KY2015-256). Written informed consents were obtained from all enrolled patients before surgery, which included a statement on the formalin-fixed paraffin-embedded (FFPE) samples and clinicopathological data for scientific research. All procedures involving human participants were in accordance with the Declaration of Helsinki.

MRI acquisition and processing

All patients underwent MRI examinations using a 3T Discovery750 MRI scanner (GE Healthcare) or a 3T Verio MRI scanner (Siemens Healthcare GmbH) randomly prior to surgical resection. Two major investigators (Y.F.Y. and Y.N.W. with 4 years of experience) blinded to clinical information manually delineated the tumor contours slice by slice on the axial section using 3D slicers. The delineated images were confirmed by another neurosurgeon (J.J Cai with 12 years of experience).

Whole side imaging (WSI) acquisition of H&E-stained slides

The H&E-stained slides of all enrolled patients were prepared by using FFPE samples. In the present study, one patient had one or multiple WSIs. The slides were scanned using the MAGSCAN-NER scanner (KF-PRO-005, KFBIO) to obtain the WSI.

Squeeze-and-excitation deep learning feature extractor (SE-DLFE) and deep learning image feature score (DLIFS)

We initially designed a feature extractor similar to a deep learning network,¹⁶ termed a DLFE, to extract the deep features from diverse medical images that efficiently reflected the prognostic capabilities. DLFE learned and extracted deep features from the input data through a neural network. Considering the differences between H&E slides and MRI scans, we first downsampled the two types of images separately to the same dimensions. Subsequently, modal fusion was performed, and a DLIFS was ultimately output through two fully connected layers. This feature score was considered to accurately reflect the prognostic index of the actual survival status of patients with glioma.

DLFE for MRI

The downsampling of the MRI scans employed a residual network incorporating an attention mechanism, similar to the Squeeze-and-Excitation Networks. The entire downsampling process consisted of five modules, each comprising two residual networks and one SE layer. The purpose of the residual network was to assist in preserving shallow-level feature information, thereby preventing gradient vanishing and enhancing the overall network performance. The SE-layer module modeled the inter-channel correlations of features at each layer and strengthened relevant features to improve accuracy. The specific implementation of the SE-layer unit is illustrated in Figure 1A.

[IMAGE OMITTED. SEE PDF]

DLFE for H&E-stained slides

The Otsu thresholding method was employed to segregate the white background from the tumor region due to the presence of a substantial amount of noninformative white background in the slides.¹⁷ Considering the nonnegligible size of the slides and the high labor cost associated with manual annotation, we aimed to capture sufficient information from the slides in the absence of predefined regions of interest (ROIs). We partitioned the slides into patches of fixed dimensions (512 × 512) at 20 × objective magnifications (0.5 µm per pixel). The patches were adjacent to one another, covering the entire WSI. Subsequently, these patches underwent processing using a pretrained ResNet50 model on ImageNet, resulting in the extraction of patch-level features characterized by dimensions of 1024 × 1. The application of the mean vector method to aggregate these patch-level features facilitated the derivation of slide-level features, manifesting dimensions of 1024 × the number of patches per slide.

Considering the potential variability in the number of slides from each patient, our analytical framework culminated in acquiring patient-level features (dimensions: 1024 × the number of slides from each patient). These features were deemed germane to the temporal aspect of survival analysis, thereby constituting a pivotal dimension within our investigative paradigm (Figure 1B). In this study, we used global average pooling as the squeeze operation. Following this, two fully connected layers formed a bottleneck structure to model inter-channel correlations, producing weights of 1/r for both output and input features. Subsequently, a sigmoid gate was employed to obtain normalized weights between 0 and 1, and a scaling operation was applied to weight each channel's feature using normalized weights. Simultaneously, we introduced a skip-connection module, adding the original data X to the feature weights obtained after the squeeze operation. This was primarily implemented to mitigate vanishing gradient problems.

In this section, we designed four Cox model sets receiving various inputs for further evaluation. The inputs included feature scores computed from various feature extractors and modalities. The inputs were as follows:

• -The feature score computed by the feature extractor constructed using MRI scans (squeeze-and-excitation DLFE [SE-DLFE]-MR].
• -The feature score computed by the feature extractor constructed using H&E-stained slides (SE-DLFE-H&E).
• -The feature scores computed by the feature extractors separately from two modalities using both SE-DLFE-MRI and SE-DLFE-H&E networks (SE-DLFE-M|H).
• -The feature score computed by the feature extractor constructed using both H&E-stained slides and MRI scans (SE-DLFE, which is the network designed in the present study).

Loss function

We employed the loss function of the CoxPH regression model for the feature extractor to ensure that the hidden layers could extract features highly relevant and robust for survival prediction and effectively handle missing data after dense connections in the intermediate layers (Figure 1C).

MRI radiomics approach

Briefly, the method first applied the package glmnet to construct radiomics features by least absolute shrinkage and selection operator (LASSO)–Cox regression. The fitted LASSO–Cox model was then applied to train and test patients, and a radiomics risk score was derived from a linear combination of model coefficient-weighted features.¹⁸ A survival analysis model was eventually constructed based on the risk score.

Pathomics approach

Based on the multi-instance learning framework, the attention mechanism was introduced to complete the feature aggregation of H&E-stained slides, which realized the transformation from the slide-level features to the patient-level features, which were then input into a convolutional neural network (CNN) to predict survival.¹⁹

Sample preparation, 5hmC-Seal profiling, and data processing

The details for preparing cfDNA, the 5hmC-Seal library establishment, the next-generation sequencing, and the data processing workflow are described in a previous study.⁹ Briefly, at least 5 mL of frozen plasma from each participant was collected from the peripheral blood, followed by cfDNA extraction and subsequent 5hmC-Seal profiling. The 5hmC-Seal count data were normalized by DESeq2,²⁰ which performed the variance-stabilizing transformation to correct the sequencing depth and library size.

Multimodal integration

In our dataset partitioning scheme, we adhered to a ratio of 8:2 for the training and test sets. Initially, we leveraged the SE-DLFE model to compute scores for image features, encompassing both pathologic and MRI data. Subsequently, we employed the LASSO–Cox analysis for patients with 5hmC sequencing to identify 25 gene sequences and calculate their corresponding gene scores. Building upon this, we constructed the support vector machine–MMP (SVM-MMP). We employed a late fusion approach to amalgamate the scores derived from these three modalities, expecting to achieve an enhanced prognostic outcome. Following model construction, we conducted a comprehensive performance analysis using pycox.evaluation and visualized the results through time receiver-operating characteristic (ROC) curves.

We applied the survival to construct our survival analysis model, whose advantage was that it could explain the complex nonlinear relationship between features and survival through kernel techniques. This allowed SVM to be described through hyperplanes, leading to extremely versatile survival SVMs applicable to a wide range of data. In the experiments, we used the HingeLossSurvivalSVM function in the pycox package, realizing its satisfactory predictive performance. Among these, the HingeLoss could obtain an optimal hyperplane by minimizing the objective function, resulting in correctly classified samples as far away from the hyperplane as possible and maximizing the interval. The output of the model was a risk score, and the higher the score, the greater the patient's risk.

Nomogram

A nomogram was delineated based on multivariate regression analysis,²¹ which integrated the risk scores calculated by the model with clinical variables (IDH1 mutation, sex, and WHO grade). Further, scaled line segments were drawn on the same plane at a certain scale to express the interrelationships between variables in the predictive model. In this study, we presented a nomogram of the multimodal survival analysis model by using the R package survival rms. We calculated the predicted value of the individual's outcome event using the relationship between the total score and the probability of outcome event occurrence.

Data availability

The imaging and clinical data could be requested by contacting the corresponding authors when necessary. The code utilized in the present study is publicly available at .

Statistical analysis

Statistics were mainly conducted by using the Python platform (version 3.6.1) and the R software (version 4.2.3). Continuous variables were compared using the independent-samples unpaired t-test to see whether they were normally distributed or using the Mann–Whitney U test to see whether they were nonnormally distributed. When comparing categorical variables, if at least one expected cell count was <5 in the 2 × C contingency table, Fisher's exact test was used; otherwise, the χ² test was performed. The Cox regression analysis was used for univariate and multivariate analyses, and the hazard ratio with a 95% confidence interval was calculated.

The follow-up duration was measured from the time of surgery to the last follow-up date, and the survival status at the last follow-up was recorded. The OS was defined as the interval between surgery and death or the last date of follow-up. The log-rank test was commonly used to assess the difference between two Kaplan–Meier (KM) survival curves. A log-rank p-value <.05 was considered statistically significant.²² The survival prediction performance of the methods was evaluated using Harrell's concordance index (C-index) for validating the prognostic performance of constructed models. This Harell's index is the most frequently used evaluation metric to assess the discriminative performance in survival analysis.²³ The prediction error curves over time and integrated Brier score on each dataset were calculated using the R package pec to assess the overall performance.²⁴ The time-dependent ROC analysis was performed to assess the prediction accuracy after 6, 12, 18, 24, 30, and 36 months for the risk score, where the area under the ROC curve (AUC) was calculated. Calibration curves were generated to evaluate performance in both the model training and test cohorts. Decision curve analysis (DCA) was conducted to assess the clinical utility of the relevant model.

RESULTS

Baseline clinicopathologic characteristics of enrolled patients with glioma

We examined the complementary prognostic information of multimodal features derived from the clinical, radiologic, histopathologic, and genomic data obtained during the clinical workup of patients with glioma. The schematic outline of the present study is demonstrated in Figure 2. Overall, 218 patients with glioma (1308 MRI scans and 396 H&E-stained slides) from the two TCGA cohorts, TCGA-GBM and TCGA-LGG datasets, were included in the present study as mentioned above. These patients were divided into the TCGA training set (n = 153, 918 MRI scans and 272 H&E-stained slides), the TCGA validation set (n = 21, 126 MRI scans and 46 H&E-stained slides), and the TCGA internal test set (n = 44, 264 MRI scans and 78 H&E-stained slides). Besides, we included patients with glioma (n = 54, 324 MRI scans and 99 H&E-stained slides) from the Huashan Hospital as an external cohort. The baseline clinicopathologic characteristics of the two cohorts are summarized in Table 1. No significant differences in age, sex, WHO grade, tumor location, and IDH mutation status were observed among the training set, validation cohort, internal test set, and external test set (two-sided Wilcoxon test or χ² test p-value >.05).

[IMAGE OMITTED. SEE PDF]

TABLE 1 Baseline clinicopathologic characteristics of enrolled patients.

	TCGA glioma cohort
Baseline characteristics	Total (n = 218)	Training set (n = 153)	Validation set (n = 21)	Internal test set (n = 44)	Huashan external test set (n = 54)	p-value
Age, year (median IQR)	52 (37–62)	53 (37–61)	50 (39–63)	50 (39–62)	42 (36–58)	.301
Male, n (%)	125 (57%)	85 (56%)	13 (62%)	27 (61%)	34 (63%)	.744
WHO grade, n (%)						.108
LGG	118 (54%)	78 (51%)	10 (48%)	30 (68%)	24 (42%)
GBM	100 (46%)	75 (49%)	11 (52%)	14 (32%)	30 (58%)
Tumor location, n (%)						.200
MLFA	70 (32%)	50 (33%)	8 (38%)	12 (27%)	25 (46%)
NMLFA	148 (68%)	103 (67%)	13 (62%)	32 (73%)	29 (54%)
IDH mutation, n (%)	126 (58%)	87 (57%)	13 (62%)	26 (59%)	17 (31%)	.070

Prognostic prediction efficacy of SE-DLFE stratification

The construction of effective feature extractors and prognostic predictors from medical images was one of the main objectives of the present study. We innovatively established SE-DLFE and DLIFS to process images. The class activation mapping based on T1-CE MRI was constructed from the feature weight of the last residual model and the SE model in the network (Supporting Information Figure S1). The features strongly associated with lesion areas were extracted from T1-CE sequences and H&E-stained slides, as well as the survival time, to obtain DLIFS of patients with glioma. All patients with glioma in the TCGA glioma training set, TCGA glioma internal test set, and Huashan external test set were obviously stratified into a high-risk group (DLIFS ≥ 0.399) and a low-risk group (DLIFS < 0.399) using median DLIFS obtained from the TCGA glioma training set as the cutoff, indicating the prognostic prediction efficiency of DLIFS (Figure 3A–C).

[IMAGE OMITTED. SEE PDF]

The C-index and Brier score were applied and calculated to evaluate the prognostic prediction efficacy of the four established SE-DLFE models (Supporting Information Table S1). We noted that SE-DLFE demonstrated the highest predictive efficacy in the TCGA glioma internal test set (C-index = 0.929; Brier score = 0.060) and Huashan external test set (C-index = 0.800; Brier score = 0.164) among the four models (Figures 3D–F and S2). In addition, we delineated 6- to 36-month time-dependent AUC and loss decay plot (Figure 3G,H), which indicated the model stability.

Application for subgroup analysis of patients with glioma

We applied the following models for subgroup analysis to confirm the survival prediction efficacy and accuracy of established SE-DLFE models for patients with glioma of multiple subtypes. First, we testified the performance of established models in patients with glioma of various WHO grades (Table 2). Four models based on SE-DLFE consistently demonstrated relatively lower efficiency for GBM, compared with LGG, which was attributed to the heterogeneity and intricacy of GBM itself. Among the four SE-DLFE-based models (Figure 4A–D), SE-DLFE possessed the highest efficiency for the TCGA-LGG internal test set (C-index = 0.953; Brier score = 0.080) and the Huashan-LGG external test set (C-index = 0.818; Brier score = 0.153), as well as a satisfactory performance for the TCGA-GBM internal test set (C-index = 0.861; Brier score = 0.062) and the Huashan-GBM external test set (C-index = 0.693; Brier score = 0.192).

TABLE 2 Prognostic prediction efficacy of the four squeeze-and-excitation deep learning feature extractor (SL-DLFE) models for patients with glioma with various World Health Organization (WHO) grades.

Prediction indicators	TCGA-LGG internal test set (n = 30)	TCGA-GBM internal test set (n = 14)	Huashan-LGG external test set (n = 24)	Huashan-GBM external test set (n = 30)
C-index
SE-DLFE-MRI	0.929	0.833	0.769	0.577
SE-DLFE-H&E	0.835	0.653	0.776	0.716
SE-DLFE-M|H	0.953	0.806	0.774	0.588
SE-DLFE	0.953	0.861	0.818	0.693
Brier score
SE-DLFE-MRI	0.177	0.108	0.126	0.251
SE-DLFE-H&E	0.194	0.138	0.164	0.162
SE-DLFE-M|H	0.191	0.105	0.141	0.252
SE-DLFE	0.080	0.062	0.153	0.192

[IMAGE OMITTED. SEE PDF]

In clinical practice, preserving neurologic function was often a primary concern, particularly for patients with functional area gliomas. Considering the impact of functional protection on the conservative extent of surgical resection, we also confirmed the prediction efficiency of constructed models in patients with MLFA gliomas (Table 3 and Figure 4E,F). We noted that SE-DLFE could maintain the highest prediction accuracy for MLFA gliomas (C-index = 0.842; Brier score = 0.095). No obvious performance difference was demonstrated between motor functional area (MFA) (C-index = 0.860; Brier score = 0.094) and language functional area (LFA) (C-index = 0.842; Brier score = 0.133) (Supporting Information Figure S3).

TABLE 3 Prognostic prediction efficacy of four SL-DLFE models for patients with gliomas in motor and language functional area (MLFA) and non-MLFA (NMLFA).

	MLFA
Prediction indicators	Total (n = 37)	MFAa (n = 24)	LFAa (n = 16)	NMLFA (n = 61)
C-index
SE-DLFE-MRI	0.832	0.884	0.825	0.785
SE-DLFE-H&E	0.810	0.837	0.775	0.773
SE-DLFE-M|H	0.824	0.884	0.783	0.783
SE-DLFE	0.842	0.860	0.842	0.855
Brier score
SE-DLFE-MRI	0.110	0.081	0.134	0.156
SE-DLFE-H&E	0.120	0.127	0.144	0.172
SE-DLFE-M|H	0.129	0.113	0.142	0.142
SE-DLFE	0.095	0.094	0.133	0.091

Prognostic value of genome-wide 5hmC profiles in cfDNA for patients with glioma

The genome-wide 5hmC profiles of our glioma cohort were previously reported and uploaded in GSE132118.⁹ We analyzed 234 differentially expressed 5hmC-related features (231 upregulated and three downregulated features) by comparing 111 patients with glioma and 111 controls (Figure 5A,B). We selected 35 genes (Supporting Information Table S2) associated with OS (p <.05) by univariate Cox regression analysis or KM plotter in 99 patients with glioma with complete 5hmC expression profiling and survival data to identify a 5hmC-related signature. Then, we obtained 25 genes as active covariates with considerable prognostic value (Supporting Information Table S3) using the LASSO regression algorithm (Figure 5C). The risk scores for the patients were calculated and demonstrated (Figure 5D). We divided the patients with glioma into high-risk subgroup (n = 49) and low-risk subgroup (n = 50) using the median risk score as a cutoff value to assess the performance of the 5hmC signature as a classifier. The results revealed a remarkable younger age and poorer OS in high-risk subgroups (Figure 5E and Table 4).

[IMAGE OMITTED. SEE PDF]

TABLE 4 Correlation between 5-hydroxymethylcytosines (5hmC) signature-based risk scores and clinicopathologic characteristics of Huashan 5hmC glioma cohort.

	Huashan glioma cohort with 5hmC profiling and survival data (n = 99)
Characteristics	Low-risk score (n = 50)	High-risk score (n = 49)	p-value
Age, year	51 (42, 61)	47 (35, 52)	.038
Male, n (%)	33 (66%)	35 (71%)	.560
WHO grade, n (%)			.122
LGG	14 (28%)	21 (43%)
GBM	36 (72%)	28 (57%)
Tumor location, n (%)			.667
MLFA	12 (24%)	10 (20%)
Other areas	38 (76%)	39 (80%)
IDH mutation, n (%)	14 (28%)	21 (43%)	.122

Multimodal prognostication and validation

Multimodal integration of advanced genetic molecular techniques, radiologic and histologic imaging, and codified clinical data presents opportunities to predict precise oncology beyond single-level modal data.⁸ We screened and obtained 57 patients with complete MR T1-CE images, H&E-stained slides, 5hmC sequencing profiles, OS, and clinical information in the Huashan glioma cohort and subsequently divided them into the training set (n = 45) and test set (n = 12; Supporting Information Table S4). Subsequently, we calculated and demonstrated the prognostic prediction efficacy of three separate multimodalities integrated with clinical data using the SVM method (Supporting Information Table S5), including the traditional MR radiomic model, H&E pathologic model, and SE-DLFE model (Figure 6A), 5hmC signature model (Figure 6B), and multimodal predictor (MMP) (Figure 6C). In addition, we also compared our models with previously conventional MR radiomic model¹⁸ and H&E pathologic model¹⁹ to confirm the accuracy superiority. We observed that SVM-MMP possessed the highest prediction efficacy (C-index = 0.867; Brier score = 0.185) for OS of patients with glioma. Time-dependent AUC also demonstrated the model accuracy in various time nodes (Figure 6D). All SVM-SE-DLFE, SVM-5hmC, and SVM-MMP could remarkably stratify the patients with glioma by OS (Figures 6E and S4A,B). In addition, the nomograms of discrete multimodal patterns were delineated to evaluate and predict the probability of patients with glioma (Figures 6F and S4C,D).

[IMAGE OMITTED. SEE PDF]

To validate the reliability of nomogram, calibration curves were demonstrated for both Huashan training set. The calibration curve of the SVM-MMP in the training cohort indicates that the 1, 2, and 3-year curves are all near the 45° diagonal line, indicating satisfactory consistency between the nomogram's predicted survival rates and the actual probabilities (Figure 7A). This suggests that multimodal risk score is significantly associated with 1, 2, and 3-year survival rates of patients. The DCA identified the satisfactory prediction efficacy and discriminatory ability of the nomogram for patients’ 1, 2, and 3-year survival status (Figure 7B). In summary, the established nomogram for predicting OS rate has strong discriminative power and high calibration accuracy.

[IMAGE OMITTED. SEE PDF]

DISCUSSION

Despite the advancements in therapeutic strategies such as surgical resection, chemoradiotherapy, and combination treatments, gliomas are generally correlated with a poor prognosis probably on account of the inherent heterogeneity across individuals. Accurate and robust survival predictions for patients with gliomas can provide valuable guidance for diagnosis, treatment planning, and outcome prediction. An improved understanding of tumor characteristics through multimodal DLs, drawing on multiple types of data from each patient, can lead to a more accurate prediction of OS, compared with unimodal models, contributing to a further optimal treatment plan, including inclusion criteria of clinical trials. Recently, a systematic review summarized the current progression of multimodal deep-learning-based prognostication in patients with glioma.⁴ Among these multimodal progressions, Braman et al. applied the most comprehensive modalities including MRI, WSI, genomic, biomarker, and clinical data to predict the OS of gliomas using a pretrained VGG-19 CNN,²⁵ achieving a median C-index of 0.788. Therefore, the prediction accuracy and efficiency have much room for improvement.

In the present study, the feature extraction process involved the direct analysis of patches cropped from H&E-stained slides, resulting in the derivation of patch-level features. However, recognizing the value of patient-level features for survival analysis, we sought to integrate information from multiple imaging modalities, including both H&E-stained slides and MRI scans. This multimodal approach presented a unique challenge, considering the inherent differences between the two imaging techniques. Our feature extractor seamlessly incorporated the information gleaned from both H&E-stained slides and MRI scans to address this challenge and create a cohesive analysis platform. The integration process involved the extraction of features from MRI scans, as well as the aggregation of features extracted from both modalities. For H&E-stained slides specifically, our method stood out for its adaptability in the absence of predefined ROIs. Instead of relying on predefined areas,²⁵ we strategically cropped each slide into multiple patches with dimensions of 512 × 512. This not only boosted computational efficiency but also proved highly effective in extracting prognostically relevant features. Our method excelled in scenarios where ROIs were undefined, allowing for a comprehensive exploration of the entire slide. By embracing this patch-based strategy, we successfully navigated the challenges posed by the absence of ROIs, ensuring a robust and insightful analysis of histologic slides with minimal manual intervention and optimal computational efficiency. Our approach combined information from distinct imaging modalities and emphasized the equal contribution of all patches or slides during the aggregation process. This was achieved by the mean vector method,¹⁹ preserving the individual impact of each patch or slide on the final patient-level features. By integrating H&E and MRI features within a unified feature extraction framework, our approach not only addressed the need for patient-level features in survival analysis but also leveraged the distinct advantages offered by both imaging modalities.

Besides improvements in imaging modalities, we also attempted to explore new biomarkers for glioma survival. The depletion of 5hmC is generally associated with the hypermethylation of gene bodies in various cancers such as glioma,^9,26 hepatocellular carcinoma,¹⁰ and colorectal cancer, implicating its predictive role in cancer diagnosis and therapy. The genome-wide 5hmC profiles of tumor tissue have been verified to exhibit a superior predictive performance, compared with IDH1 mutation status for GBM prognosis,²⁶ but it also possessed several limitations such as precision loss caused by the heterogeneity of sampling sites. Recently, the sensitivity and accuracy of liquid biopsy have been improved by applying high-throughput sequencing technologies. Modern liquid biopsy techniques have not only enabled the early detection and diagnosis of various cancers but also provided a sensitive modality for the noninvasive evaluation of disease status and prognosis. Technically, the potential value of the cfDNA 5hmC-Seal approach has emerged in cancer diagnosis,^9,10,12 treatment response prediction,²⁷ and prognostic implications.¹⁴ We previously reported a noninvasive diagnostic strategy through an integrative analysis of genome-wide 5hmC profile using a highly sensitive 5hmC-Seal technique in cfDNA samples,⁹ which has been confirmed in other studies.²⁸ In this study, we first selected 25 prime 5hmC-related signatures and validated the prognostic prediction capacity of cfDNA-based genome-wide 5hmC signatures in patients with glioma.

No study to date has integrated radiology, pathology, and cfDNA-based liquid biopsy within a deep learning framework for outcome prediction or patient stratification. In this study, we first implemented the integration of multimodal data including T1-CE images of MRI, histologic images, cfDNA-based genome-wide 5hmC signature, and clinical variables to achieve satisfactory performance for glioma survival prediction (C-index = 0.897; Brier score = 0.118). We also established SE-DLFE to process medical images in large numbers to achieve multimodal integration; data fusion schemes must be highly efficient in learning complex multimodal interactions. Second, we applied the SVM method for integration as it possesses more compatibility for small samples and nonlinear data.²⁹

This study had several limitations. First, the sample size of the enrolled cohort was relatively small mainly because of the challenge of incomplete modality fusion especially for 5hmC-sequencing profiles, which was also a potential barrier to the clinical implementation of multimodal deep learning methods. Therefore, it is imperative to popularize the application of 5hmC-related signature panel or delve into the survival analysis of patients lacking a specific modality, a facet that warrants further discussion.³⁰ Second, the glioma region was manually identified and delineated in MRI scans because the current algorithms could not automatically recognize the precise margin of gliomas. This time-consuming procedure led to a limited number of eligible samples. Third, the adjuvant therapy scheme was not included in our established model, which could partially influence postoperative outcomes and survival, especially for LGG.^31,32 In addition, the follow-up duration of the Huashan external test set was relatively shorter than TCGA glioma cohorts. This is limited by data property because the glioma patients from the TCGA dataset are enrolled in a relatively earlier period.

In conclusion, we constructed a deep learning framework (CoxPH-MMP) integrating radiologic, histologic, cfDNA-based epigenetics, and clinical variables into a fused prognostic prediction model to stratify patients with glioma by OS for further clinical decision-making, such as clinical trial enrollment and application of multidisciplinary strategies.

AUTHOR CONTRIBUTIONS

Yifan Yuan and Yining Wang collected clinical imaging and pathologic slides, conducted statistical analysis, and wrote the manuscript. Xuan Zhang and Hongyan Li processed the images and data for algorithm construction. Zengxin Qi, Zunguo Du, and Danyang Feng helped in data and sample collection. Qingguo Xie and Ying-Hua Chu provided protocol guidance of this study. Jie Hu, Yuqing Liu, and Jianping Song sponsored the research. Jiajun Cai, Yuqing Liu, and Jianping Song conceived the project and revised the manuscript.

ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China under grant number 62131009, the University Synergy Innovation Program of Anhui Province (GXXT-2022-032), Fujian Province Science and Technology Innovation Joint Fund (No. 2021Y9135), Science and Technology Innovation Plan of Shanghai Science and Technology Commission (21Y21900600) and Shanghai Municipal Alliance for Clinical Competence Improvement and Advancement in Neurosurgery (SHDC22021303).

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

References

H. Sung, J. Ferlay, R. L. Siegel, M. Laversanne, I. Soerjomataram, A. Jemal, F. Bray, CA Cancer J. Clin. 2021, 71, 209.

T. A. Dolecek, J. M. Propp, N. E. Stroup, C. Kruchko, Neuro. Oncol. 2012, 14, v1.

P. Yang, Y. Wang, X. Peng, G. You, W. Zhang, W. Yan, Z. Bao, Y. Wang, X. Qiu, T. Jiang, J. Neurooncol. 2013, 113, 259.

K. Alleman, E. Knecht, J. Huang, L. Zhang, S. Lam, M. Decuypere, Cancers 2023, 15, 545.

K. M. Boehm, E. A. Aherne, L. Ellenson, I. Nikolovski, M. Alghamdi, I. Vázquez‐García, D. Zamarin, K. Long Roche, Y. Liu, D. Patel, A. Aukerman, A. Pasha, D. Rose, P. Selenica, P. I. Causa Andrieu, C. Fong, M. Capanu, J. S. Reis‐Filho, R. Vanguri, H. Veeraraghavan, N. Gangai, R. Sosa, S. Leung, A. Mcpherson, J. Gao, Y. Lakhman, S. P. Shah, Nat. Cancer 2022, 3, 723.

J. Zhou, A. Foroughi Pour, H. Deirawan, F. Daaboul, T. N. Aung, R. Beydoun, F. S. Ahmed, J. H. Chuang, EBioMedicine 2023, 94, [eLocator: 104726].

J. Goddard, J. Castle, E. Southworth, A. Fletcher, S. Crosier, I. Martin‐Guerrero, M. García‐Ariza, A. Navajas, J. Masliah‐Planchon, F. Bourdeaut, C. Dufour, O. Ayrault, T. Goschzik, T. Pietsch, M. Sill, S. M. Pfister, S. Rutkowski, S. Richardson, R. M. Hill, D. Williamson, S. Bailey, E. C. Schwalbe, S. C. Clifford, D. Hicks, Acta Neuropathol. 2023, 145, 651.

K. M. Boehm, P. Khosravi, R. Vanguri, J. Gao, S. P. Shah, Nat. Rev. Cancer 2022, 22, 114.

J. Cai, C. Zeng, W. Hua, Z. Qi, Y. Song, X. Lu, D. Li, Z. Zhang, X. Cui, X. Zhang, Z. Yang, J. Zhang, K. Quan, W. Zhu, J. Cai, C. He, S. Y. Cheng, W. Zhang, Y. Mao, Neurooncol Adv. 2021, 3, [eLocator: vdab049].

J. Cai, L. Chen, Z. Zhang, X. Zhang, X. Lu, W. Liu, G. Shi, Y. Ge, P. Gao, Y. Yang, A. Ke, L. Xiao, R. Dong, Y. Zhu, X. Yang, J. Wang, T. Zhu, D. Yang, X. Huang, C. Sui, S. Qiu, F. Shen, H. Sun, W. Zhou, J. Zhou, J. Nie, C. Zeng, E. K. Stroup, X. Zhang, B. C.‐H. Chiu, W. Y. Lau, C. He, H. Wang, W. Zhang, J. Fan, Gut 2019, 68, 2195.

G. D. Guler, Y. Ning, C.‐J. Ku, T. Phillips, E. Mccarthy, C. K. Ellison, A. Bergamaschi, F. Collin, P. Lloyd, A. Scott, M. Antoine, W. Wang, K. Chau, A. Ashworth, S. R. Quake, S. Levy, Nat. Commun. 2020, 11, 5270.

J. Shao, S. Wang, D. West‐Szymanski, J. Karpus, S. Shah, S. Ganguly, J. Smith, Y. Zu, C. He, Z. Li, Sci. Rep. 2022, 12, [eLocator: 12410].

C. Zeng, E. K. Stroup, Z. Zhang, B. C. Chiu, W. Zhang, Cancer Commun. 2019, 39, 12.

B. C.‐H. Chiu, Z. Zhang, Q. You, C. Zeng, E. Stepniak, P. M. Bracci, K. Yu, G. Venkataraman, S. M. Smith, C. He, W. Zhang, Blood Adv. 2019, 3, 2790.

Z. Cai, J. Zhang, Y. He, L. Xia, X. Dong, G. Chen, Y. Zhou, X. Hu, S. Zhong, Y. Wang, H. Chen, D. Xie, X. Liu, J. Liu, ESMO Open 2021, 6, [eLocator: 100021].

J. Yan, Q. Sun, X. Tan, C. Liang, H. Bai, W. Duan, T. Mu, Y. Guo, Y. Qiu, W. Wang, Q. Yao, D. Pei, Y. Zhao, D. Liu, J. Duan, S. Chen, C., W. Wang, Z. Liu, X. Hong, X. Wang, Y. Guo, Y. Xu, X. Liu, J. Cheng, Z.‐C. Li, Z. Zhang, Eur. Radiol. 2023, 33, 904.

N. Otsu, IEEE Trans. Syst. Man Cybern. 1979, 9, 62.

B. J. Geraghty, A. Dasgupta, M. Sandhu, N. Malik, P. J. Maralani, J. Detsky, C.‐L. Tseng, H. Soliman, S. Myrehaug, Z. Husain, J. Perry, A. Lau, A. Sahgal, G. J. Czarnota, J. Neurooncol. 2022, 156, 579.

J. Yao, X. Zhu, J. Jonnagaddala, N. Hawkins, J. Huang, Med. Image Anal. 2020, 65, [eLocator: 101789].

M. I. Love, W. Huber, S. Anders, Genome Biol. 2014, 15, 550.

G. H. Bell, A. C. Paterson, Arch. Dis. Child. 1930, 5, 369.

V. Bewick, L. Cheek, J. Ball, Crit. Care 2004, 8, 389.

F. E. Harrell, K. L. Lee, R. M. Califf, D. B. Pryor, R. A. Rosati, Stat. Med. 1984, 3, 143.

U. B. Mogensen, H. Ishwaran, T. A. Gerds, J. Stat. Softw. 2012, 50, 1.

N. Braman, J. W. H. Gordon, E. T. Goossens, C. Willis, M. C. Stumpe, J. Venkataraman. Deep Orthogonal Fusion: Multimodal Prognostic Biomarker Discovery Integrating Radiology, Pathology, Genomic, and Clinical Data. Springer International Publishing, Cham 2021.

C. Zeng, X. Song, Z. Zhang, Q. Cai, J. Cai, C. Horbinski, B. Hu, S.‐Y. Cheng, W. Zhang, Acta Neuropathol. Commun. 2023, 11, 133.

Q. Li, C. C. Huang, S. Huang, Y. Tian, J. Huang, A. Bitaraf, X. Dong, M. T. Nevalanen, J. Zhang, B. J. Manley, J. Y. Park, M. Kohli, E. M. Gore, D. Kilari, L. Wang, medRxiv 2023 [DOI: https://dx.doi.org/10.1101/2023.10.13.23296758]

S. Zhang, J. Zhang, X. Hu, S. Yin, Y. Yuan, L. Xia, F. Cao, X. Yan, Z. Yan, Q. Mao, D. Xie, Y. Liu, Int. J. Cancer 2023, 152, 1707.

V. V. Belle, K. Pelckmans, J. A. K. Suykens, S. V. Huffel, Artif. Intell. Med. 2011, 53, 107.

C. Cui, H. Liu, Q. Liu, R. Deng, Z. Asad, Y. Wang, S. Zhao, H. Yang, B. A. Landman, Y. Huo, Survival Prediction of Brain Cancer with Incomplete Radiology, Pathology, Genomic, and Demographic Data. Springer Nature, Switzerland 2022.

X.‐P. Liu, X. Jin, S. Seyed Ahmadian, X. Yang, S.‐F. Tian, Y.‐X. Cai, K. Chawla, A. M. Snijders, Y. Xia, P. J. Van Diest, W. A. Weiss, J.‐H. Mao, Z.‐Q. Li, H. Vogel, H. Chang, Neuro. Oncol. 2023, 25, 68.

J. Hayes, Y. Yu, L. E. Jalbert, T. Mazor, L. E. Jones, M. D. Wood, K. M. Walsh, H. Bengtsson, C. Hong, S. Oberndorfer, T. Roetzer, I. V. Smirnov, J. L. Clarke, M. K. Aghi, S. M. Chang, S. J. Nelson, A. Woehrer, J. J. Phillips, D. A. Solomon, J. F. Costello, Neuro. Oncol. 2018, 20, 632.