Prognostic and predictive value of an immune infiltration signature in diffuse lower-grade gliomas

Gliomas are the most common primary malignant brain tumors, and can be classified into 4 grades (grades I–IV) according to the WHO grading system (1, 2). Diffuse low-grade and intermediate-grade gliomas (WHO grades II and III), including astrocytomas, oligodendrogliomas, and mixed oligoastrocytomas, are considered lower-grade gliomas (LGGs) (3). Surgical resection combined with chemoradiotherapy is the most frequent treatment (4). Owing to their highly invasive nature, complete surgical resection is difficult, and the residual tumor can result in recurrence and even malignant progression (4, 5). Survival outcome ranges widely across patients. Some LGGs may progress to glioblastoma quickly, while some remain stable for a long time.

Although the histopathological classification of LGGs is widely recognized, it cannot adequately predict survival. Therefore, clinicians tend to rely on genetic classifications to guide treatment. Recurrently mutated genes like IDH1, IDH2, TP53, EGFR, and ATRX are well-recognized factors for the prognosis of patients with LGGs in clinical practice (6–8). Other molecular markers, including 1p/19q codeletion and MGMT promoter methylation, are also important prognostic factors for LGGs (9). However, these clinicopathological and genetic factors fail to evaluate survival outcomes accurately. Patients with the same risk factors might have conflicting outcomes. Consequently, a more comprehensive study is needed to increase the prognostic and predictive accuracy of the current assessment system.

Numerous studies have provided evidence that cancer progression and recurrence are not only driven by the tumor’s underlying genetic changes, but also by the tumor microenvironment (TME) (10–12). Increasing evidence has confirmed that immune cells in the TME are involved in tumor progression and recurrence (13, 14). Furthermore, the effects of infiltrating immune cells on clinical outcome have been widely reported (15–17). Thorsson et al. conducted immunogenomics analyses of more than 10,000 tumors and identified 6 immune subtypes to define immune response patterns impacting prognosis (18), prompting that the immunogenomics may serve as a resource for understanding tumor-immune interactions and improving clinical management. Recently, the prediction of long-term outcome based on gene expression profile analysis has shown remarkable prospect. Charoentong et al. analyzed 8243 TCGA samples from 20 solid tumors and defined a set of 782 pan-cancer metagenes for 28 immune cell subpopulations to evaluate the intratumoral immune infiltration landscape (19). Based on Charoentong’s findings, the specific involvement of immune infiltration–related genes in LGGs was identified with high-throughput technology and used to build a model to more accurately predict survival outcome. In addition, the predictive model could be utilized to guide the postoperative chemoradiotherapy. Patients predicted with high risk may be more suitable for close imaging monitoring and radical postoperative adjuvant therapy.

Integrating multiple gene markers into a single model would significantly improve the accuracy and robustness of prediction compared with using a single marker. Chai et al. analyzed the RNA processing genes in LGGs and identified a 19-gene risk signature, which had better prognostic value than the traditional factors (20). In this study, we developed a risk signature based on intratumoral immune infiltration–related genes with the least absolute shrinkage and selection operator (LASSO) Cox regression model and established a nomogram that incorporated the immune-related risk score and clinical factors to predict the survival of LGG patients. In addition, we assessed the prognostic and predictive accuracy of this model in 1 internal validation set and 4 external validation sets.

