All transcriptomic and clinical characteristics of enrolled samples were extracted from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases (Brat et al., 2015; Zhao et al., 2021). A total of 952 primary LGG samples with detailed clinical information were enrolled in our study, in which 508 samples extracted from the TCGA database were defined as the training set; 444 samples extracted from the CGGA database were defined as the validation set. Normal or glioblastoma (GBM, grade IV glioma) samples were excluded.

CIBERSORT algorithm was employed for evaluating the percentage of 21 human hematopoietic cell phenotypes, including seven T cell types, naïve and memory B cells, plasma cells, natural killer cells, and myeloid subsets (Newman et al., 2015). The associations between infiltrating abundance and overall survival (OS) in TCGA LGG cohort were evaluated using univariate Cox analysis, and the survival curves were correspondingly established by Kaplan–Meier analysis.

Considering that the abundance of these infiltrating cells was mostly low, we used the mean abundance to divide the entire LGG population into high- and low-infiltrating groups for each infiltrate. For the immune infiltrates that were significantly correlated to the outcomes of LGG patients, we conducted a differential expression analysis using the R package “limma” and obtained fold change and p-value for each gene. Subsequently, for those infiltrating cells that were detrimental to patient prognosis, we selected genes ranked in the top 500 by fold change; whereas for those with a beneficial effect on patient prognosis, we selected genes whose fold change ranked in the bottom 500. Finally, we got five gene sets with 500 gene in each one.

Hub genes are derived from the intersection of the above gene sets, which were screened by using univariate Cox regression analysis. Thereafter, we used the R package “glmnet” to conduct Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis (with the penalty parameter estimated by 10-fold cross-validation), we developed an immune infiltration-related signature (IIRS) for the LGG patients. The risk score calculating formula is:

where “n” means the number of genes included in the model, “β_i” means the LASSO coefficients, “xi” is the expression level of each model gene.

Risk scores were subsequently computed for all patients included in our study. For both cohorts, the patients were divided into high- and low-risk groups according to the median risk score. Then risk plots, scatter diagrams, heatmaps, survival curves, and time-dependent receiver operating characteristic curves (ROC) were plotted using the R package ‘‘ggplot2.’’ The relationships between risk signature and survival as well as other clinicopathological characteristics were also assessed. We calculated the correlations among the signature genes, as well as the correlations between individual gene and OS and progression free survival (PFS) in LGG populations. A principal components analysis (PCA) analysis, based on Gene Expression Profiling Interactive Analysis (GEPIA2¹) webtool (Tang et al., 2019), was performed to examine the resolving power of the IRRS.

Twenty-four samples were collected and then immediately stored in liquid nitrogen. Total RNA was extracted from the tissues using TRIzol (Invitrogen, Carlsbad, CA, United States) following the instructions. The mRNA library was then constructed after quantification using NanoDrop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, MA, United States). Total RNA was purified and fragmented into small pieces for cDNA synthetization. The cDNA fragments were further amplified by polymerase chain reaction after incubating with A-tailing mix and RNA Adapter Index for end repair. The qualified double-stranded PCR products were then used to construct the final library. Eventually, the 24 qualified glioma samples were further sequenced on a BGISEQ-500 platform (BGI-Shenzhen, China). The gene expression levels were calculated using RSEM (v1.2.12). The sequencing and clinical data of these samples were summarized in Supplementary Table 1.

The efficacies of four therapies (radiotherapy, chemotherapy, targeted therapy, and immunotherapy) in high-risk and low-risk groups were evaluated. We used Tumor Immune Dysfunction and Exclusion (TIDE²) algorithm to assess the ability of the IIRS in predicating response to the immunotherapy.

Two resources for therapeutic biomarker discovery in cancer cells, including Cancer Therapeutics Response Portal (CTRP) and Genomics of Drug Sensitivity in Cancer (GDSC) (Yang et al., 2012; Rees et al., 2016), were employed to evaluate the relationship between drug sensitivity (IC50) and mRNA expression.

The Gene Set Variation Analysis (GSVA) score represents the integrated level of the expression of model gene set. The GSVA score of each patient was calculated using R package “GSVA” (Hänzelmann et al., 2013).

Data of reverse phase protein array (RPPA), a high-throughput antibody-based technique with the procedures similar to that of western blots, were used to calculate pathway activity score of 10 cancer related pathways, including TSC/mTOR, RTK, RAS/MAPK, PI3K/AKT, Hormone ER, Hormone AR, EMT, DNA Damage Response, Cell Cycle, Apoptosis pathways.

The infiltrates of 24 immune cells were evaluated through Immune Cell Abundance Identifier (ImmuCellAI) (Miao et al., 2020).

