The available transcriptome and clinical data of GBM patients from the Chinese Glioma Genome Atlas (CGGA), Gene Expression Omnibus (GEO), and The Cancer Genome Atlas (TCGA) were used as the datasets for finding prognostic biomarkers and constructing prediction models. The samples, whose histology grade was grade IV or diagnosed with GBM, were downloaded and kept for analysis. There were 237 GBM samples from the CGGA dataset (Zhao et al., 2021), 150 GBM samples from GSE16011 (Gravendeel et al., 2009), and 159 GBM samples from the TCGA-GBM project. Another independent GBM dataset (GSE83300) that contained 50 GBM patients was used as the validation dataset for evaluating the performance of prediction models (Feng et al., 2017). There was a median survival time of 12, 8.7, 9, and 16.8 months in the CGGA, GSE16011, TCGA-GBM, and GSE83300 datasets, respectively. All the accession IDs of used samples from the four datasets are displayed in the Supplementary Table S1.

The expression values of genes in each dataset were normalized by a min–max normalization method using the following formula: z_i = (x_i–min(x))/(max(x)—min(x)). Here, z_i is the normalized value of the gene expression, x_i is the expression value before normalization, min(x) is the minimum value of the gene expression in the dataset, and max(x) is the maximum value of the gene expression in the dataset. After normalization, the gene values range from 0 to 1. We used a univariate Cox regression analysis in CGGA, GSE16011, and TCGA-GBM to examine the links between ECM gene (1,936 unique ECM-related genes from 47 ECM gene sets) expression and overall survival (OS). The prognostic ECM genes were defined as the genes with p-value<0.05 in the univariate Cox regression analysis. The prognostic ECM genes were overlapped among the three datasets. Then, CGGA was set as the training dataset, GSE16011 and TCGA-GBM were set as the internal testing dataset, and GSE83300 was set as the independent validation dataset. The LASSO method from the “glmnet” package was adopted in CGGA to select the prognostic genes with a high importance value (Friedman et al., 2010). The ECM index was calculated by applying a multivariate Cox regression analysis to calculate the regression coefficient of genes.

The ECM index was calculated using the following formula: ECM index = β1* Gene_1 +…+ βn*Gene_n. In the formula, β is the regression coefficient and Gene_n is the expression value of genes. In the training, testing, and validation datasets, the patients were classified as high or low index based on the median value of the ECM index. Subsequently, we calculated the log-rank test and plotted the Kaplan–Meier curve to evaluate the difference in OS between the high– and low–ECM index groups. Then, the OS prediction ability of the ECM index was assessed by the area under the curve (AUC) at 12, 18, 24, and 36 months.

The expression values of the selected genes were compared between the GBM and normal brain samples from the GEPIA web server (Tang et al., 2017). The GEPIA web server contains 163 GBM from the TCGA-GBM project and 207 normal brain samples from the GTEx project. For comparison of LGG and GBM, we downloaded and used the expression profiles of 159 GBM samples and 500 LGG samples from the TCGA-GBM and TCGA-LGG projects, respectively.

The MCP-counter which is capable of estimating the absolute abundance of eight immune and two stromal cell populations in the tumor samples by the expression data was used in this study (Becht et al., 2016). The differences in the immune and stromal cell populations were compared between the high– and low–ECM index groups using t-test analysis.

In CGGA, GSE16011, and TCGA-GBM, the differences in the expression levels of protein coding genes between the high–ECM index group samples and low–ECM index group samples for each ECM gene were analyzed using the linear models from the “limma” R package (Ritchie et al., 2015). We defined the cutoff criteria as p-value<0.05 and log2 fold change>0.5 to filter the statistically significant differentially expressed genes (DEGs) in each dataset. The DEGs were overlapped among the three datasets to obtain the robust DEGs. By using the Connectivity Map (CMap) database, we were able to examine the potential drugs with a close correlation to the diseases. The robust DEGs between the low– and high–ECM index groups were submitted to the Cmap database. The connectivity value from the Cmap represents the ability of the drug to reduce ECM formation, and the optimal connectivity value should be −1.

