From TCGA, The RNA-Seq data and clinical information of LGG and GBM samples were collected. Molecular subtype and treatment details were retrieved from Ceccarelli et al.’s work (22), while survival status and demographic information were accessed using the RTCGAclinical R package (version: 1.16.0). RNA-Seq data was measured as Fragments Per Kilobase of transcript per Million mapped reads (FPKM) and log2-based transformation. Then, the sva package (23) was utilized for the normalization of RNA expression profiles and to remove the batch effects between TCGA-LGG samples and TCGA-GBM samples. For external validation, normalized RNA-Seq data via RNA-Seq by Expectation Maximization (RSEM) (24) and clinical information were obtained from the CGGA data set. All data sets were downloaded on June 1st, 2020.

Stromal Score (SS) and Immune Score (IS) were calculated by the Estimation of Stromal and Immune cells in Malignant Tumors using Expression data (ESTIMATE) algorithm (18) in R software loaded with the estimate package. The scores represent the stromal and immune components in TME respectively. Glioma samples in TCGA were divided into high-score and low-score groups with the median of SS and IS respectively. Samples with positive values of survival time were used for survival analysis between groups with survival and survminer packages (25).

The Linear Models for Microarray Data (LIMMA) package (26) in R software was used to replace replicates of probes with their average, while the Wilcoxon test was applied to extract DEGs between groups as defined by medians of SS and IS respectively. P < 0.05 and |log2FC| > 1 were set as the threshold for DEGs. The heatmap of the 100 DEGs with the most significant P values was plotted. Feature DEGs were identified as unanimously upregulated or downregulated DEGs in both the high-SS and high-IS groups. DEGs in independent cohorts (TCGA-LGG/TCGA-GBM) were also identified using the same methods.

The clusterProfiler package (27) was utilized to perform the Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for feature DEGs. Three categories were included in the GO enrichment analysis, i.e., biological process (BP), cellular component (CC), and molecular function (MF); whereas, KEGG revealed enriched pathways related to these feature DEGs. Subsequently, the PPI network was constructed via Search Tool for the Retrieval of Interacting Genes (STRING) (28) to identify hub genes. Visualization of the network was performed in Cytoscape 3.8.0 (29, 30).

In the context of feature DEGs, univariate cox regression was conducted to identify prognostic DEGs, which was subsequently utilized in gene selection by the least absolute shrinkage and selection operator (LASSO) regression analysis (31) in the TCGA-LGG. And then we calculated the individualized risk score with coefficient-weighted gene expressions and constructed a PROMISE model with the following formula. Samples were divided into high-score and low-score groups based on the median risk score. Then the clinical relevance was validated using survival analysis between groups with the threshold of p < 0.05. The PROMISE model was then validated in the TCGA-GBM samples and the TCGA combined set with the same risk score formula and cutoff value. The receiver operating characteristic (ROC) analysis was performed, and we calculated the area under the curve (AUC) to evaluate the predictive accuracy (32) in the TCGA samples. Subsequently, a risk plot was presented that includes a heatmap presenting expression profiles of included genes, a distribution plot of risk scores based on the model, and a dot plot showing the survival status of patients in different risk groups.

The PROMISE model was then validated using survival analysis in CGGA-LGG, CGGA-GBM, and CGGA data sets respectively. The Log-rank p-value was calculated while the Kaplan-Meier graph was plotted. P-values less than 0.05 were considered statistically significant. The prognostic performance was validated in the CGGA with ROC analysis as well as a risk plot.

Clinicopathological factors were collected from TCGA data set and integrated with transcriptome profile derived from TCGA data set. The univariate and multivariate Cox regressions were conducted to determine whether the PROMISE model was independent of clinicopathological factors as well as to identify other independent prognostic factors. Based on independent prognostic factors, we formulated a nomogram to predict the survival probability of individual patients using methods reported in previous literature (33, 34). The ROC analysis for nomogram-based prediction was performed in both TCGA and CGGA data sets.

We performed Gene Set Enrichment Analysis (GSEA) between high-risk and low-risk groups as separated by the PROMISE model via clusterProfiler and enrichplot packages. The gseKEGG function was applied to identify the enriched pathways in KEGG.

To characterize the immune TME in different risk groups as defined by the PROMISE model, the Cibersort (35) method was adopted in our study. Immune cells with significant differences between risk groups were identified in TCGA data set and CGGA data set, respectively. We highlighted cells that were unanimously upregulated or downregulated in the high-risk groups in both data sets. Subsequently, Spearman correlation analysis was performed using the abundance of these cells and the PROMISE risk score.

The workflow of our study is shown in Figure 1 . We obtained the RNA-Seq data and clinical information of 529 LGG samples and 169 GBM samples from TCGA, out of which 604 cases were recorded with positive values of survival time. Likewise, normalized RNA-Seq data and clinical information for 443 LGG samples and 249 GBM samples were collected from the CGGA. Among them, 657 cases had survival time>0. The distribution of cases concerning the clinicopathological factors was presented in Table 1 .

With the ESTIMATE algorithm, we calculated SS and IS for all these TCGA samples. To explore the potential association of survival rate in glioma patients with the stromal and immune components in TME, we classified these samples into high-score and low-score groups with the median scores. Survival analysis revealed that the survival time of samples with high SS was shorter than that of samples with low scores, yet no significant difference was detected ( Figure 2A , p=0.109). The survival analysis presents the same trend between groups as separated by the median of IS with a statistical difference ( Figure 2B , p<0.001).