We calculated the degree to which the inputted gene set is overrepresented at the top or bottom of all genes ranked by gene expression fold change between high- and low-risk groups. The GSEA calculation was performed based on R package ‘‘fgsea.’’ The mutational profiles of model genes were assessed in cBioPortal website.³

For the differential expressed genes (DEGs) between high- and low-risk groups, protein–protein interaction enrichment analysis has been carried out with the Metascape databases integrating STRING (Szklarczyk et al., 2019), BioGrid (Stark et al., 2006), OmniPath (Türei et al., 2016), and InWeb_IM (Li et al., 2017). The resultant network contains the subset of proteins that form physical interactions with at least one other member in the list. The Molecular Complex Detection (MCODE) algorithm (Bader and Hogue, 2003) has been applied to identify densely connected network components.

Kaplan–Meier curve and log-rank test were used to compare the survival between various subgroups. The Student’s t-test was used to compare the risk scores between pairs of subgroups based on the following clinicopathologic features: age at initial pathologic diagnosis (≤40 vs. > 40 years old), gender (male vs. female), WHO grade (II vs. III), histological type (astrocytoma, oligoastrocytoma, and oligodendroglioma), and IDH1 status (mutant vs. wild-type). Wilcoxon test and Kruskal–Wallis test were used for comparison between two groups, and for comparison among more than two groups, respectively. p < 0.05 was the significance threshold in most analyses. The statistical analyses were achieved by using R language (version 4.0.3).

We obtained data from two databases (TCGA, n = 508; CGGA, n = 444) to evaluate the abundances of immune infiltrates in LGGs. The characteristics of patients in the two cohorts were summarized in Table 1. The distribution and correlations of immune infiltrates in the TCGA cohort was displayed as Figures 1A,B. Among the 21 infiltrates, the abundances of five immune infiltrates were correlated with OS in grade II and III glioma patients (Figures 1C–G). Specifically, higher abundances of naïve CD4+ T cells (p = 0.002), activated mast cells (p = 0.015), and monocytes (p = 0.014) were related to better prognosis, while higher abundances of resting memory CD4+ T cells (p = 0.015) and M1 macrophages (p = 0.020) correlated with poorer OS.

According to the mean value of abundances of the five immune infiltrates, DEGs were obtained, respectively (Figures 2A–E). We profiled 500 candidate genes in these five groups, respectively. And 44 hub genes were finally obtained when overlapping the five sets containing 2,500 genes (Figure 2F). All 44 genes were significantly associated with patient outcome, and 10 genes were profiled by LASSO regression analysis to construct the IIRS (Figures 2G,H). Detailed descriptions, LASSO coefficients, and hazard ratios for all model genes are summarized in Table 2. And we exhibited the detailed expression levels as well as prognostic curves of these genes in Supplementary Figure 1. Briefly, all ten genes except ACTN1 were differentially expressed between tumor and normal tissues. Among the remaining nine genes, FABP and PLAT were decreased in tumor tissues, while the others exhibited the opposite distribution (Supplementary Figure 1A). Regarding survival, these model genes were risk factors for OS, PFS, and disease specific survival (DSS) in patients with LGG (Supplementary Figure 1B). And the specific Kaplan–Meier survival curves (for OS) were displayed in Supplementary Figure 1C.

Risk plots, survival distributions, and model gene expressions were plotted in Figure 3A. Kaplan–Meier survival curve indicated that LGG patients with higher risk scores had significantly worse outcomes in the training set (p < 0.0001, Figure 3B). The time-dependent ROC curve demonstrated a promising ability of the model to predict OS in the training cohort (1-year AUC = 0.66, 3-year AUC = 0.69, 5-year AUC = 0.78; Figure 3C). The results were similar in the external CGGA cohort (Figure 3D). Higher risk scores also indicated poorer OS (p < 0.0001, Figure 3E). The risk model retained stable and high predication ability (1-year AUC = 0.69, 3-year AUC = 0.70, 5-year AUC = 0.74; Figure 3F). These results showed that the IIRS had a robust and stable OS-predictive ability for LGG patients. Furthermore, we performed a stratification analysis and found that the risk model maintained the ability to predict OS in most subgroups in both cohorts (Figures 3G,H).

Sankey diagrams were displayed showing the distribution of risk scores and clinicopathologic characteristics among LGG patients (Figures 4A,B). In the TCGA cohort, LGG patients with higher WHO grade had higher risk scores, while the risk score was not associated with gender. Besides, an individual patient would have a higher risk score if he had a pathologic type of astrocytoma (Figure 4C). In the CGGA cohort, risk scores were higher in patients with 1p/19q codeletion (Figure 4D). Importantly, LGG patients with wild-type IDH1 would have higher risk scores in both cohorts. Further Cox analyses showed that higher age, higher grade, pathological type (astrocytoma), and higher risk score were significantly associated with worse survival in both cohorts (Figures 4E,F).