Machine learning methods were used in the construction of models to predict the survival rates of GBM patients at 6, 12, 18, 24, 30, and 36 months after diagnosis. Four machine learning algorithms, including support vector machine (SVM), random forest (RF), generalized linear models (GLM), and artificial neural network (ANN) from the caret R package were considered. To identify the best parameters for model creation, five-fold cross-validation and grid search were utilized. Data merging was performed on the CGGA, GSE16011, and TCGA-GBM datasets, and the samples were separated into the training dataset (60%) and the testing dataset (40%). GSE83300 was set as the independent validation dataset. After setting the best parameters, the models were trained on the training dataset and evaluated on the testing dataset and an independent validation dataset. The receiver operating characteristic (ROC) curve, a plot of sensitivity versus (1–specificity), was used to estimate the performance of models. We used the R package “shiny” to develop a web-based survival rate calculator and survival curve predictor using these machine learning models.

In this study, 47 ECM gene sets containing 1,936 unique ECM-related genes were obtained (Supplementary Table S2). Among these 1,936 ECM-related genes, by the univariate Cox regression analyses, 134 genes were found to be positively correlated with prognosis in CGGA, 118 genes were found to be positively correlated with prognosis in GSE16011, and 74 genes were found to be positively correlated with prognosis in TCGA-GBM. These genes were identified as protective genes. Similarly, 223 genes were negatively correlated with prognosis in the CGGA dataset, 205 genes were negatively correlated with prognosis in the GSE16011 dataset, and 211 genes were negatively correlated with prognosis in the TCGA-GBM dataset. These genes were selected as risky genes. The intersection of the results showed that there were two common protective genes and 37 common risky genes (Figures 1A,B).

The LASSO technique was used on these 39 prognostic ECM genes in the CGGA, and the results indicated that AEBP1, F3, FLNC, IGFBP2, and LDHA were the best combination of genes to construct the model (Figures 1C,D). A multivariate Cox regression analysis calculated the regression coefficient. The ECM index was obtained by the following formula: [1.0391×AEBP1] + [0.6346×F3] + [0.5396×FLNC] + [2.1276×IGFBP2] + [2.7396×LDHA]. All the five selected genes (AEBP1, F3, FLNC, IGFBP2, and LDHA) were the risky genes since their regression coefficients were positive. The median value of the ECM index was selected to divide the patients into high– and–low ECM index groups. In Figure 2A, the ECM index distribution and overall survival data in the CGGA dataset were displayed and ranked according to the index value. The survival analysis indicated that the low–ECM index group had an apparent better prognosis than the high–ECM index group (p-value<0.05; Figure 2B). In Figure 2C, gene expression profiles of the high and low index groups were illustrated. The AUC values for 12-, 18-, 24-, and 36-month OS predictions were 0.63, 0.67, 0.69, and 0.66, respectively, which indicate its high ability to predict the prognosis of GBM (Figure 2D).

In CGGA, GSE16011, and TCGA-GBM, the GBM patients were classified into two groups by the median gene expression of the five selected genes (AEBP1, F3, FLNC, IGFBP2, and LDHA). Then, survival analysis was used to assess the relationship between the gene expression value and OS. There were significant negative correlations between the expression levels of AEBP1, F3, FLNC, and LDHA in CGGA (p-value<0.05; Supplementary Figure S1). The findings were then validated in the GSE16011 and TCGA-GBM datasets. In GSE16011, the gene expression levels of AEBP1, FLNC, and LDHA were negatively correlated with OS (Supplementary Figure S2). In TCGA-GBM, there were significant negative correlations between the expression levels of AEBP1 and FLNC with OS (p-value<0.05; Supplementary Figure S3).

