RNA-seq data (level-3, HTseq-FPKM) and clinical information of lower-grade glioma (LGG, grade II–III) samples were obtained from The Cancer Genome Atlas (TCGA) dataset. A total of 481 samples were finally selected after removing samples with no survival state, no WHO grade, recurrent tumor, reduplicated sequencing, and whose survival time were less than 1 day. In addition, RNA-seq and clinicopathological data were obtained from the Chinese Glioma Genome Atlas (CGGA) website as the validation set. A total of 332 samples with complete survival information were chosen from the CGGA-693 dataset (CGGA-LGG-1) and 162 samples were obtained from the CGGA-325 dataset (CGGA-LGG-2). The batch effect correction was performed by the R package termed “SVA.”

Single-sample Gene Set Enrichment Analysis (ssGSEA) was conducted in the three datasets based on the expression level of 29 immunity-associated signatures by the R package “GSVA.” Consensus clustering was then applied to define the immune subgroups based on the ssGSEA scores by the “ConsensusClusterPlus” package in R. The K-means clustering algorithm was performed with 100 resampling iterations by random selection of 80% of the total samples to ensure the clustering stability. The best cluster number was determined by the consensus matrix (CM) heat maps, cumulative distribution function (CDF) curves, and delta area score of CDF curves. A principal component analysis (PCA) was used to illustrate the reliability of optimal methods. Thorsson et al. (2018) had identified six immune function subtypes by an extensive immunogenomic analysis: wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depleted (C4), immunologically quiet (C5), and TGF-β response (C6). A Sankey plot was applied to visualize the relationships between our four clusters and the six identified immune functional subtypes mentioned previously.

The Estimation of Stromal and Immune cells in Tumors using Expression data (ESTIMATE) algorithm was used to evaluate the LGG microenvironment (Yoshihara et al., 2013). Immune scores and stromal scores were calculated to reveal the abundance of infiltrating immune and stromal cells. ESTIMATE scores were calculated for reflecting non-tumor composites. Tumor purity was inferred by the aforementioned scores. The Kruskal–Wallis test was used to compare differences in multiple clusters. Heat maps were drawn by the “pheatmap” package in R.

The Microenvironment Cell Populations-counter (MCP counter) algorithm (Becht et al., 2016) was used to quantitate the abundance of immunocytes in heterogeneous tissues by the “MCPcounter” R package. The absolute abundances of two stromal cells and eight immune cells were evaluated by immune cell scores, including T cells, CD8 T cells, cytotoxic lymphocytes, B lineage, NK cells, monocytic lineage, myeloid dendritic cells, neutrophils, endothelial cells, and fibroblasts.