Patient characteristics. The study flowchart is presented in Figure 1. All the included patients were pathologically diagnosed as LGG. The characteristics of patients in the training set and validation sets were presented in Figure 2 and Supplemental Tables 1 and 2 (supplemental material available online with this article; https://doi.org/10.1172/jci.insight.133811DS1). A total of 1216 LGG patients with a mean age of 42.0 ± 12.4 years were entered according to the screening criteria. There were 434 events (deaths) over a median follow-up time of 3.2 years (range, 1.2 months to 20.7 years).

Figure 1

Study flowchart. LGGs, lower-grade gliomas; DEIRs, differentially expressed immune-related RNAs; TPM, transcripts per million.

Figure 2

Association of clinical factors with overall survival based on univariate Cox analysis in the training set and validation sets. A, astrocytoma; NA, not available; NA*, not applicable; OA, oligoastrocytoma; OD, oligodendroglioma.

Immune infiltration landscape in LGGs. To analyze the immune infiltration landscape in LGGs, single-sample gene set enrichment analysis (ssGSEA) was applied to evaluate the relative abundance of each immune cell subpopulation using RNA sequencing (RNA-Seq) data of 469 LGGs from The Cancer Genome Atlas (TCGA). Infiltrating immune cells were found to be significantly heterogeneous in different LGGs. Of the 28 cell subpopulations, 22 were correlated with prognosis (P < 0.1) in univariate Cox survival analysis (Supplemental Table 3). For further characterization, unsupervised clustering was implemented to categorize the 469 patients into 2 infiltration subgroups termed as type A (n = 253) and type B (n = 216) infiltration based on the 22 immune cell subpopulations (Figure 3). Kaplan-Meier analysis for overall survival (OS) showed that patients with type B infiltration had a worse prognosis (P = 0.0028, Supplemental Figure 1); this was in accordance with the fact that infiltrating immune cells in the TME would affect the clinical outcome.

Figure 3

The immune infiltration landscape of LGGs. Unsupervised clustering of 469 LGGs from TCGA using ssGSEA scores based on the 22 survival-related immune cell subpopulations. Age, sex, grade, IDH mutation status, 1p/19q codeletion status, and histological type are annotated in the lower panel. Two distinct immune infiltration subgroups termed as type A infiltration and type B infiltration were defined.

The characteristics of the 469 patients by immune infiltration type are summarized in Supplemental Table 4. An interesting finding was that LGGs with type A and type B infiltration differed significantly in 1p/19q status and histological type. More than 80% of LGGs with 1p/19q codeletion were of type A infiltration. In the correlation analyses between clinical variables, genetic variables, and immunological subtype (Supplemental Table 5), type A infiltration presented notable positive correlations with 1p/19q codeletion (Spearman’s rho = 0.397, P < 0.001) and histological oligodendroglioma (Spearman’s rho = 0.330, P < 0.001).

Selection of an immune infiltration signature and building of a risk score. Of the 782 immune metagenes, 273 were identified as differentially expressed through the analysis of 469 LGG samples from TCGA and 469 normal brain specimens from the Genotype-Tissue Expression (GTEx) Portal. Of the 273 differentially expressed immune-related RNAs (DEIRs), 179 and 143 were respectively filtered from 469 LGGs of TCGA and 405 LGGs of Chinese Glioma Genome Atlas (CGGA) by the criteria of maximal expression > 2 transcripts per million (TPM) and univariate Cox analysis P < 0.05. A total of 120 DEIRs were common to both data sets (Supplemental Table 6). After learning in the training set, the 120 DEIRs were reduced to the 20 most powerful prognostic markers with nonzero coefficients in the LASSO Cox regression model (Supplemental Figure 2 and Table 1). The formula for the risk score is presented in Supplemental Figure 3.

Table 1

Information of the risk score model based on the 20 survival-related immune metagenes

Enrichment analysis of the 120 survival-related DEIRs identified overrepresented biological processes in gene ontology (Supplemental Figure 4). Most of the biological processes were related to the activation and proliferation of immune cells, indicating that intratumoral immune infiltration played an important role in the prognosis of LGGs.

Prognostic value of the risk score. The predictive accuracy of the 20-DEIR–based risk score was assessed with time-dependent receiver operating characteristic (ROC) analysis on OS at 3 and 5 years in all data sets, and the results revealed that the risk score was a valuable predictor (Figure 4). To compare the predictive accuracy of the risk score with that of traditional clinical factors, time-dependent ROC analysis was applied to the 1216 LGGs of the whole data set (Supplemental Figure 5). The risk score gave the highest 3-year and 5-year AUC (0.814 and 0.765, respectively), indicating that it was a better predictor of survival for LGG patients.

Figure 4

Kaplan-Meier survival analyses and time-dependent ROC curves based on the risk score in the training and validation sets. (A–F) Training set (A), internal validation set (B), external validation set 1 (C), external validation set 2 (D), external validation set 3 (E), and external validation set 4 (F). Patients with a risk score of more than 0.6907 were assigned to the high-risk group and the rest were assigned to the low-risk group. P values were calculated based on the log-rank test. AUCs at 3 and 5 years were used to assess predictive accuracy.

To further assess the prognostic value of the risk score, high-risk and low-risk patients were divided based on the optimum cutoff score generated by the ROC curve for predicting 5-year survival in the training set, and Kaplan-Meier analysis was performed in the 2 groups. Patients with a risk score of more than 0.6907 were assigned to the high-risk group, and the rest were assigned to the low-risk group. As expected, patients in the high-risk group had a worse OS (P < 0.001; Figure 4A). To confirm whether the risk score had stable prognostic value in different populations, we applied it to an internal validation set and 4 external validation sets using the same score model and cut-off value. The results are presented in Figure 4.

In addition, the 2016 WHO classification of central nervous system tumors categorized LGGs into 3 distinct subtypes based on their molecular features: IDH WT, IDH mutation and 1p/19q codeletion, and IDH mutation and 1p/19q non-codeletion (1). Therefore, we performed a subsequent Kaplan-Meier analysis of the patients with high- and low-risk LGGs based on the molecular subtypes. The results revealed that the 20-DEIR–based risk score remained effective in stratifying survival within different subtypes (Figure 5). Moreover, Kaplan-Meier analyses of the patients with high- and low-risk LGGs based on clinical factors including age, sex, grade, and histological diagnosis were also performed. The results further confirmed the robust stratification ability of the 20-DEIR–based risk score (Figure 6). Kaplan-Meier survival analyses stratified by immune-related risk score, age, grade, IDH mutation, 1p/19q status, and histological type were also performed in the whole data set. The hazard ratio (HR) value of immune-related risk was the highest, indicating that its stratification ability was superior to other prognostic factors (Supplemental Figure 6).

Figure 5

Kaplan-Meier curves based on the IDH and 1p/19q status in the whole data set. (A–F) Kaplan-Meier curves of OS of patients with high- and low-risk LGGs based on the IDH and 1p/19q status for 764 patients with IDH mutation (A), 257 patients with WT IDH (B), 349 patients with 1p/19q codeletion (C), 690 patients without 1p/19q codeletion (D), 276 patients with IDH mutation and 1p/19q codeletion (E), and 382 patients with IDH mutation and without 1p/19q codeletion (F). Log-rank test was used for statistical findings.

Figure 6

Kaplan-Meier curves based on clinical factors in the whole data set. Kaplan-Meier curves of OS of patients with high- and low-risk LGGs based on age, sex, grade, and diagnosis for 604 patients with age > 40 years (A), 610 patients with age ≤ 40 years (B), 621 male patients (C), 459 female patients (D), 608 patients with grade III gliomas (E), 607 patients with grade II gliomas (F), 455 patients diagnosed with astrocytomas (G), 320 patients diagnosed with oligodendrogliomas (H), and 441 patients diagnosed with oligoastrocytomas (I). Log-rank test was used for statistical findings.

Also, the individualized risk score may be utilized to guide the postoperative chemoradiotherapy. Based on the immune-related risk score, patients were divided into high- and low-risk groups. Multivariate Cox analysis was conducted in high- and low-risk groups to identify the prognostic factors. Radiotherapy was found to be a favorable prognostic factor in the high-risk group (HR, 0.44; 95% CI, 0.21–0.91; P = 0.027), but it made no sense in the low-risk group (HR, 1.05; 95% CI, 0.57–1.92; P = 0.882; Supplemental Table 7).

Furthermore, the relationship between the risk score and biomarkers for checkpoint inhibitor immunotherapy was evaluated. The risk score was positively correlated with the expression levels of PD-L1 and TGF-β, especially in grade III gliomas (Supplemental Figures 7–9).

Biological processes associated with the risk score in LGGs. Based on the immune infiltration–related risk score, LGGs were significantly stratified into high- and low-risk groups. For further exploring the biological processes associated with the risk score, we performed GSEA to identify the enriched cancer hallmarks in TCGA and CGGA RNA-Seq data sets (Figure 7). We found that the essential pathways correlated with cell proliferation and invasion were enriched in the both data sets, such as G2/M checkpoint, E2F targets, epithelial-mesenchymal transition, and angiogenesis. Moreover, the pathways correlated with immune and inflammatory response were significantly enriched, including IFN-γ response, IFN-α response, allograft rejection, inflammatory response, TNF-α signaling via NF-κB, IL-6/JAK/STAT3 signaling, and IL-2/STAT5 signaling. These findings indicate that the intratumoral immune infiltration plays an important role in the malignant progression of LGGs.

Figure 7

Biological processes associated with the risk score in LGGs. (A and B) Use of the GSEA to identify the biological processes positively correlated with the risk score in TCGA (A) and CGGA RNA-Seq (B) data sets. The pathways were arranged according to the normalized enrichment scores (NES). padj, FDR-adjusted P value.

Development of a predictive nomogram for OS. After backward elimination of variables using Akaike’s information criterion (AIC) as a stopping rule, risk score, age, and grade were included finally, and a nomogram incorporating the 3 factors was generated to predict the 3-year and 5-year OS in the training set (Figure 8A). The calibration curves suggested that the nomogram was well-calibrated across all data sets (Figure 8, B–G). The C-indices (0.873 for training set, 0.881 for internal validation set, 0.781 for external validation set 1, 0.765 for external validation set 2, 0.721 for external validation set 3, and 0.753 for external validation set 4) indicated the good discriminative ability of our model.

Figure 8

Predictive nomogram and calibration curves. (A) Nomogram to predict the 3-year and 5-year OS. (B–G) Calibration plots for the nomogram model in the training set (B), internal validation set (C), external validation set 1 (D), external validation set 2 (E), external validation set 3 (F), and external validation set 4 (G). The 45-degree dotted line represents a perfect prediction, and the solid lines reflect the predictive performance of the nomogram.

In this study, we found that the gene expression pattern of intratumoral immune infiltration was associated with the malignancy of LGGs and identified a 20-DEIR–based risk signature significantly correlated with the survival of LGG patients. We further built a nomogram that incorporated the 20-DEIR–based risk score, patient’s age, and tumor grade to predict the OS of patients with LGGs after surgical resection. Lower-grade brainstem gliomas were excluded from our study due to their unique clinical characteristics, low morbidity, and few cases. The performance of the nomogram was verified in all data sets, which guaranteed the repeatability of our model. Compared with a previously proposed long noncoding RNA–based (lncRNA-based) signature (21), the number of cases included in this study was larger and the predictive AUC was higher, indicating that our immune-related prognostic signature was more accurate. Also, several predictive nomograms were established in previous studies to predict the long-term OS of LGGs (22, 23). In comparison, the more adequate validation with consistently high C-indices of our nomogram indicated the improvement of reliability.

In the analysis of immune infiltration landscape, we found that LGGs had 2 immune infiltration subtypes. LGGs with type A and type B infiltration differed significantly in 1p/19q status and histological type. In the correlation analysis, type A infiltration presented notable positive correlations with 1p/19q codeletion (Spearman’s rho = 0.397, P < 0.001) and histological oligodendroglioma (Spearman’s rho = 0.330, P < 0.001). As we know, the 1p/19q codeletion has been proposed as a diagnostic standard for oligodendroglioma in the 2016 WHO classification of tumors of the central nervous system (24), which is in accordance with our analysis. In other words, the infiltrating immune cells in oligodendroglioma are different from those in astrocytoma. Oligodendroglioma tends to type A infiltration, while astrocytoma tends to type B infiltration, which may contribute to the difference in prognosis. However, the driving factors of the difference in the infiltrating immune cells are still unknown, which requires further study.

Using LASSO Cox regression, 120 candidate survival-related DEIRs were reduced to the 20 most powerful prognostic predictors, and a risk score model based on the 20 DEIRs was then generated. The prognostic value of the risk score was verified in time-dependent ROC analyses of all data sets. Furthermore, the 20-DEIR–based risk score was identified to have a higher accuracy in predicting 3-year and 5-year OS compared with traditional stratifying factors (age, WHO grade, IDH mutation, 1p/19q status, and histological type). In the multivariate Cox analysis, after adjusting for clinical factors, the risk score remained an independent prognostic factor (HR, 1.65 per 1 score increase; 95% CI, 1.46–1.86; P < 0.001) in the whole data set. In addition, the Kaplan-Meier analysis showed that the risk score still remained effective in stratifying patients within different molecular subtypes of LGGs based on the 2016 WHO classification, further clarifying its clinical value. Also, Kaplan-Meier survival analyses stratified by immune-related risk score, age, grade, IDH mutation, 1p/19q status, and histological type showed that the immune-related risk gave the highest HR value, indicating that its stratification ability was superior to other prognostic factors. Our findings highlight the role of intratumoral immune infiltration in the outcomes of LGGs, and the 20-DEIR–based risk signature is promising to be an effective supplement for the diagnostic criteria of LGGs in further stratifying the prognosis of patients. Moreover, the prognostic signature could be utilized to guide the postoperative chemoradiotherapy. Based on the individualized risk score, patients were divided into high- and low-risk groups. High-risk patients may be more suitable for close imaging monitoring, radical postoperative radiotherapy, and chemotherapy, while low-risk patients could consider radiotherapy alone or just observation. In order to verify the clinical value of our model, multivariate Cox analysis was conducted in high- and low-risk groups to identify the prognostic factors. Radiotherapy was a favorable prognostic factor in the high-risk group, but it made no sense in the low-risk group. Chemotherapy was not an independent prognostic factor in either group. This indicated that high-risk patients were more likely to benefit from the radiotherapy. However, further verification by prospective studies is necessary.

In addition, our risk score model might also help to assess the immune microenvironment and guide immunotherapy. Checkpoint inhibitor immunotherapy has been receiving increasing attention due to its impressive success in the treatment of various solid tumors, and it is likely to be a promising choice of treatment for LGGs in the future. TGF-β secreted by glioma cells and innate immune cells can block the antitumor response of T cells and maintain the immunosuppressive microenvironment (25, 26). As the risk score is positively correlated with the expression level of TGF-β, a high risk score is likely to reflect a more immunosuppressive glioma microenvironment. The ligand PD-L1 on the surface of tumor cells can prevent T cells from identifying cancer antigens through the binding of PD-1 (27, 28). In the TME, this interaction is an important mechanism for the immune escape of tumor cells (28). The overexpression of PD-L1 is a widely used predictive biomarker to determine response to checkpoint inhibitor immunotherapy (29). As the risk score is positively correlated with the expression level of PD-L1, LGG patients with a high risk score would be more likely to benefit from anti–PD-1/PD-L1 immunotherapy. However, resistance to immune checkpoint inhibitors occurs frequently in clinical treatment. One possible reason is the increasing immunosuppressive status in the TME owing to the activation of TGF-β pathway (30). Inhibiting TGF-β has been found to enhance the efficacy of anti–PD-1/PD-L1 treatment in numerous studies (30–32). M7824, a bifunctional fusion protein targeting PD-L1 and TGF-β, has shown encouraging signs of efficacy in a Phase I clinical trial (32). Our risk score model, which is positively correlated with the expression levels of both TGF-β and PD-L1, is expected to be a potential biomarker for the cotargeting immunotherapy. Due to the stronger correlation between the risk score and the expression levels of TGF-β and PD-L1 in grade III gliomas, grade III patients with high risk scores are more likely to benefit from the cotargeting immunotherapy. It should be pointed out that IHC is used to assess the expression of PD-L1 at the protein level in clinical immunotherapy (33). The value of PD-L1 detection at the gene level needs to be further verified.

Based on the stratification of immune infiltration–related risk, GSEA was performed. We identified the essential biological processes that were associated with the 20-DEIR–based risk score, including not only cell proliferation and invasion, but also immune and inflammatory response. Not surprisingly, 7 pathways correlated with immune and inflammatory response were enriched, further supporting the crucial role of immunity in the prognosis of LGGs.

We noted that adjuvant chemotherapy and radiotherapy were not independent prognostic factors for OS in the multivariate Cox analysis of the whole data set and were excluded from the predictive nomogram. Moreover, chemotherapy and radiotherapy were both significant risk factors in the univariate Cox analysis of the whole data set. Even so, the importance of chemoradiotherapy in the prognosis of LGGs cannot be denied. In the past, the choice and timing of adjuvant therapies for LGGs were controversial (34). Chemotherapy and radiotherapy were more likely to be recommended to patients with factors of poor outcome after surgery. In the correlation analysis of clinical factors, tumor grade was significantly correlated with chemotherapy (Spearman’s rho = 0.36, P < 0.001) and radiotherapy (Spearman’s rho = 0.23, P < 0.001). Patients with grade III gliomas were more likely to receive chemotherapy and radiotherapy. It was acknowledged that patients with grade III gliomas had worse prognosis compared with those with grade II gliomas, which might partly explain why chemotherapy and radiotherapy were not independent prognostic factors for OS in the multivariate Cox analysis.

The limitations of this study should be noted. Firstly, the retrospective nature of data collection, the lack of some data, and the failure to incorporate some recognized prognostic factors (e.g., MGMT promoter methylation, tumor size, and extent of surgical resection; refs. 9, 35, 36) indicate the need for prospective data collection and high incorporation of comprehensive factors to increase the credibility of the model. Secondly, the microarray data from the 3 data sets were used to validate the model trained from sequencing data. Due to the differences in technological principles between microarray and sequencing, systematic errors existed, although all data were normalized to [0,1]. However, from another point of view, the model also showed good accuracy in the microarray data sets, indicating the robustness of our model. Thirdly, the CGGA data set in our study was not an absolutely external validation set because it was used to obtain the survival-related DEIRs in the training phase. From another perspective, since the samples in TCGA were mainly from White people and those in CGGA were mostly from people of Asian descent, the model constructed by our method could reduce the influence of racial differences on the predictive accuracy; this supported the extensive use of this predictive model, regardless of race and region. Moreover, the prognostic value of our model was confirmed by the inclusion of 3 additional external validation sets based on microarray data. Finally, the intratumoral immune infiltration was only measured using ssGSEA at the molecular levels. The relationship between infiltration detected at the molecular levels and infiltration of actual immune cells was not evaluated in this study. If possible, we will verify it on the histological aspects in the future.

In conclusion, our study highlights the prognostic value of intratumoral immune infiltration in LGGs and developed a 20-DEIR–based risk signature with favorable predictive value and stratifying ability. Additionally, our nomogram incorporating the 20-DEIR–based risk score, patient’s age, and tumor grade performed well, with good calibration and discrimination for predicting the survival of LGG patients after surgical resection. Furthermore, the 20-DEIR–based risk score is a potential biomarker for the PD-L1 and TGF-β cotargeting immunotherapy.

Patients and data sets. RNA-Seq data and clinical information for 510 LGGs were retrieved from TCGA data portal. The most comprehensive LGG data set from TCGA was used as the training set in our study. Clinical and molecular information from the CGGA data base contained RNA-Seq data of 630 LGGs and microarray data of 177 LGGs with approval by the ethics committee of Beijing Tiantan Hospital. The 2 data sets were used as external validation sets 1 and 2. Microarray data for 109 LGGs based on GSE16011 and 137 LGGs based on GSE61374 were downloaded from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/), the clinical information of which were found in previous studies (37, 38). The 2 data sets were used as external validation sets 3 and 4.

