Machine learning derived development and validation of extracellular matrix related signature for predicting prognosis in adolescents and young adults glioma

Integrated construction of MLDPS

The workflow of our study was illustrated in Fig. 1. With |log₂FC| > 1 and adjusted P value < 0.05, 6,401 differentially expressed genes (DEGs) were screened out between AYAs glioma and control. After intersecting with 1,026 ECM-related genes, 508 ECM-related DEGs were identified. Next, we obtained 361 overlapped genes existed in ECM-related DEGs, TCGA, CGGA-693 and CGGA-325 cohorts. Different from previous studies, we hypothesized that each cohort had the potential to generate the optimal prognostic model when treated as the training cohort. Therefore, our study proposed an innovative circuit training and validation procedure to the machine learning workflow, which means when one cohort utilized to training the model, others were used for validation. We identified 104, 184 and 185 prognostic genes in TCGA, CGGA-693 and CGGA-325 cohorts, respectively, by univariate Cox analysis. Subsequently, we applied 65 machine learning algorithm combinations to develop the prognostic models with ten-fold cross-validation and calculated C-index for each algorithm in all cohorts. The highest average C-index in CGGA-693 training, CGGA-325 and TCGA training cohorts were 0.828 (0.823 in TCGA, 0.840 in CGGA-693 and 0.821 in CGGA-325, Fig. 2A), 0.818 (0.809 in TCGA, 0.774 in CGGA-693 and 0.871 in CGGA-325, Fig. 2B) and 0.829 (0.909 in TCGA, 0.754 in CGGA-693 and 0.824 in CGGA-325, Fig. 2C), respectively. Top five average C-index in each training cohort were shown in Fig. 2D-F. Obviously, overfitting was detected in the highest mean C-index in TCGA training session, with C-index of 0.909 in TCGA, while a C-index less than 0.8 (0.754) in external validation CGGA-693 cohort. Consequently, we selected the second highest average C-index derived from Ridge algorithm in CGGA-693 training session as the optimal model on account of all C-index were more than 0.80 in this model and defined it as MLDPS. The prognostic genes included in MLDPS and the formula for calculating MLDPS score were summarized in supplementary Table 3 and Table 4, respectively.

Robust and consistent predictive performance of MLDPS

In clinical practice and management, several clinical and molecular characteristics such as grade, IDH status, age, gender and 1p/19q status have been applied for designing treatment regimens and prognosis prediction in glioma. Therefore, we contrasted the predictive performance of MLDPS with these features in the training and validation cohorts. As shown in Fig. 2G-I, the C-index of MLDPS was obviously higher than other features, suggesting MLDPS had obvious improved accuracy in prognosis prediction.

In addition, patients were dichotomized into high and low MLDPS groups according to the median MLDPS score (supplementary Table 5). The Kaplan-Meier curves showed that patients in high MLDPS group had obviously dismal OS compared with low MLDPS group in CGGA-693 cohort (p < 0.0001, Fig. 3A). Similar results were observed in CGGA-325 (p < 0.0001, Fig. 3B) and TCGA cohorts (p = 2e-04, Fig. 3C). Additionally, univariate Cox analysis revealed that MLDPS was a risky prognostic factor in CGGA-693 cohort (HR: 4.609 [3.647–5.824], p < 0.001, Fig. 3D). After adjusting for common clinical characteristics like grade, IDH, 1p/19q status and age (p < 0.05), MLDPS still remained a remarkably risky factor for prognosis in CGGA-693 cohort (HR: 5.091 [3.779–6.868], p < 0.001, Fig. 3D). Consistently, the analyses in CGGA-325 and TCGA cohorts demonstrated that MLDPS could serve as an independent prognostic factor for AYAs glioma (Fig. 3E-F). Furthermore, we also evaluated the predictive performance of MLDPS using ROC analysis. The areas under ROC curve (AUC) for 1-, 3- and 5-year survival were 0.851, 0.898 and 0.934, respectively, in CGGA-693, indicating that MLDPS owns robust performance in training cohort (Fig. 3G). Additionally, similar results were witnessed in two validation cohorts, including 0.869, 0.904 and 0.924 in CGGA-325 cohort (Fig. 3H) and 0.984, 0.928 and 0.731 in TCGA cohort (Fig. 3I). Taken together, these results indicated MLDPS possessed a robust and stable performance in predicting prognosis across different independent AYAs glioma cohorts.

The clinical significance of MLDPS

