The datasets for TCGA GBMLGG, Chinese Glioma Genome Atlas (CGGA) #325, and CGGA #693 cohorts were obtained from the University of California Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/datapages/) (18) and the CGGA data portal (http://www.cgga.org.cn/) (19). After excluding cases with missing clinical data, the study finally included 647, 309, and 656 patients, respectively (Supplementary Tables S1–S3). Three GBM and three low-grade glioma (LGG) tissue samples were collected in the Zhengzhou Central Hospital Affiliated to Zhengzhou University. Clinical information for all tissue samples used in the study was listed in the supplementary material (Supplementary Table S4). The study was approved by the ethics committee of Zhengzhou Central Hospital Affiliated to Zhengzhou University, and informed consent was obtained from all patients.

Abundance of tumor-infiltrating immune and stromal cells in gliomas was calculated by the EPIC method on the TIMER2 platform (http://timer.cistrome.org/) (20). The differentially expressed genes (DEGs) (p < 0.05 and |log₂FC| ≥ 1) were screened using the R package “limma” in TCGA cohort, CGGA #325 cohort, and CGGA #693 cohort, respectively (21). Gene Ontology (GO) analysis of DEGs was performed using R package “clusterProfiler”. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was conducted using the KOBAS-i (http://bioinfo.org/kobas) (22).

DEGs were analyzed using the R software package “WGCNA,” and genes with high connectivity in the tumor-infiltrating fibroblast significant module were identified as hub genes. Univariate, LASSO, and multivariate regression analyses were then performed to screen for positive hub genes significantly associated with overall prognosis. The risk score was calculated as follows:

After baking, dewaxing, rehydration, antigen retrieval, and blocking, paraffin sections were incubated overnight with anti-PLAUR primary antibody (Abcam, ab218106, 1:250) at 4°C. After washing three times, the sections were incubated with biotin-labeled secondary antibody for 20 min at room temperature, followed by avidin-labeled Horseradish peroxidase (HRP) for 20 min at room temperature, and then stained with diaminobenzidine. Finally, the sections were counterstained with hematoxylin.

One-way ANOVA, Wilcoxon test, and t test were used to analyze the significance of differences in gene expression and immune cell infiltration. Univariate regression, LASSO regression, multivariate regression, and Kaplan–Meier analyses were performed using the R packages “glmnet” and “survival”. sinlge Receiver operating characteristic (ROC) curve was drawn using the R package “pROC”. All statistical analyses were performed using GraphPad Prism 8, R software, and SPSS, and p values <0.05 were considered statistically significant.

Previously, seven methods including CIBERSORT, CIBERSORT abs, quanTIseq, MCP-counter, TIMER, xCell, and EPIC have been evaluated for their accuracy in estimating different immune and stromal cells based on gene expression data. Considering the robust overall performance of the EPIC method in predicting TAF infiltration (23), we analyzed RNA Sequencing (RNA-seq) data from the CGGA #325, CGGA #693, and TCGA GBMLGG cohorts by this method to characterize TAFs in glioma (Figures 1, 2). There were significant differences in TAF infiltration in different WHO grades of gliomas, with higher TAFs in high-grade gliomas (Figure 2A). Similarly, the infiltration of TAFs was significantly increased in GBM relative to LGG (Figure 2B). In addition, fibroblast infiltration was significantly associated with the prognosis of gliomas, with high infiltration predicting poorer overall survival (Figure 2C).

The DEGs (|log2FC| ≥ 1 and adjusted p < 0.05) between LGG and GBM were identified using R package “limma”. In this study, 3,868, 9,555, and 1,831 DEGs were screened out in TCGA GBMLGG cohort, CGGA #325 cohort, and CGGA #693 cohort, respectively, of which 1,159 DEGs were identifiable in all cohorts (Figures 3A, B). GO analysis of 1,159 target DEGs was performed using R package “clusterProfiler” (Figure 3C), and then KEGG analysis was conducted using the KOBAS-i (Figure 3D). GO analysis showed that the DEGs between LGG and GBM were mainly involved in the organization and remodeling of ECM (Figure 3C). KEGG analysis showed that the DEGs were significantly enriched in the ECM–receptor interaction pathway, focal adhesion pathway, Phosphatidylinositol 3 kinase-protein kinase B (PI3K-Ak) signaling pathway, Ras signaling pathway, cytokine–cytokine receptor interaction, and chemokine signaling pathway (Figure 3D). Therefore, ECM remodeling may be involved in glioma progression.

We performed WGCNA to identify gene modules that are highly synergistically varied and to identify candidate hub genes based on the interconnectivity of gene modules and the association between gene modules and phenotypes (Figure 4). Four modules were identified from the co-expression network (Figure 4B), among which the red module was most associated with TAF infiltration (Figures 4C, D). Based on the cutoff criteria (|MM| > 0.9 and |GS| > 0.1), 10 genes with high connectivity in the clinically significant module were identified as hub genes. Univariate Cox regression analysis indicated that all the 10 hub genes were significantly associated with prognosis of glioma in TCGA GBMLGG cohort (Figure 5A). Subsequently, we performed LASSO and multivariate Cox regression analyses to analyze the 10 hub genes in TCGA GBMLGG cohort (Figures 5B, C). Finally, S100A4, PLAUR, and EMP3 were identified as prognostic-related hub genes (Figure 5D) and were selected to construct the risk score signature in TCGA GBMLGG cohort (Figure 6A), CGGA #325 cohort (Figure 6D), and CGGA #693 cohort (Figure 6G), respectively. The survival of patients was analyzed using the R package “survival,” and finally, we observed a significant prognostic difference between the high-risk and low-risk groups (Figures 6B, E, H). We performed ROC analysis using the R package “pROC” at 1-, 3-, and 5-year time points to assess the sensitivity and specificity of risk score in predicting the survival of glioma patients and found that the predictive accuracy of the risk score was very high (Figures 6C, F, I).

We then analyzed the relationship between the risk score and clinical phenotypes. A Sankey diagram was created to display the distribution of the survival status, WHO grade, Isocitrate dehydrogenase (IDH) status, risk score, and TAF infiltration of glioma patients (Figure 7A). As shown in Figure 7B, the risk score of glioma in GBM was significantly higher than that of the corresponding LGG subtype, and the IDH wild-type group had a higher risk score relative to the IDH mutant group (Figure 7C). In addition, the risk score signature also showed high prognostic value across different WHO grades and IDH status subtypes, and significant prognostic differences between high- and low-risk groups were observed in all subtypes (Figures 7D–G).

As the CGGA #693 cohort contained many recurrent gliomas (Supplementary Table S3), we analyzed the expression of prognostic genes and the risk scores in primary and recurrent subgroups, respectively. The expression of S100A4, EMP3, and PLAUR was increased in recurrent gliomas compared with primary gliomas (Supplementary Figures S1A–C) as were the risk scores of recurrent gliomas (Supplementary Figure S1D). Furthermore, recurrent gliomas had worse overall survival than primary gliomas (Supplementary Figure S1E), and the risk score signature also showed a high prognostic value in both primary and recurrent subgroups (Supplementary Figures S1F, G).

We analyzed the relationship between risk score and TAF infiltration and found that the high-risk group had increased levels of fibroblast infiltration in all three cohorts (Figure 8A). In addition, there is a significant correlation between the risk score and TAF infiltration (Figure 8B). Prognostic-associated genes (S100A4, PLAUR, and EMP3) that constructed the risk signature were also significantly associated with TAF infiltration in gliomas (Figure 8B).

To investigate whether risk score is an independent prognostic factor for gliomas in TCGA GBMLGG cohort, we performed univariate and multivariate Cox regression analyses, respectively. The results showed that the risk score was significantly correlated with prognosis and was identified as an independent prognostic factor for glioma (Figures 9A, B). We then built a survival nomogram prediction model based on independent prognostic parameters for gliomas in TCGA GBMLGG cohort (Figure 9C). The calibration curves displayed excellent agreement between observation and prediction at 1-, 3-, and 5- year time points (Figure 9D).

To validate the expression of prognostic-associated genes that constructed the risk signature in glioma, we collected tissue samples resected from patients undergoing surgical treatment. Similar to the mRNA transcriptome results of TCGA GBMLGG cohort, CGGA #325 cohort, and CGGA #693 cohort, PLAUR expression was increased in GBM relative to LGG (Figure 10). Probably due to the limitation of the sample size collected in this study, only one gene was validated by immunohistochemical staining. Therefore, more samples may be required to verify the expression of S100A4 and EMP3.