The exclusion criteria for the downloaded data were as follows: (a) brainstem glioma, (b) recurrent glioma, (c) patients who underwent biopsy alone without tumor resection, and (d) patients with survival data ≤ 30 days. After filtering, 469 LGGs from the TCGA RNA-Seq data set, 405 LGGs from the CGGA RNA-Seq data set, 118 LGGs from the CGGA microarray data set, 88 LGGs from GSE16011, and 136 LGGs from GSE61374 were included in this study.

Evaluation of immune infiltration. A previously described method was used to evaluate the immune infiltration in LGGs (39). The 782 metagenes for 28 immune cell subpopulations involved in innate immunity and adaptive immunity were obtained from Charoentong et al. (19). The immune infiltration levels were quantified using enrichment scores calculated by ssGSEA. The ssGSEA applied gene signatures expressed by the 28 immune cell subpopulations to individual tumor samples (39). Unsupervised clustering of 469 LGGs from TCGA was performed using the calculated ssGSEA scores.

Data processing and risk score building. Differentially expressed RNAs were analyzed in 469 LGG samples from TCGA (counts value) and 469 normal brain specimens from the GTEx Portal (counts value) using the R package DESeq2. Simultaneously, to obtain the DEIRs, the 782 immune-related metagenes (19) were used as criteria for further screening. The expression differences were characterized by log₂ fold change (log₂FC) and the associated FDR-adjusted P value (adj.P). Immune-related RNAs of |log₂FC| > 1 and adj.P < 0.05 were identified as being differentially expressed and selected as DEIRs.