The results of subgroup analyses indicated that patients aged in 30–39 had lower MLDPS compared with younger patients in CGGA-693 (Fig. 4A) and CGGA-325 cohorts (Fig. 4B), while no differences were found in TCGA cohort (Fig. 4C). Besides, patients with IDH-wildtype, higher grade and 1p/19q non-codeletion had higher MLDPS in all cohorts (Fig. 4A-C). However, there were no differences in gender between high and low MLDPS groups (Fig. 4A-C). In addition, we also performed stratification survival analysis in different subgroups using Kaplan-Meier method except for the subgroups with only few patients like age 15–19 subgroups. As shown in Fig. 4D, in different groups, such as age 20–29, 30–39, male, female, IDH mutant, IDH-wildtype, WHO II/III, WHO IV, 1p/19q codeletion and non-codeletion subgroups, patients with high MLDPS had significantly worse OS than low groups in CGGA-693 cohort (all p < 0.05). The Kaplan-Meier curves in CGGA-325 and TCGA cohorts were similar with the above results (supplementary Fig. 2). These findings indicated that high MLDPS was associated with worse clinical behavior in AYAs glioma.

MLDPS outperforms previous 89 published prognostic signatures

The rapid development in high-throughput sequencing has greatly facilitated the precise treatment and stratified management for patients with cancer. Numerous prognostic signatures in glioma have been developed through different algorithms such as univariate Cox analysis and Lasso algorithm based on RNA-seq or microarray data among different cohorts. However, to the best of our knowledge, we found no studies focused on prognostic signatures in AYAs glioma. The available published prognostic signatures in glioma always focused on patients across all age groups. Considering that the authors concluded that their models could predict the prognosis of patients with glioma, which included the AYAs group in our study, we decided to collect these published models focused on both lower grade gliomas and glioblastoma multiforme. Finally, we comprehensive collected 89 published mRNA prognostic signatures to compare their predictive performance in AYAs glioma with MLDPS. The detailed information of these signatures was summarized in supplementary Table 6.

We calculated C-index for each published signature across all cohorts and made comparison with MLDPS via compare C package. Our MLDPS exhibited the highest mean C-index 0.828 than other signatures across the three cohorts (Fig. 5A). Additionally, the C-index of MLDPS ranked first in both CGGA-693 (Fig. 5B) and CGGA-325 cohorts (Fig. 5C). In TCGA cohort, MLDPS ranked the fifth among all signatures, while the top four displayed poor performance in CGGA-693 and CGGA-325 cohorts (Fig. 5D). For example, the signatures presented by Liu et al.⁵¹ and Zhang et al.⁵² ranked the first and second in TCGA cohort, whereas their C-index in CGGA-693 were both less than 0.7. This overfitting in models might weaken the generalization power for clinical practice^30,34. Furthermore, we also compared the AUC values of 1-, 3- and 5-year between MLDPS and 89 published signatures. In CGGA-693 cohort, MLDPS exhibited the highest AUC values in predicting overall survival for 1-year (Fig. 5E), 3-year (Fig. 5F) and 5-year (Fig. 5G). Moreover, the AUC values of MLDPS ranked the first or top in CGGA-325 and TCGA cohorts (supplementary Fig. 3). In summary, the above results demonstrated that MLDPS possessed a distinctly superior performance and better extrapolation potential than other prognostic signatures.

Moreover, given the robust predictive performance of MLDPS in AYAs glioma, we additionally evaluated its prognostic value across pan-cancer level. The process included the evaluation of MLDPS in AYAs cancers and the evaluation of MLDPS in cancers with all age groups. After preprocessing, merging and eliminating patients without survival information, OS less than 30 days (these patients may die due to lethal complication such as severe infection and hemorrhage^53,54, overlapped patients with our AYAs glioma cohort, 9,062 patients from 32 cancer types were included for survival analysis. There were six cancer types with more than 50 patients in AYAs including breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), pheochromocytoma and paraganglioma (PCPG), skin cutaneous melanoma (SKCM), testicular germ cell tumors (TGCT) and thyroid carcinoma (THCA). Interestingly, there were no evident association between MLDPS and the prognosis in theses cancers (Supplementary Fig. 4A-F), which might be due to the small sample sizes in these cancer types. Additionally, in cancers with all age groups, the Kaplan-Meier curves indicated that patients in high MLDPS groups exhibited dismal prognosis in adrenocortical carcinoma (ACC, p = 0.0078, Fig. 6A), BRCA, (p = 0.0034, Fig. 6B), colon adenocarcinoma (COAD, p = 0.0057, Fig. 6C), glioma at other age groups (p < 0.0001, Fig. 6E), mesothelioma (MESO, p = 0.012, Fig. 6F), sarcoma (SARC, p = 0.011, Fig. 6G) and THCA (p = 0.048, Fig. 6H). While intriguingly, the opposite trend with strong tendency was observed in patients with acute myeloid leukemia (LAML, p = 0.05, Fig. 6D), which was probably due to the compositions of ECM in blood cancer differed significantly from those in solid cancers. The above results suggested that MLDPS had potential for generalization to other cancer types.

