Gene transcriptome data for low-grade gliomas (LGGs) and glioblastomas (GBMs) were collected from The Cancer Genome Atlas (TCGA), GSE108474 was downloaded as background corrected and RMA-normalized data, and the IMvigor 210 cohort was used as an immunotherapy cohort to verify clinical treatment outcomes. Among them, TCGA-LGG included a total of 515 samples’ transcriptome and clinical data, and TCGA-GBM included a total of 599 samples’ transcriptome and clinical data.

All transcriptome profiles were converted to transcripts per million (TPM) prior to analysis. Transcriptome data obtained from different platforms were corrected using the normalizeBetweenArrays function. Finally, transcriptome matrix data were screened to calculate the average value of duplicated genes.

CIBERSORTx is a suite of machine learning tools for assessing cell abundance and cell-type-specific gene expression patterns from a large set of tissue transcriptome profiles (Steen et al., 2020). QuanTIseq is a method to quantify the tumor immune background as determined by the type and density of tumor-infiltrating immune cells (Finotello et al., 2017). TIMER explores the association between immune infiltration and a wide range of factors by performing predictive calculations of the levels of tumor-infiltrating immune subsets of cancer types with a comprehensive study of the molecular characteristics of tumor-immune interactions (Li et al., 2017). CIBERSORTx, QuanTiseq, and TIMER were employed to analyze the cell infiltration in TCGA samples, and the immune cell infiltration from TCGA data was investigated to identify the relationship between infiltrated immune cells, especially CD4⁺ and CD8⁺ T cells, and overall survival (OS) of glioma patients, with p < 0.05 being the entry for further analysis.

The Gene Set Enrichment Analysis (GSEA) algorithm was applied to perform biological pathway enrichment between the two groups with Hallmark, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) reference gene sets.

Weighted correlation network analysis (WGCNA) can be used to describe the correlation between candidate features and gene sets. WGCNA was conducted using the R package WGCNA, and to ensure scale-free topological networks, the power value was set to β = 3 and the soft threshold parameter of scale-free networks to R2 = 0.9. Seven modules were retrieved as a result, with the green module displaying the closest relationship to the clinic for further study.

Random forest analysis was first performed on the selected genes, with the initial nTree = 1,000. When the error rate reached the minimum at nTree = 859, seven genes with relative importance greater than 0.2 were extracted for analysis. LASSO regression analysis was applied to the model equation, as follows: r i s k   s c o r e = C o e f 1 × G e n e   e x p r e s s i o n 1 ＋ C o e f 2 × G e n e   e x p r e s s i o n 2 ＋ ⋯ C o e f n × G e n e   e x p r e s s i o n n

Coef represents the importance index for each gene in the analysis random forest analysis. Gene expression values represent the expression values of the corresponding genes. For the assigned risk score, patients above the median risk score were defined as the high-risk group, whereas those below the median risk score were defined as the low-risk group. ROC and KM curves were used to assess model performance.

GBM cell lines (U251) and normal human astrocyte cell lines (HEB) were obtained from Dr. Cai and Dr. Chen, respectively. The cells were seeded in RPMI-1640/DMEM supplemented with 10% fetal bovine serum (Gibco, China) at 37°C in a 5% CO₂ atmosphere.

Cells were treated with TRIzol reagent (Takara, Japan). We then extracted all RNA and reverse-transcribed it into cDNA. qRT-PCR was used to analyze the relative expression of HSPB1, and data were normalized to GAPDH. The primers used are listed in Supplementary Table S1.

All data analysis results in our investigation were processed using R software. For continuous variables with normal distribution, we conducted Student’s t-test, while non-normally distributed continuous variables were calculated using the Mann–Whitney U test. Differential expression analysis was conducted with a threshold set at |logFC| > 0.5 and p < 0.05.

First, the three pathways of IL2, IFN-γ, and TNFα were assessed for their Gene Set Variation Analysis (GSVA) scores according to the HALLMARK gene set, based on which unsupervised clustering analysis was conducted. The delta results showed that grouping into four clusters yielded the best results (Figure 1A). The heat map illustrates the expression levels of IL2, IFN-γ, and TNFα pathways in different subtypes, and we found that the pathways of IL2, IFN-γ, and TNFα all presented a state of high expression in cluster B (Figure 1B). Subsequently, the proportion of glioma and LGG patients in different TEXcluters was examined, with GMB patients accounting for the highest proportion (50%) in Cluster B and only 1% in Cluster D (Figure 1C). Meanwhile, the predictive efficacy of their TEX groupings was tested, all of which showed that patients in Cluster D had the best survival status, whereas patients in Cluster B had the worst survival (Figure 1D). The immune cell infiltration component of the tumor microenvironment of the four cluster patients was assessed using the CIBERSORTX algorithm, and Cox regression analysis was performed to investigate the prognostic value of a wide range of cells in each TEXCluster patient. Among them, resting T cell memory was found to be a protective factor in both TEXb and TEXc (Figure 1E). Therefore, resting T cell memory was grouped according to its level of infiltration and tested to predict clinical outcomes in glioma patients. Intriguingly, patients with high T cell memory resting infiltration had poorer survival prognosis (Figure 1F). Therefore, further analysis was conducted.