Using samples from Xiangya hospital, we found the consistent results that glioma patients with histological type of astrocytoma, wild-type IDH1, and unmethylated MGMT had higher risk scores (Figure 5A). Although no differences in tumor purity were observed between low- and high-risk groups (Figure 5B), we found that, in our samples, high risk patients had higher abundances of M0 macrophages and activated mast cells, and lower abundances of naïve B cells, monocytes, M2 macrophages, resting mast cells, and neutrophils (Figure 5C).

Furthermore, we investigated the efficacies of multiple treatments in low- and high- risk groups (Figures 5D,E). Chemotherapies demonstrated high efficacies in both groups, while radiation and targeted therapies did not improve patients’ outcome. Intriguingly, patients in the low-risk group uniquely responded well to immunotherapy. To better evaluate the potential of our IIRS in predicting patients’ responses to immunotherapies, we employed TIDE algorithm to compare our model with other published biomarkers in immunotherapy response prediction (Figure 5F). Compared with recognized signatures (or genes) including Merck18, TIDE score, microsatellite instability (MSI) score, tumor mutation burden (TMB), CD274, CD8, IFNG, T clonality, and B clonality, the IIRS showed robust ability in predicting the response to immunotherapies in patients with different cancers. Besides, the correlations between drug sensitivity and model gene mRNA expressions were exhibited in Supplementary Figures 2A,B.

To better understand how the IRRS participates in oncologic processes, we analyzed our model both intrinsically and extrinsically. Strong and significant correlations existed among the model genes (Figure 6A), all of which were risk factors for the OS in the LGG cohort (Figure 6B). Besides, the expressions of all genes, except for LGALS3, were correlated with poorer PFS in patients with LGG (Figure 6B). The comprehensive investigation in the mutational profiles of model genes was displayed as Supplementary Figure 3, where we found BGN was most likely to mutate in LGGs. These mutations were significantly related to molecular status including IDH, 1p, and 19q. A subsequent PCA analysis suggested the IRRS could distinguish LGG from normal cerebral cortex tissues and glioblastomas (Figure 6C).

Next, we evaluated the involvement of individual genes in various cancer related pathways, and found that all genes participated in activated “epithelial-mesenchymal transition (EMT)” pathway while five of them were involved in inhibited “DNA damage” pathway (Figure 6D). Furthermore, GVSA score was significantly and positively correlated with EMT pathway (r = 0.27, p = 9.6e-9) while negatively correlated with “Hormone androgen/androgen receptor (AR)” pathway (r = −0.16, p = 5.1e-4) (Figure 6E).

We then focused on the association of the IRRS and multiple immune cell types (Figure 6F). Strong and positive associations were observed between GVSA scores and activated T lymphocytes including cytotoxic T lymphocytes (r = 0.61, p = 3.6e-60) and type 1 helper T cells (r = 0.51, p = 1.3e-39). In contrary, the GVSA score was strongly and negatively correlated with the abundance of neutrophil (r = −0.47, p = 2.4e-33).

We obtained DEGs between low- and high-risk groups (Figure 7A), which were displayed as volcano plot (Figure 7B). Functional analysis indicated these DEGs were mainly located in extracellular matrix and were involved in binding functions such calcium ion, kinase, and cell adhesion molecule binding (Figure 7C). As for biological processes, they participated in responses to wounding and inflammatory responses (Figure 7D). Furthermore, the gene sets were mainly enriched in several Reactome pathways including “cytokine signaling in immune system,” “hemostasis,” and “neutrophil degranulation” (Figure 7E). To further capture the relationships between the terms, a subset of enriched terms has been selected and rendered as a network plot, where terms with a similarity above 0.3 are connected by edges. We select the terms with the best p-values from each of the 20 clusters and visualized the network as Supplementary Figures 2C,D (Shannon et al., 2003).

The MCODE networks identified for individual gene lists have been gathered and are shown in Figures 7F,G. Pathway and process enrichment analysis has been applied to each MCODE component independently, and the three best-scoring terms by p-value have been retained as the functional description of the corresponding components: (1) GO:0002399, MHC class II protein complex assembly (log₁₀p = −15.7); (2) GO:0002503, peptide antigen assembly with MHC class II protein complex (log₁₀p = −15.7); and (3) GO:0009611, response to wounding (log₁₀p = −15.4).