RNA-Seq data (fragments per kilobase per million mapped reads [FPKM] value) of the 469 LGGs from TCGA and 405 LGGs from the CGGA were converted to TPM values using a formula described in previous studies (40–42). As a considerable expression abundance is necessary for gene function, the selected DEIRs with a maximal expression of > 2 TPM (43) were chosen for the survival analysis. To reduce systematic error across different data sets, the gene expression data (TPM value) was normalized to [0,1] within every DEIR. Furthermore, the DEIRs that were not significantly correlated with the OS (P ≥ 0.05) in the univariate Cox analysis were filtered out in the 2 data sets. The related calculation formulae and R codes are listed in the Supplemental Methods.

LASSO is an acknowledged method for regression with high-dimensional data (44). It has been extensively applied to the Cox proportional hazard regression model for survival analyses (45, 46). The survival-related DEIRs that were common to the 2 RNA-Seq data sets were subsequently analyzed using LASSO Cox regression to select the most powerful prognostic markers. We first randomly selected 70% of TCGA patients for training (n = 329) and the remaining 30% of TCGA patients for internal validation (n = 140). A formula that combined the relative expression of the DEIRs (Exp_i) and their respective LASSO coefficients (L_i) was constructed to calculate a risk score (RS) for each patient:

(Equation 1)