Differential expression analysis was performed on patients with TEXCluster 4, and concatenation was performed for WGCNA analysis. The WGCNA method was used to screen for clinically relevant genes, and a weighted gene co-expression network was constructed using the following parameters: power of β = 4 and scale-free R2 = 0.9. Consequently, seven color modules were finally obtained by merging similar modules, and the MEgreen group was found to be the most relevant to survival-related information, after which MEgreen was selected for subsequent analysis (Figures 2A,B).

Based on the selected gene sets and corresponding clinical outcomes, we performed survival random forest analysis, an algorithm was used to process right-censored survival data. When Ntree = 859, the error rate was the lowest, and the seven genes most associated with clinical outcomes were HSPB1, HOXD10, HOXA5, SEC61G, H19, ANXA2P2, and HOXC10 (Figure 2C). The importance of these genes in predicting the clinical outcome of patients with glioma is shown in Figure 2D. The LASSO regression algorithm was used to construct the model. r i s k   s c o r e =   ″ 0.4956 × H S P B 1 + 0.1231 × H O X D 10 + 0.0697 × H O X A 5 + 0.0037 × SEC ⁡ 61 G + 0.0911 × H 19 + 0.145 × A N X A 2 P 2 + 0.064 × H O X C 10

As shown by the Kaplan–Meier curve, patients with a low TEXScore had a significantly higher OS rate than those with a high TEXScore (p < 0.001), suggesting that the overall survival time of patients would decrease with increasing risk score (Figure 2E), and the ROC curve indicated that its performance in predicting survival rate at 1, 3, and 5 years was 0.870, 0.887, and 0.843, respectively (Figure 2F). It was found by the multivariate COX regression analysis that TEXSore was an independent factor for prognosis with a Hazard ratio of 2.014 (Figure 3A). An external dataset was used to validate the predictive performance of the TEXScore model, and the Kaplan–Meier survival results were in agreement with the training set (Figure 3B). We then detected the expression of HSPB1, the most important gene, in the normal human astrocyte cell line (HEB) and GBM cell line (U251), and found significant differences. HSPB1 was expressed at higher levels in U251 cells than in HEB cells (Supplementary Figure S1). These results indicated that HSPB1 plays an important role in glioma progression.

Based on the HALMARK, KEGG, and GO datasets, we explored the correlation between TEXScore and biological function in patients with glioma. First, we found that the TEXScore was correlated with tumor pathways, including positive regulation of cell activation, focal adhesion, and the JAK-STAT signaling pathway. The TEXScore was also highly correlated with toll-like receptor signaling, granulocyte migration pathway, adaptive immune response, and other immune pathways (Figures 3C–E). Both the HALMARK and KEGG gene sets confirmed these results (Figure 4).

The levels of pathway enrichment in patients with high and low TEXScore are illustrated in the heatmap (Figure 5A). Based on the above GSEA results, several gene sets were selected for GSVA analysis, and a negative correlation was found between TEXScore and tumor formation and other pathways, such as the ERBB signaling pathway and WNT signaling pathway, but a positive correlation with mismatch repair, antigen processing, and presentation (Figure 5B).

To check the effectiveness of the TEXScore in predicting the clinical response to immunotherapy in glioma patients, the expression of immune checkpoints in patients with high and low TEXScore was examined. The results showed that immune checkpoint expression was higher in patients with a high TEXScore than in those with a low TEXScore (Figure 5C). Subsequently, the degree of immune cell infiltration in patients in the high and low TEXScore groups was calculated using quanTlseq and TIMER algorithms. The quanTlseq algorithm revealed that the M2 macrophage component was significantly enriched in the high TEXScore patients, whereas immune effector cells were predominant in the low-risk group (Figure 6A).

The TIMER algorithm showed a higher proportion of both immune effector cells and M macrophages in the high-risk group, which sheds light on the subsequent analysis of whether the immune effector cells in the high-risk group were in a state of immune depletion (Figure 6B). The mutations in LGG and glioma are shown in Figure 6C. Figures 6D,E present the mutation landscape of the high and low TEXScore patients, respectively, with the high TEXScore group having a lower mutation rate than the low TEXScore group. The overall mutation landscape is illustrated in Figure 6F, including variant classification, mainly composed of missensemutation, variant type, mainly composed of SNPs, SNV class, mainly composed of C > T, and variants per sample, Mediant:32.

Data from IMvigor210 were then obtained to validate the performance of the TEXScore model in predicting the effect of immunotherapy. First, the patients were categorized into high- and low-risk groups according to the TEXScore formula. The Kaplan–Meier survival model proved the effectiveness of TEXScore in predicting the clinical outcome of patients, with patients with a high TEXScore experiencing worse survival status than those with a low TEXScore (Figure 7A). Following this, the bar chart demonstrates that patients with low TEXScore responded more significantly to immunotherapy (Figure 7B), suggesting that the lower the TEXScore, the better the patients responded to immunotherapy, with PR and CR patients having the lowest TEXScore (Figures 7C,D).