The SubMap analysis (Hoshida et al., 2007; Roh and Chen, 2017) from Gene Pattern (https://www.genepattern.org/) was used to predict the response to immune checkpoint blockade (anti-PD-1 and anti-CTLA-4 immunotherapy). In addition, the chemotherapy response was predicted based on the public pharmacogenomics database termed “Genomics of Drug Sensitivity in Cancer” (GDSC, http://www.cancerrxgene.org). The half-maximum inhibitory concentration (IC50) of each patient was estimated by the R package “pRRophetic” with Ridge’s regression, and the accuracy of the prediction was estimated by a 10-fold cross-validation. The IC50 of each sample in TCGA dataset was calculated based on the prediction models of bleomycin and doxorubicin, and cisplatin and gemcitabine.

The GSEA algorithm was used to investigate the biological functions and pathways of clusters A and D, with C2:CP KEGG gene sets from MSigDB as the reference gene sets. False discovery rate (FDR) < 0.05 was the screening threshold.

A weighted gene co-expression networks analysis (WGCNA) algorithm was used to mine the synergistically expressed gene modules. Immune-related genes from the ImmPort dataset (https://www.immport.org/) were classified into different modules which were significantly correlated with the four immune clusters by the R package “WGCNA.” Samples were clustered by a hierarchical clustering algorithm implemented in the R function “hclust.” The soft thresholding power β = 3 was selected by the R function “pickSoftThreshold” (scale free R² = 0.85). The expression matrix was converted into the adjacent matrix and then into the topological matrix for gene clustering. An average linkage hierarchical cluster approach was used to cluster genes into a dendrogram. The STRING database (Szklarczyk et al., 2011) was explored to construct the protein–protein interaction (PPI) network. In the PPI gene network of the target module, hub genes were the top ten genes ranked by the MCC algorithm of “cytoHubba” plugin in Cytoscape 3.8.0. In addition, the survival curves based on the best cut-off value of hub genes were drawn by the “survminer” package in R.

Student’s t-test was applied for normal distributions, and the Mann–Whitney U-test was performed for nonparametric distribution. Chi-square or Fisher’s exact tests were used for categorical data. Kaplan–Meier curves and log-rank tests were used to evaluate the survival time of different immune clusters. The nonparametric Kruskal–Wallis (KW) test was used to analyze the difference in IC50 in different clusters. The Benjamini–Hochberg procedure was applied to control the false discovery rate (FDR) for multiple testing. p < 0.05 was considered statistically significant (* represented p < 0.05, ** referred p < 0.01, and *** referred p < 0.001).

Unsupervised consensus clustering was applied to explore a novel immune classification of LGGs based on the ssGSEA scores of TCGA dataset. The optimal clusters number was found to be four with maximal consensus within clusters and minimal ambiguity among clusters (Figures 1A–C). PCA verified that the ssGSEA scores could completely be distinguished into four subtypes which were referred to as cluster A, cluster B, cluster C, and cluster D in TCGA dataset (Figure 1D). The clustering results were the same in CGGA-1 (Supplementary Figure S1A) and CGGA-2 datasets (Supplementary Figure S1B). The Sankey diagram revealed the immune function characteristics of the four clusters (Figure 1E). The majority of clusters A and B were related to the C5 function of “Immunologically Quiet.” Cluster D was related to the C4 function of “Lymphocyte Depleted.”

We identified four immune clusters and their immune characteristics in TCGA, CGGA-1, and CGGA-2 datasets (Figures 2A, 3A, 4A) were shown in the heat maps. Cluster D showed the highest degree of immune infiltration and it was followed by clusters C, B, and A. Cluster D was considered the “immune-rich” phenotype with the highest enrichment scores and the least tumor purity while cluster A was the opposite which was regarded as an “immune cold” phenotype (Figures 2A–D, Figures 3A–D, Figures 4A–D). Apart from that, clusters A and B could be seen as a “low-immune infiltration” subgroup whereas clusters C and D were considered as a “high-immune infiltration” subgroup. The expression of the immune checkpoint genes PDCD1 (PD-1), CD274 (PD-L1), PDCD1LG2 (PD-L-2), CTLA-4, HAVCR2, and LAG3 which played a vital role in the oncogenesis and progression of LGG were expressed in the following order: D > C > B > A (Figures 2E–J, Figures 3E–J, Figures 4E–J).

To evaluate the clinical features and gene mutation characteristics among the four immune clusters, age, gender, tumor grade, IDH1 (R132) status, IDH2 R172 status, PTEN status, EGFR status, ATRX status ,and TP53 status in TCGA dataset were counted (Table 1; Figure 5). In clusters A, C, and D, people aged more than 40 years accounted for more proportion, whereas those younger than 40 years were more common in cluster B. WHO grade II glioma tended to be common in the “low-immune infiltration” subgroup (clusters A and B), whereas the “high-immune infiltration” subgroup (clusters C and D) counted more in WHO grade III glioma. The frequency of IDH1 (R132) mutation was much higher in the “low-immune infiltration” subgroup than the “high-immune infiltration” subgroup. The frequency of PTEN and EGFR mutations was significantly higher in cluster D, which had the highest immune infiltration. In clusters B and C, which had mild immune infiltration changes, ATRX mutation frequencies were higher than those in clusters with extreme immune infiltration changes (clusters A and D). TP53 mutation was common in cluster B. Gender and IDH2 (R172) mutation status were not significant covariates in the immune classification. In addition, the “low-immune infiltration” subgroup (clusters A and B) showed longer overall survival than the “high-immune infiltration” subgroup (clusters C and D) (Figure 6), indicating that immune infiltration of LGG played a negative role in the prognosis.

To explore the differential distribution of immunocytes and stromal cells in tumor immunity clusters, the MCPcounter algorithm was used to calculate the contents of two stromal cells and eight immune cells in the four clusters in TCGA, CGGA-1, and CGGA-2 datasets (Figures 7A–C). Immune cell scores of CD8 T cells, B lineage, NK cells, myeloid dendritic cells, endothelial cells, and fibroblasts in the “high-immune infiltration” subgroup (clusters C and D) were significantly higher than those in the “low-immune infiltration” subgroup (clusters A and B). Then, the correlation landscape of immunocytes was characterized to compare the relative subpopulations of infiltrating immune cells and immune scores among the four cluster patterns (Figure 7D). Cox regression analysis of the 10 immune cells in TCGA, CGGA-1, and CGGA-2 datasets are shown in Supplementary Figure S2, revealing the prognostic risk factors of infiltrating immunocytes.

We compared the expression profiles of the four immune clusters in TCGA datasets by the Subclass Mapping algorithm which assessed the response to anti-PD-1 and anti-CTLA-4 therapies. A significant correlation was observed when comparing cluster D with the PD-1 response group (Bonferroni-corrected p = 0.001, Supplementary Figure S3A). It revealed that cluster D might have a better response to anti-PD-1 therapy while no significant correlation of anti-CTLA-4 therapy was observed in all the clusters. The “pRRophetic” algorithm was applied to evaluate the sensitivity of four common chemical drugs: cisplatin, bleomycin, doxorubicin, and gemcitabine for the four immune clusters. A lower IC50 value would indicate a better sensitivity in the prediction models. For bleomycin and doxorubicin, the “low-immune infiltration” subgroup (clusters A and B) was more sensitive than the “high-immune infiltration” subgroup (clusters C and D). For cisplatin and gemcitabine, cluster D was the most sensitive and cluster A was the least sensitive compared with the other clusters (Figures 8A–D). Moreover, to compare the accuracy of the four immune clusters, prognosis signatures in other references were used to compare the C-index. The results were also exciting: in the C-index for predicting the LGG survival possibility, our immune clusters showed better predictive value than other signatures (Maimaiti, Aierpati et al., 2022; Maimaiti, Aierpati et al., 2021) (0.813 > 0.774 > 0.712 > 0.662, Supplementary Figure S3B).

A GSEA analysis was performed to screen the correlated biological pathways in immune clusters A and D. Cluster D was enriched in Allograft rejection, complement and coagulation cascades, cytokine–cytokine receptor interaction, graft versus host disease, antigen processing and presentation, cell adhesion molecules (CAMs), ECM-receptor interaction, and focal adhesion in TCGA, CGGA-1, and CGGA-2 datasets. Enrichment results of cluster A were not significant under the strict FDR <0.05 threshold in all the three datasets (Figure 9).

Tumor mutational burden of coding errors is reported to have a certain correlation with the tumor immune microenvironment. We explored this correlation of different immune clusters in TCGA-LGG datasets (Figures 10A–D). The frequency of IDH1 missense mutations in the “low-immune infiltration” subgroup was higher than that in the “high-immune infiltration” subgroup (80 and 87% in clusters A and B, 61 and 43% in clusters C and D). TP53 mutations were higher in clusters B (52%) and C (49%) than those in clusters A (29%) and D (36%). Meanwhile, most of them were missense mutations. The CIC missense mutation was high in cluster A. ATRX mutations including missense mutations, nonsense mutations, and multi-hit were at a high frequency in clusters B and C. TTN, EGFR, and ATRX mutations were common in the “high-immune infiltration” subgroup. PTEN, KEL, and PLK3CA mutations were higher in cluster D. the TP53 pathway was highly affected in cluster A and RTK-RAS and TP53 pathways were affected in cluster D (Figures 10E,F).

The WGCNA networks of immuno-related genes with immune infiltrating clusters were constructed. The optimal soft thresholding power β was selected (Figures 11A, Figure 9) and three modules were obtained (Figures 11B,C). Red, blue, and green modules were positively correlated with the “high-immune infiltration” subgroup (clusters C and D) and negatively correlated with the “low-immune infiltration” subgroup (clusters A and B), p < 0.05 was used as the threshold. The turquoise module was positively correlated with clusters C and D and negatively correlated with cluster A. The brown module was positively correlated with cluster A and negatively correlated with clusters B and D. Pink and yellow modules were positively correlated with cluster B and negatively correlated with cluster A. In the blue module, which was most correlated with cluster D, the top ten hub genes selected by the MCC algorithm in the PPI network were CD28, CD8A, CSF2, GZMB, IFNG, IL15, IL2, IL2RA, IL7R, and PRF1(Figure 11D). Hub genes were positively correlated with most of the immune cells and immune functions, such as HLA, CCR, etc., (Figure 11E). A survival analysis showed that high expression of CD28, CD8A, IFNG, IL2RA, IL7R, IL15, and PRF1 had a poor prognosis whereas a better prognosis was found in IL2 and GZMB (Figure 12).