In addition to the 405 LGGs from the CGGA RNA-Seq data set, a further 118 LGGs from the CGGA microarray data set, 88 LGGs from GSE16011, and 136 LGGs from GSE61374 were used as external validation sets. Primary microarray data from GSE16011 and GSE61374 were processed in the same way as the microarray data from the CGGA (6). Multiple probe sets were mapped to a single gene using the median value of the signals. To reduce the systematic error in the different data sets, the gene expression data of microarray was normalized to [0,1] in GSE16011, GSE61374, and CGGA. The normalization formula and R codes are listed in the Supplemental Methods.

The risk score was calculated for each patient in the validation sets using the formula constructed in the training set. The predictive accuracy of the risk score was assessed using time-dependent ROC curve analysis in all the data sets (47).

Construction and assessment of prediction model. Multivariate Cox regression analysis began with the risk score and the following clinical factors: age, grade, IDH mutation, 1p/19q status, histological type, chemotherapy, and radiotherapy. Step-wise backward elimination was applied with AIC to select the best variables to be included in the model (48). A nomogram was then constructed based on the multivariate Cox regression. As age and risk score were continuous variables, 3-knot restricted cubic splines were used (49).

The calibration of the nomogram was assessed using calibration curves. Harrell’s C-index was calculated to assess the discrimination.

Statistics. All statistical analyses were performed using the R software version 3.5.2. Univariate and multivariate Cox regression analyses were used to evaluate the prognostic value of factors. Pearson’s χ² test was used to compare the distribution of patient characteristics. The Kaplan-Meier analysis with a 2-sided log-rank test was used to compare the OS of patients. Spearman’s correlation test was performed to assess the correlation. Statistical significance was set at P < 0.05 unless specified otherwise.

Study approval. This retrospective study was approved by the ethics committee of Beijing Tiantan Hospital, and informed consent was waived.