Identification of a novel cuproptosis-related gene signature and integrative analyses in patients with lower-grade gliomas

The gene expression data and corresponding clinical information of LGG samples were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/↗) and Chinese Glioma Genome Atlas (CGGA) (http://cgga.org.cn/↗) databases. A total of 529 samples with a gene expression profile, copy number variation (CNV), single-nucleotide variant (SNV), and relevant clinicopathological data were downloaded from TCGA-LGG. The fragments per kilobase million values of TCGA-LGG were transformed into transcripts per kilobase million. For all the included RNA-seq data, normalization and log2 transformation were performed. The loss and gain levels of copy-number changes have been identified using segmentation analysis and the GISTIC algorithm. The SNV data were further analyzed by R package “maftools” and visualized by R package “oncoplot”. We further downloaded the TPM- normalized GTEx RNAseq data of 1,152 normal human brain samples from the GTEx data portal (https://xenabrowser.net/datapages/↗). Two CGGA cohorts that contained 182 and 174 LGG samples (CGGA1, mRNAseq_325, RNA-seq; CGGA2, mRNA-array_301, Microarray) were obtained as external validation cohorts.

A total of 13 CRGs were extracted from the previous study (6), and the details of genes were shown in Supplementary Table 1. Based on the expression levels of CRGs, we performed consensus unsupervised clustering analysis to classify patients into distinct cuproptosis-related clusters (CRG clusters) using the R package “ConsensusClusterPlus” (22), with the parameters of reps = 1000 and pItem = 0.8. Principal component analysis was conducted to show the classification of CRG clusters. Then, we compared the OS probability of CRG clusters using the R package “survival” and “jskm”. A landmark time of 9 years was chosen. Chemoradiotherapy sensitivity–related genes (CRSGs) and immune checkpoint genes (ICGs) were further retrieved (23–28), and their expression levels between CRG clusters were analyzed.

R package “ESTIMATE” can calculate TME scores including the stromal score, immune score, and estimate score using gene expression profiles (29). TME scores for LGG patients were evaluated and compared between CRG clusters. A single sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the immune infiltration degree of immune cells in the LGG TME.

To investigate the difference of the biological function between CRG clusters, gene set variation analysis (GSVA) was performed with “c2.cp.kegg.v7.5.symbols” and “c5.go.bp.v7.5.symbols” using R package “GSVA”. R package “pheatmap” was applied to visualize the results. GSEA was performed by R package “clusterProfiler” to determine whether the prior-defined functional sets of genes differ significantly between CRG clusters with the hallmark gene set (“h.all.v7.2.symbols”) from the MSigDB database.

The differentially expressed genes (DEGs) between different CRG clusters were identified using R package “Limma”. The significance criteria for identifying DEGs was set as |log2 (FoldChange)| > 0.5 and adjusted P-value< 0.05. To explore the biological functions of CRG cluster-related DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted by applying the “clusterProfiler” package (30).

We performed univariate Cox regression analysis for CRG cluster–related DEGs to identify DEGs that were related to OS (OS-related DEGs). According to the expression levels of OS-related DEGs, we performed consensus unsupervised clustering analysis using R package “ConsensusClusterPlus”, with the parameters of reps = 1000 and pItem = 0.8. The TCGA-LGG patients were divided into distinct cuproptosis gene clusters, and the OS time was compared through Kaplan–Meier analysis.

Then, patients in the TCGA-LGG cohort were randomly divided into the training cohort and internal testing cohort at a ratio of 1:1 using R package “caret”. Based on OS-related DEGs, the least absolute shrinkage and selection operator (LASSO) Cox regression were performed to reduce the dimension of high-latitude data using R package “glmnet”. Ten-fold cross-validation was employed to avoid the overfitting problem and select the penalty parameter (λ) according to the minimum criteria. We conducted a multivariate Cox regression analysis to determine genes from candidate genes and further performed GSEA analyses based on a single gene expression, respectively. Next, we construct the cuproptosis-related predictive model in the training cohort. We calculated the CRG score for each sample using the following formula:, with Coef indicating the coefficient and Exp referring to the expression level of each CRG. C R G s c o r e C o e E x _ = × ∑ i = 1 n f i p i

The prognostic scoring system for LGG patients was established, and the median value of the predicted CRG scores was regarded as the cut-off. Then, patients were divided into high-risk (CRG score > median value) and low-risk (CRG score< median value) groups accordingly. R package “survival” and “survminer” was applied to compare the survival probability between the two groups via Kaplan–Meier analysis. The R package “timeROC” was employed to perform 1-, 3- and 5- year receiver operating characteristic (ROC) analysis and calculate the value of the area under the curve (AUC). The calibration plots were further conducted to better validate the advantage of the CRG_score. The expression levels of CRGs were analyzed between different CRG score groups. The internal testing cohort and two external cohorts CGGA1 and CGGA2 were employed to verify the cuproptosis-related prognostic model. The CRG score was calculated for LGG patients in each cohort, and samples were divided into different risk groups. Similarly, they were subjected to Kaplan–Meier analysis, ROC analysis, and calibration analysis. In CGGA1 and CGGA2 cohorts, patients were also stratified into four risk-treatment subgroups by the CRG score and treatment with temozolomide (TMZ) or radiotherapy. In addition, survival analyses among risk-treatment subgroups were conducted.

In order to identify the gene sets of statistical differences between high- and low-risk groups, GSEA was performed. The annotated gene sets “h.all.v7.2.symbols” and “c5.bp.v7.2.symbols” from the MSigDB database were adopted in our analysis. The enrichments of gene sets with an adjusted P-value<0.05 were regarded to be significant. We employed R package “ESTIMATE” to evaluate the TME score levels between high- and low-risk groups. Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) is a developed algorithm that uses a set of reference gene expression matrices to evaluate 22 immune cell type proportions from bulk tumor sample expression data based on the principle of linear support vector regression (31). We processed the TCGA-LGG RNA-Seq data (TPM normalized) to calculate the relevant abundance of immune cells. We analyzed the Spearman correlation between the abundance of infiltrating immune cells and the CRG score. Furthermore, we downloaded the RNAss file named “StemnessScores_RNAexp_20170127.2.tsv”. The tumor stem cell characteristics were extracted from the transcriptome and epigenetics of the samples and then used to evaluate the stem cell-like features of tumors. We performed a correlation analysis the investigate the association between the CRG score and Cancer Stem Cell (CSC) index.