The potential biological functions of MLDPS

GSEA analysis was utilized to elucidate the potential signaling pathways and biological processes between different MLDPS groups. As shown in Fig. 6I-J, tumor aggressiveness-related biological pathways were enriched in the high MLDPS group, such as DNA replication initiation, collagen biosynthetic process and formation, NF-kB pathway, cell cycle and focal adhesion. While the low MLDPS group was significantly correlated with metabolism-related functions like creatine metabolism, purine catabolism and oxidative phosphorylation (Fig. 6K-L). Additionally, we also investigated the potential biological functions related to MLDPS. The results indicated that MLDPS was remarkably enriched in cell migration, cell motility, cell adhesion, angiogenesis, extracellular structure organization processes, PI3K-Akt signaling pathway, TGF-beta signaling pathway, ECM-receptor interaction and focal adhesion (supplementary Fig. 4G-H).

Immune landscape correlated with MLDPS

To elucidate the connection between immune landscape and MLDPS in AYAs glioma, we explored the relationship between MLDPS and different immunity-related indexes. Firstly, we used the ESTIMATE algorithm to calculate immune scores, stromal scores and tumor purity. As shown in Fig. 7A, the high MLDPS group possessed higher ESTIMATE, immune and stromal scores compared to the low MLDPS group in CGGA-693 cohort. The same results were observed in CGGA-325 (p < 0.05, Fig. 7B) and TCGA cohorts (p < 0.05, Fig. 7C), respectively. These results implied that high MLDPS group owned more infiltration in both stromal and immune cells, which contributes to more complex TME. Additionally, the analyses in tumor purity revealed that high MLDPS groups had remarkably lower tumor purity than low MLDPS groups in CGGA-693 (p < 0.05, Fig. 7A), CGGA-325 (p < 0.05, Fig. 7B) and TCGA cohorts (p < 0.05, Fig. 7C), respectively, which were in harmony with the analyses in above TME scores.

Subsequently, we utilized the ssGSEA algorithms to assess tumor-infiltrating immune cells toward patients with different grades in CGGA-693 cohort. As shown in Fig. 7D, patients in high MLDPS group harbored significantly higher gamma delta T cells, natural killer T cells, T follicular helper cells, neutrophils, activated CD4 T cells and myeloid-derived suppressor cells (MDSCs) than low MLDPS groups in WHO II subgroup. Similar results were observed in patient with WHO III grade (Fig. 7E). Moreover, patients with WHO IV grade in high MLDPS group manifested with obviously elevated memory B cells, immature dendritic cells, activated CD4 T cells, central memory CD4 T cells, gamma delta T cells, which is consistent with the results in WHO II and WHO III subgroups (Fig. 7F). While CD56dim natural killer cells, eosinophils and monocytes were prominently increased in low MLDPS group (Fig. 7F). In addition, the CIBERSORTx analyses revealed that high MLDPS group had obviously fewer resting NK cells than low MLDPS group in WHO II subgroup, while neutrophils and M2 macrophages were evidently higher in high MLDPS group in WHO III subgroup (supplementary Fig. 5). These results depicted a distinctive immune cell infiltration landscape toward high and low MLDPS groups among patients with different grades.

Predictive value of MLDPS in immunotherapy

Given that the high and low MLDPS groups possessed different tumor immune microenvironment, we speculated that there might be differences in immunotherapy response between patients with high and low MLDPS. As shown in Fig. 7G, the Kaplan-Meier curve demonstrated that in PRJNA482620, a cohort that included glioblastoma patients who received anti-PD-1 therapy, patients in low MLDPS group had better prognosis (p = 0.0019). Similarly, advanced urothelial carcinoma patients in low MLDPS group exhibited better prognosis after PD-L1 therapy in the IMvigor 210 cohort (p = 0.012, Fig. 7H). As expected, low MLDPS group of patients received anti-PD-1 therapy had better prognosis in two advanced melanoma cohorts including GSE91061 and GSE78220 (p < 0.05, Fig. 7I-J). Furthermore, the stacked histogram indicated that patients in low MLDPS group were more likely to respond to immunotherapy in PRJNA482620 cohort (Fig. 7K), IMvigor 210 cohort (Fig. 7L), GSE91061 cohort (Fig. 7M) and GSE78220 cohort (Fig. 7N). Overall, the above results implied that MLDPS had the potential to be a prognostic signature for evaluating prognosis in patients receiving anti-PD-1/PD-L1 therapy including glioblastoma.