According to the threshold, we identified DEGs based on SS or IS scores respectively. The heatmap of the 100 DEGs with the most significant P values was plotted in Figures 2C, D . In order to identify key DEGs on different TME, we obtained 220 unanimously upregulated DEGs and 42 unanimously downregulated DEGs in both high-IS and high-SS groups ( Figure 2E ). The subsequent analysis was based on these feature DEGs, which can be accessed in Supplementary File S1 . In parallel, we also identified DEGs in TCGA-LGG and TCGA-GBM, respectively ( Supplementary Files S2 and S3 ).

The GO analysis ( Figure 3A ) demonstrated that these feature DEGs were significantly enriched the BP that includes regulation of immune effector process, leukocyte migration, cell chemotaxis, and so on. The CC analysis indicated that these feature DEGs were present on the external side of the plasma membrane, collagen-containing extracellular matrix, secretory granule membrane, etc. The feature DEGs were enriched in MF, such as G protein-coupled receptor binding, cytokine receptor binding, and cytokine activity. Furthermore, the KEGG analysis indicated that the genes were mainly involved in Cytokine-cytokine receptor interaction, which was closely related to the immune response ( Figure 3B ). Then we constructed a PPI network to further explore the interplay among the feature DEGs and obtained 40 hub genes with the most interactions ( Figures 3C, D ).

Out of the 262 feature DEGs, a total of 155 prognostic DEGs were identified using univariate cox regression ( Supplementary File S4 ). Based on these prognostic DEGs, we selected a total of 4 specific DEGs (CD86, ANXA1, C5AR1, and CD5) in the TCGA-LGG with Lasso regression ( Figures 4A, B ). After extracting the coefficient values, we calculated the PROMISE risk scores with coefficient-weighted expression levels of 4 hub genes with the following formula:

LGG samples were divided into high-risk and low-risk groups based on the median risk score (0.20). Significant difference was shown between groups with survival analysis ( Figure 4D , HR = 4.3, 95% CI = 2.8 - 6.5, p < 0.001). Subsequently, we calculated individualized risk scores in TCGA-GBM samples and divided them into high-risk and low-risk groups with the same cutoff value. A significant difference was also presented between groups with survival analysis in TCGA-GBM ( Figure 4C , HR = 1.8, 95% CI = 0.97–3.4, p=0.028). Then we repeated the same validation method in the TCGA combined set and yielded distinct separation in survival outcomes with p<0.001(HR = 2.6, 95% CI = 1.8–3.5) ( Figure 4E ). ROC curve analysis of the PROMISE model achieved AUC of 0.685, 0.619, and 0.621 at 1-, 3-, and 5-year ( Figure 4F ). Figure 4G shows the risk plot in TCGA glioma cohort, including the distribution of risk score and accurate classification of survival outcomes across different risk groups with red dots being dead and blue ones being living cases. The expression levels of four hub genes demonstrated a trend of upregulation of ANXA1 along with downregulation of CD86, C5AR1, and CD5 in the high-risk group as classified by the PROMISE model ( Figure 4G ).

In order to determine if the PROMISE model was independent of data sets, we validated the model in the CGGA cohort. Survival analysis showed significant differences in CGGA-LGG, CGGA-GBM, and CGGA combined set with p<0.001, p=0.015, and p<0.001 respectively ( Figures 5A, B, C ). ROC curve analysis of the PROMISE model in CGGA Glioma samples achieved AUCs of 0.677, 0.737, and 0.769 at 1-, 3-, and 5-year ( Figure 5D ). Likewise, the risk plot of CGGA encompassing the distribution of the PROMISE risk score, survival status of cases as well as gene expression of 4 hub genes was presented with a similar expression pattern to TCGA data set ( Figure 5E ).

Multivariate analysis revealed that the PROMISE model based on 4 hub genes could be an independent prognostic marker after clinical characteristics were adjusted. Likewise, age and tumor grade were also independent prognostic factors ( Figure 6A ). Based on the independent prognostic factors, we formulated a nomogram to predict the survival probability of individual samples ( Figure 6B ). Individualized survival chance could be predicted intuitively with specific age, tumor grade, Isocitrate dehydrogenase (IDH) mutation status, and the PROMISE risk score. ROC curve analysis showed adequate discrimination with an AUC of 0.917 and 0.817 in TCGA Glioma samples ( Figure 6C ) and CGGA Glioma samples ( Figure 6D ) respectively, which outperformed all clinicopathological factors including classifications of molecular subtype.

To explore the underlying mechanism, we performed GSEA between groups based on the PROMISE model to identify the enriched KEGG pathway. Figure 7A shows the results of the KEGG analysis in TCGA Glioma samples. Antigen processing and presentation, cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, phagosome, T cell receptor signaling pathway, and Th17 cell differentiation were significantly enriched. In CGGA Glioma samples ( Figure 7B ), antigen processing and presentation, cell adhesion molecules, cytokine-cytokine receptor interaction, ECM- receptor interaction, cell differentiation of Th1, Th2, and Th17 were significantly enriched.

Distribution levels of 22 immune cells based on risk groups in TCGA Glioma samples and TCGA Glioma samples were shown in Figures 8A, C . Immune cells with significant differences between risk groups were shown in Figures 8B, D respectively. Cells that were unanimously upregulated or downregulated in both data sets were identified. The specific immune cells were marked in red. It is shown that T cells CD8 (p=0.002, p<0.001) and Neutrophils (p<0.001, p<0.001) level were significantly enriched in the high-score group while lower T cell follicular helper (p=0.004, p=0.008) and natural killer (NK) cells resting (p<0.001, p<0.001) were observed in high-risk groups in both TCGA and CGGA data sets. Spearman correlation analysis indicated that these cells were correlated with the PROMISE model in both data sets ( Supplementary File S5 ).