TMB and ICGs were associated with patients’ response rate to immunotherapy. We extracted the mutation annotation format (MAF) from the TCGA database with the “maftools” R package to identify the mutational landscape of LGG patients between different CRG score groups. The TMB score was also calculated for each LGG patient in the entire TCGA cohort. Then, we evaluated the correlations of ICGs with the CRG score and five genes in the cuproptosis-related prognostic model using Spearman’s rank correlation coefficient.

The immunotherapy response for LGG patients was estimated through the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu↗), which can help doctors select patients who are more suitable for immunotherapy (32). Furthermore, GSE126044↗, GSE78220↗, checkmate cohort (33), and IMvigor210 cohort (34) were used to validate the predictive ability of the cuproptosis-related prognostic model in immunotherapy response.

We compared the expression levels of CRSGs between different risk groups and performed the correlation between CRG scores and gene expression levels. The calcPhenotype function of the “oncoPredict” R package was applied to estimate drug sensitivity scores for common drugs in the LGG therapy regimen including TMZ, procarbazine, teniposide, and vincristine in the TCGA cohort. The lower-imputed drug sensitivity represents more sensitivity to the drug.

We obtained the clinical characteristics including the age, grade, and IDH mutation status of LGG patients in the entire TCGA cohort and two CGGA cohorts. In combination with the CRG score, these variables were analyzed in univariate and multivariable Cox regression analyses.