The expression values of AEBP1, F3, FLNC, IGFBP2, and LDHA were compared among normal brain tissue, low-grade gliomas (LGG), and GBM. The results from the GEPIA dataset showed that the expression values of these genes were higher in the GBM samples than in the normal brain samples (p-value<0.0001; Supplementary Figure S4). The expression values of these genes were also significantly higher in the GBM samples than in the LGG samples (p-value<0.0001; Supplementary Figure S5).

To further evaluate the robustness of the ECM index, the ECM index of the samples in the GSE16011 and TCGA-GBM datasets was obtained using the same formula. The median value of the ECM index was also used to stratify the samples into the high– and low–ECM index groups. The ECM index distribution, overall survival, and ECM gene expression profiles in these two validation sets are plotted in Figure 3A, D. The ROC analysis results in the GSE16011 and TCGA-GBM datasets suggested that the ECM index had a powerful ability to predict GBM survival (Figures 3B,E). Compared to the high–ECM index patients, the low–ECM index patients had better overall survival prognosis (GSE16011: p-value<0.0001; TCGA-GBM: p-value = 0.0048; Figures 3C,F).

Another independent dataset (GSE83300), which contains 50 GBM patients, was used to validate the results. The ECM index distribution, overall survival, and ECM gene expression profiles in GSE83300 are plotted in Figure 4A. The ROC analysis results suggested that the ECM index had a powerful ability to predict GBM survival (Figure 4B). Compared to high–ECM index patients, the low–ECM index patients had better overall survival prognosis (p-value = 0.046; Figure 4C), which is in concordance with the results from the CGGA, GSE16011, and TCGA-GBM datasets.

We calculated the association between the ECM index and population abundances of the immune and stromal cells in the combined CGGA, GSE16011, and TCGA-GBM datasets. As shown in Figure 4D, the high index group samples were found to comprise higher endothelial cells and fibroblast.

Data profiles from CGGA, GSE16011, and TCGA-GBM were used to find the DEGs between the high– and low–ECM index groups. DEGs calculated from the three datasets were overlapped to obtain the robust DEGs, and a total of 286 robust DEGs were identified. To predict the potential drugs for GBM, CMap was implemented on the robust DEGs between the low– and high–ECM index groups. The 20 drugs with high and significant associations with robust DEGs are listed in Supplementary Table S3. Among these drugs, podophyllotoxin, lasalocid, MG-262, and nystatin revealed higher negative correlations and the potential to reverse the high ECM index tumor status.

The combination of genes, age, and gender was used to predict the survival status of GBM patients. The survival status of 6, 12, 18, 24, and 30 months was predicted by the machine learning methods. First, CGGA, GSE16011, and TCGA-GBM were combined into one dataset. Then, the combined dataset was randomly divided into training (60%) and testing (40%) datasets. The commonly used machine learning classifiers, including GLM, ANN, SVM, KNN, and RF, were constructed to predict the survival status of patients at 6, 12, 18, 24, 30, and 36 months after the treatment, respectively. After setting the best parameter via a five-fold cross-validation strategy, we calculated the AUC value in the testing dataset to characterize the ability of the model to distinguish between the dead and alive. KNN, RF, ANN, ANN, RF, and RF demonstrated the highest AUC values in predicting the survival status at 6, 12, 18, 24, 30, and 36 months (Figures 5A–E, F), respectively. Then, these constructed machine learning models were also validated by an independent dataset (GSE83300). RF, SVM, GLM, GLM, ANN, and GLM demonstrated the highest AUC values in predicting the survival status at 6, 12, 18, 24, 30, and 36 months (Figures 6A–F), respectively.

To make these constructed machine learning models accessible to GBM researchers, we developed a web tool to predict the survival probability of GBM patients called OSPG. OSPG can be used in five steps 1) visiting the website https://ospg.shinyapps.io/OSPG/;(2) inputting the values of five genes, including AEBP1, F3, FLNC, IGFBP2, and LDHA (gene expression values range 0–1); 3) inputting the values of age and gender (male or female); 4) selecting the machine learning models for each time point (6, 12, 18, 24, 30, and 36 months), and 5) clicking the “submit” button. An example is illustrated in Figure 7.