To individualize the predicted LGG patients’ survival probability, we developed a nomogram using clinical characteristics and the CRG score. The R package “rms” and “regplot” were employed. Time-dependent ROC analysis was conducted to assess the predictive accuracy for 1-, 3-, and 5-year survival probability. The calibration plots were applied to compare model-predicted probability with observed outcomes in the TCGA-LGG cohort and two validation cohorts (CGGA1 and CGGA2).

To further validate the potential roles of signature genes in LGG, five paired LGGs and matched normal adjacent tissue samples were collected from the SYSUCC. Ethical approval was confirmed by the ethical committee of the hospital. Associated clinicopathological features were further confirmed as listed in. Tissue specimens were frozen in liquid nitrogen and stored at –80°C until used. 1

Total RNA was extracted with a TRIzol Reagent (ThermoFisher: #15596018), and the concentration was calculated by the A260/A280 ratio. The PrimeScript RT reagent kit (EZBioscience: #A0010CGQ) and SYBR-Green PCR reagent (EZBioscience: #A0012-R2-L) were used to perform cDNA synthetization and further conduct RT-qPCR based on the LightCycler ^® 480 System (Roche). The housekeeping gene GAPDH was used as an endogenous control. The 2^−ΔΔCT cycle threshold method was used to calculate the relative expression. Supplementary Table 3 lists the primers used in this study.

The protein expression levels of crucial CRGs were confirmed by Western blotting. Tissues were treated with RIPA lysis buffer (Fdbio: #FD009) containing phosphatase and protease inhibitors. The BCA protein detection kit (ThermoFisher: #23227) was applied to detect the protein concentration. Equivalent protein was then separated by 10% Tris-Tricine SDS-PAGE and transferred onto polyvinylidene fluoride (PVDF) membranes. After blocking with skimmed milk in TBST for 2 h, the membrane was further probed using antibodies against GAPDH (Proteintech: #60004-I-Ig), C21orf62 (Signalway Antibody, SAB: #C08364H), DRAXIN (Proteintech: #26342-1-AP), ITPRID2 (Proteintech: #14157-1-AP), MAP3K1 (Proteintech: # 19970-1-AP), and MOXD1 (Bioss: bs-17733R) overnight at 4°C. The membranes were subsequently washed with Tris-buffered saline containing Tween and then incubated with an HRP‐conjugated anti-rabbit antibody at 37°C for 1 h. Finally, the bands on the membranes were observed with a ChemiDoc™ Imaging System.

All statistical analyses were performed using R software (Version 4.1.2). Correlation coefficients were evaluated by Spearman analysis. To compare variables between two groups, we employed the independent sample t-tests for normally distributed continuous variables, and Mann–Whitney U tests for non-normally distributed continuous variables. One-way ANOVA and Kruskal–Wallis tests were used to perform the difference comparisons of three or more groups. The survival analysis was conducted via the Kaplan–Meier method, and log-rank tests were employed to identify the significance of differences. The statistical significance was defined as P< 0.05.

TThe workflow of the study is outlined in Figure 1. A total of 13CRGs were included in our study The genetic mutation landscape in LGG patients is shown in Figures 2A–D. Of the 506 LGG patients in the TCGA cohort, 488 (96.44%) had genetic mutations and IDH1 had the highest mutation frequency (77%), followed by TP53, ATRX, CIC, and TTN. However, only seven samples had genetic mutations in CRGs (Supplementary Figure 1). We investigated the frequencies of the CNVs of 13 CRGs in LGG. DLD exhibited the highest amplification frequency, while ATP7B and DLST had a widespread frequency of CNV loss (Figure 2E). Figure 2F shows the location of the CNV alterations of 13 CRGs on 23 chromosomes.

We then explored the expression levels, molecular interactions, and prognostic values of 13 CRGs. Twelve CRGs were upregulated in tumor samples including FDX1, LIPT1, LIAS, DLD, DBT1, GCSH, DLST, DLAT, PDHA1, PDHB, SLC31A1, and ATP7A (P< 0.001), whereas only ATP7B was downregulated (P< 0.001) (Figure 2G). Figure 2H exhibits the molecular interactions between CRGs. Nine prognostic CRGs were identified by Kaplan–Meier analysis and univariate Cox regression analysis (Supplementary Figure 2). The result of multivariate Cox regression analysis further revealed that three prognostic CRGs (FDX1, GCSH, and ATP7B) were independent prognostic factors (Table 1).

In order to investigate the expression features and potential biological characteristics of CRGs in LGG, a consensus clustering algorithm was utilized to classify LGG patients in the TCGA cohort. Based on the expression of 13 CRGs, patients were categorized into CRG cluster A (n=219) and CRG cluster B (n=292) (–). The PCA plot demonstrated an obvious different distribution between CRG clusters (). 1 1 1

Figure 3A shows the different expressions of CRGs and clinicopathological characteristics between CRG cluster A and B. CRG cluster A was preferentially associated with higher expression levels of CRGs, higher grade (G3), and more death events. As the Kaplan–Meier survival curves crossed, we employed landmark analysis to compare the difference between CRG clusters (Figure 3B). The result of landmark analysis showed a longer OS in LGG patients in CRG cluster B within 9 years (P = 0.003). Nevertheless, no significant difference was found in the survival probability beyond 9 years (P = 0.149).

Then, we explored the correlation of CRG clusters with ICGs, CRSGs, and the TME. We found that CRG cluster A was associated with a higher expression of ICGs (Figure 3C). In addition, CRSGs, including AKR1C1, EGFR, EZH2, HOXA9, HGMT, SOX2, and TBX5, were differentially expressed between two CRG clusters (Figure 3D). To explore the potential function of CRGs in the immune infiltration of LGG, we compared the TME score and the relevant abundance of immune cells between two CRG clusters using “ESTIMATE” and “ssGSEA” algorithms. Patients in CRG cluster B had higher immune scores than those in CRG cluster A (Figure 3E). We observed that the immune infiltration levels of the activated CD4+ T cell, gamma delta T cell, and Type 2 T helper cell were significantly higher in CRG cluster A than those in CRG cluster B, while the CD56 dim natural killer cell and monocyte had significantly lower infiltration levels in CRG cluster A than those in CRG cluster B (Figure 3F).

To further explore the functional annotation between CRG cluster A and B, GSVA and GSEA were performed. The results of the GSVA of the gene ontology biological process (GOBP) showed that CRG cluster A was significantly enriched in transportation-related processes including Golgi vesicle transport, nuclear pore organization, and vesicle targeting (Figure 4A). In addition, the GSVA of KEGG terms showed that CRG cluster A was abundant in metabolism-related pathways (citrate cycle TCA cycle, glycosylphosphatidylinositol GPI anchor biosynthesis), cancer-related pathways (small cell lung cancer, endometrial cancer), cell cycle–related pathways (cell cycle), and genomic stability–related pathways (mismatch repair, nucleotide excision repair) (Figure 4B). GSEA indicated that CRG cluster A was predominantly associated with the cell cycle, tumorigenesis, and metastasis including the G2M checkpoint hallmark, E2F target hallmark, epithelial mesenchymal transition hallmark, and angiogenesis hallmark (Figures 4C–E) (Supplementary Tables 4–5).

A total of 1,966 CRG cluster–related DEGs were identified between CRG clusters using R package “Limma”. Consistent with GSVA and GSEA, the result of GO and KEGG showed that these DEGs were associated with the cell cycle, genomic stability, and cancer, which revealed that cuproptosis plays a significant role in tumorigenesis and metastasis (Figures 4F–H) (Supplementary Table 6).

Subsequently, the prognostic values of the above 1,966 CRG cluster-related DEGs were assessed by univariate Cox regression analysis and a total of 1,424 genes associated with OS (OS-related DEGs) were screened out (P< 0.05). Based on the expression levels of OS-related DEGs, a consensus clustering algorithm was employed and LGG patients were classified into three gene clusters, termed gene cluster A (n=215), gene cluster B (n=88), and gene cluster C (n=208) (Supplementary Figures 4A–H). PCA analysis revealed that gene clusters could be identified clearly (Supplementary Figure 4I). The expression profiles and clinical information of OS-DEGs in different gene clusters are shown in Figure 5A. We found that gene cluster B was correlated with high gene expression levels, an advanced grade (G3), and more death events. Patients in gene cluster B were proven to be related to worse prognosis than those in gene cluster A and C by Kaplan–Meier analysis (P< 0.001) (Figure 5B). Furthermore, CRGs were differentially expressed among three gene clusters, with the highest expression levels in gene cluster B and the lowest expression levels in gene cluster C. (Figure 5C).

We performed LASSO and multivariate Cox regression analysis for 1,424 OS-related DEGs to establish the cuproptosis-related prognostic model. Patients in the TCGA cohort were divided into training and test cohorts at a ratio of 1:1. In the TCGA training cohort, 11 genes were obtained, followed by LASSO Cox regression analysis, and were subjected to multivariate Cox regression analysis (Supplementary Figure 5). Ultimately, five key genes remained, including Chromosome 21 Open Reading Frame 62 (C21orf62), Dorsal Inhibitory Axon Guidance Protein (DRAXIN), ITPR Interacting Domain Containing 2 (ITPRID2), Mitogen-Activated Protein Kinase Kinase Kinase 1 (MAP3K1) and Monooxygenase DBH Like 1 (MOXD1). Further GSEA analyses based on a single gene expression indicated that B-cell-mediated immunity, complement cascade, interferon-gamma response, chromosome segregation, mitotic spindle checkpoint, G2M checkpoint, T-cell receptor complex, CD22-mediated BCR regulation, epithelial–mesenchymal transition, cell cycle checkpoints, voltage-gated potassium channels, and other biological processes related to immunity or cancer were significantly enriched (Supplementary Figure 6). Next, we constructed a five-gene cuproptosis-related prognostic model and the associated CRG score can be calculated as follows: CRG_score= 0.164859* C21orf62 + 0.293187* DRAXIN + 0.882099* ITPRID2 + 0.625577* MAP3K1 + 0.256801* MOXD1. We also explored the relationship between CRG clusters, gene clusters, and CRG scores. We found that CRG scores in CRG cluster A were significantly higher than that in CRG cluster B, and the expression of CRGs were upregulated in the high-risk group, which suggested that a high CRG score and the high expression of CRGs were associated with tumorigenesis and metastasis (Figures 5D, E). Meanwhile, the rank order of the CRG score in gene clusters was B > A > C (Figure 5F). The Sankey diagram showed subgroup distributions in groups with different CRG scores and survival outcomes (Figure 5G).

The heatmap showed the different expressions of C21orf62, DRAXIN, ITPRID2, MAP3K1, and MOXD1 between the high- and low-risk groups in TCGA training and test cohorts (Figure 6A). In addition, the risk plot of the CRG score revealed that patients with a high CRG score were related to more cases of death and shorter survival time (Figures 6B, C). Kaplan–Meier analysis showed that patients in the high-risk group had a worse OS than low-risk patients in both TCGA training and testing cohorts (P< 0.001) (Figures 6D, E). The 1-, 3-, and 5-year survival probability of the CRG score was represented by the AUC values of 0.876, 0.863, and 0.844, respectively, in the TCGA training cohort (Figure 6F). In the TCGA testing cohort, AUC values for predicting 1-, 3-, and 5-year OS were 0.852, 0.856, and 0.805. Meanwhile, the calibration curve also manifested a satisfactory agreement between predictive and observational values at the probabilities of 3- and 5-year survival (Supplementary Figures 7A, B).

To further verify the prognostic performance of the model, we applied two external validation cohorts (CGGA1 and CGGA2) (Supplementary Table 7). Patients were also classified into high- and low-risk groups according to CRG scores. Kaplan–Meier analysis showed a significantly better prognosis in the low-risk group compared to that in the high-risk group (Figures 7A, B). Meanwhile, the model also demonstrated a high AUC value in external validation cohorts (Figure 7C, D), and the calibration curve also exhibited a satisfactory agreement between predictive and observational values at the probabilities of 3- and 5-year survival (Supplementary Figures 7C, D).

Then, patients were stratified into four groups by the CRG score and treatment strategies. Kaplan–Meier analysis showed that patients in low-risk and no-radiotherapy groups had the best prognosis (Figures 7E, F). However, the predictive ability of the model was not impacted by TMZ treatment. Whether patients received TMZ or not, the low-risk group always showed a strong survival advantage (Figures 7G, H).

GSEA was conducted to explore the potential biological functions between high- and low-risk groups (Figures 8A, B) (Supplementary Table 8). We observed that a high CRG score was mainly related to the cell cycle (negative regulation of metaphase–anaphase transition of the cell cycle, G2M checkpoint hallmark, E2F target hallmark), tumor progression (epithelial–mesenchymal transition hallmark, angiogenesis hallmark), and immunity (interferon gamma–mediated signaling pathway, T-cell activation via T-cell receptor contact with antigen bound to the MHC molecule on antigen-presenting cell, inflammatory response hallmark). TME scores including the stromal score, immune score, and ESTIMATE score were significantly higher in the high-risk group (Figure 8C). We assessed the correlation of immune infiltration with the five genes cuproptosis-related prognostic model (Figure 8D). Macrophage M0 and Macrophage M1 were positively correlated with the CRG score, while mast cells activated, monocytes, neutrophils, and NK cells activated were negatively correlated with the CRG score (Figures 8E–K).

Cancer stem cells are a group of cells with the features of self-renewing, being multipotent, and tumor-initiating, which can drive the growth and recurrence of tumors and are resistant to many current treatments. The CRG score and CSC index were synthesized to evaluate their relationship. A mild but significant negative correlation (R = -0.24, P = 4.2e-08) was detected (Figure 8L).

Previous studies demonstrated that high TMB scores are associated with increased treatment response to immunotherapy. Here, patients in the high-risk group (64%) had a markedly lower mutation incidence of IDH1 than those patients in the low-risk group (91%) (Figures 9A, B). The TMB score was significantly higher in the high-risk group (P = 4.8e-0.8), and the CRG score was positively correlated with the TMB score (R = 0.31, P = 2.8e-12) (Figures 9C, D). Considering that the expression levels of ICGs have been reported to correlate with the clinical benefit of checkpoint blockade immunotherapy (35, 36), we further explored the relationship between the CRG score and ICGs. We found that most ICGs were significantly correlated with the five genes in the model (Figure 9E). The expression levels of 61 ICGs including CD276, BTN2A2, PDCD1LG2, CD274, and CD40LG increased with the increasing CRG score (Figure 9F). Only BTNL9 was negatively correlated with the CRG score (R = -0.24, P = 2.6e−08) (Supplementary Table 9). Then, we found that patients with a high CRG score and high expression levels of CD276, CD274, BTN2A2, PDCD1LG2, and CD40LG were associated with poor OS than others (Supplementary Figure 8).

The TIDE algorithm was applied to estimated immunotherapy response for LGG patients based on the transcriptomic data. The results showed that the TIDE score of the high-risk group was significantly higher than that of the low-risk group, indicating that the patients in the low-risk group could benefit more from immunotherapy (Figure 10A). Patients in the low-risk group were associated with a higher dysfunction score and lower exclusion score (Figures 10B, C). Moreover, patients were classified into no responders and responders by the TIDE algorithm. We found that immunotherapy responders were correlated with a lower CRG score (Figures 10D, E). Patients with a combination of a high CRG score and high TIDE score group showed an association with the worst prognosis (Figure 10F). Then, we validated the predictive value of the cuproptosis-related prognostic model in immunotherapy response. The clinical response to immunotherapy in patients with non-small cell lung cancer or metastatic melanoma in the low-risk group compared to those in the high-risk group were confirmed (Figures 11A, B). However, no differences were identified in renal cell carcinoma patients and urothelial cell carcinoma patients (Figures 11C–E). Figure 11F shows the subtype distributions in groups with different CRG scores and immunotherapy response in the IMvigor210 cohort. Furthermore, in the IMvigor210 cohort, the CRG scores were different among different immune phenotypes, tumor cell (TC) levels, and immune cell (IC) levels (Figures 11G–I).

Additionally, we observed a strong correlation between the expression levels of CRSGs and the expression levels of C21orf62, DRAXIN, ITPRID2, MAP3K1, and MOXD1 (Figure 12A). As the CRG score increased, the expression levels of CPZ, EGFR, EZH2, and HOXA9 increased, but the expression levels of AKR1C1 decreased, which revealed a potential association between the CRG score and chemoradiotherapy (Figures 12B–G). To explore the values of the CRG score as a biomarker to predict the chemotherapeutic response in LGG patients, we performed drug sensitivity analysis using the “oncoPredict” R package. A higher imputed sensitivity score represented lower sensitivity to the drug. We found patients in the high-risk group had higher imputed sensitivity scores of TMZ and procarbazine, while the imputed sensitivity scores of teniposide and vincristine were lower in patients in the high-risk group (Figures 12H–L).

To determine if the CRG score could be used as an independent prognostic predictor for OS, we combined clinical characteristics and the CRG score to conduct univariate and multivariate Cox regression analyses. In the TCGA cohort, age, grade, IDH mutation status, and CRG score demonstrated significant differences (Figures 13A, B). The prognostic value of the CRG score was also verified in external CGGA cohorts, and the results showed that CRG score is an independent prognostic predictor for LGG patients (Supplementary Figure 9).

Then, we combined the CRG score and clinical chrematistics (age, grade, IDH mutation status) to establish a nomogram in the TCGA cohort, which exhibited a quantitative method to generate personalized predictions for LGG patients (Figure 13C). The AUC values of the model were estimated, and calibration analyses were performed to assess the predictive ability and accuracy for prognosis. Figure 13D showed that 1-, 3- and 5-year AUC values of the model in the TCGA cohort were 0.893, 0.887, and 0.828, respectively. The model also had good predictive ability in the CGGA1 cohort (AUC > 0.75) and CGGA2 cohort (AUC > 0.7) (Figures 13F, H). Subsequently, the calibration plots demonstrated good agreement between model-predicted probability and the observed outcomes (Figures 13E, G, I).

To verify the expression level of signature genes in LGG, we collected five paired cancer- and adjacent normal tissues from SYSUCC. As shown in Figures 14A–D, qRT-PCR showed that the expression of the DRAXIN, ITPRID2, and MAP3K1 were significantly upregulated while MOXD1 was downregulated in tumor samples. Nevertheless, there was no significant difference in the expression of C21orf62 (Figure 14E). Further investigation indicated heterogeneity in the expression of C21orf62. Briefly, it was upregulated in two patients but downregulated in the remaining patients (the relative expression was 5.33, 1.68, 0.36, 0.07, and 0.04, respectively). The results of WB also demonstrated the heterogeneity in C21orf62, MAP3K1, and MOXD1 but consistently a significant upregulation in DRAXIN and ITPRID2 at the protein level (Figure 14F).