The consensus clustering results showed that the best clustering was achieved when k = 2. All glioma samples were classified into cluster 1 (n = 477) and cluster 2 (n = 217) (Figures 1A,B). In addition, the overall survival (OS) of cluster 1 patients was worse (p < 0.05). From the above, it is clear that NAD+ metabolism-related genes (NMRGs) can classify glioma patients into two subtypes with statistical difference in OS, indicating that NMRGs can affect the OS of glioma patients (Figure 1C).

5,998 DEGs were obtained between glioma and normal samples (Figure 2A and Supplementary Tables S1,S2). 11 NMR-DEGs were obtained after crossover analysis (Figure 2B and Table 1), and they were enriched in metabolic processes such as NAD+ metabolism, pyridine nucleotide metabolism, nucleotide metabolism, nicotinate and nicotinamide metabolism, as well as NAD+ biosynthesis process, pyridine nucleotide biosynthesis, and nucleotide biosynthesis. And NMR-DEGs were associated with NAD+ nucleotidase, circulating ADP-ribose-gluconeogenesis, NAD(P)+ nucleosidase activity, NAD+-dependent protein deacetylase activity, NAD+ binding transferase activity, NAD+ nucleosidase activity, NAD + ADP ribosyl transferase activity, transferase activity-transferring pentanediyl group and other pathways (Figures 2C,D).

11 NMRGs were obtained (Figure 3A). When lambda min = 0.0313, PARP9, BST1, NMNAT2, and CD38 were obtained (Figure 3B and Table 2). Then the glioma cases were classified into a high-risk group (n = 347) and a low-risk group (n = 347) (Median Risk score = 0.8125), and cases in the high-risk group had a lower OS (Figures 4A,B). The area under curves (AUCs) of the ROC curves in the training cohort were all greater than 0.8 (Figure 4C). In the validation set, the OS was worse in the high-risk group, and the AUCs of the 1, 3, 5 years in the ROC curves were >0.7, consistent with the results of TCGA (Figures 4D–F).

In the TCGA-glioma dataset, risk Scores were significantly correlated with age, Grade, survival status, IDH status, MGMT status, ATRX status and BCR status. And risk Scores for age, Grade, BCR status, IDH status, MGMT status, ATRX status, and BCR status were significantly different (Table 3). In the CGGA-glioma dataset, risk Scores for gender were significantly different (Table 4). In addition, the risk Scores for age (>50 vs. ≤50), Grade (G2/G3/G4), BCR status (NCH vs. IGC), IDH status (WT vs. Mutant), MGMT status (Unmethylated vs. Methylated) and ATRX status (WT vs. Mutant) were significantly different (Figure 5). Moreover, the risk Score was higher in tumor tissue (Figure 6A). The ROC curves demonstrated that the risk Score could significantly distinguish between tumor and normal samples (AUC = 0.81, Figure 6B), and it had a strong ability to distinguish between high-grade glioma (G3 + G4) and low-grade glioma (G2) (AUC = 0.802, Figure 6C). Besides, the risk Score could also distinguish between GBM and non-GBM (AUC = 0.922, Figure 6D).

The results showed that age, Grade, IDH status, ATRX status, BCR Status, and risk Score in Cox analysis were all p < 0.05 (Figures 7A,B). The nomogram in TCGA glioma samples was shown in Figure 7C, and the C-index was 0.8617, which indicated that the model had a better prediction effect. The correction curve of nomogram was shown in Figure 7D, and the closer the slope was to 1, the more accurate the prediction was. The accuracy of this risk model in predicting the OS of cases was high, indicating that the constructed prediction model could be a valid model.

309 pathways were significantly different, 40 pathways were connected with risk Score, 12 of which were positively connected with risk Score and 28 were negatively connected with risk score (Supplementary Table S3). The heat map of risk score and KEGG pathway was shown in Figure 8. As could be seen from the figure, all the pathways were positively connected with risk Score except for the Lysosome pathway, which was negatively connected with risk Score. Moreover, neuroactive ligand-receptor interactions, aldosterone synthesis and secretion, GnRH secretion, GABAergic synapses, adrenergic signaling pathway in cardiac myocytes, insulin secretion, and glutamatergic synapses were significantly correlated with risk Score.

Immune score, stromal score and ESTIMATE score were higher in high-risk group (p < 0.05) (Figure 9A). The CIBERSORT algorithm calculated the proportion of each immune cell and the corresponding statistical values. After excluding samples with p > 0.05, 37 samples remained in the high-risk group and 3 samples remained in the low-risk group. The box line plots and stacked bar graphs plotted for the abundance of immune cells were shown in Figures 9B,C. Mast cells activated, NK cells activated, and B cells naive were less in the high-risk group (Figure 10A). All 24 HLA family genes were higher in the high-risk group (Figure 10B). Except for TIGIT, all the immune checkpoints were higher in the high-risk group (Figure 10C).

CancerSEA can provide insight into the model genes in individual glioma cell. NMNAT2 was positively connected with stemness and proliferation, and negatively connected with hypoxia (Figure 11A). PARP9 was significantly positively connected with EMT, cell metastasis, cell differentiation and negatively connected with DNA damage (Figure 11B). CD38 was negatively connected with both cell invasion and apoptosis (Figure 11C). BST1 was negatively connected with cell invasion (Figure 11D).

We acquired the expression data and clinical data of glioma from the TCGA and CGGA databases (Table 5). The TCGA database contains 697 glioma and 5 normal samples, of which 694 glioma samples had survival information. The CGGA 693 and CGGA 325 datasets contain 657 and 313 glioma samples with survival information, respectively. In addition, NMRGs were acquired from the Kyoto Encyclopedia of Genes and Genomes and the Reactome databases (Table 5). Finally, 49 NMRGs were obtained by combined genes obtained from the two databases (47).

Consensus clustering is a common method for classifying cancer subtypes (48). This study used “ConsensusClusterPlus” package (version 1.54.0) (49) to execute consensus clustering based on 42 NMRGs (Seven genes were not expressed in the Glioma samples) in the TCGA-glioma dataset. Then the survival analysis of different subtypes was performed by “survival” package (version 3.2–11).

The differentially expressed genes (DEGs) between glioma and normal samples were acquired by “limma” package in TCGA-glioma database (version 3.46.0) (50). Next, the DEGs were intersected with the NMRGs using the “VennDiagram” package (version 1.6.20) (51) to obtain NMR-DEGs. Then, functional enrichment analysis of NMR-DEGs was performed based on gene ontology (GO) and KEGG databases using the “clusterProfiler” package (version 3.18.0) (52). The enrichment results were also visualized by plotting bar graphs using the “enrichplot” package (version 1.10.2).

In this study, the 694 glioma cases from the TCGA database were used as the training cohort and the 657 cases from the CGGA database were used as the validation cohort to construct and validate the risk model. First, we extracted the expression data of NMR-DEGs in the training cohort for univariate Cox analysis. Then the factors with p < 0.05 were subjected to the least shrinkage and selection operator (Lasso) regression analysis. The risk Scores of glioma cases were counted using the following formula.

The median value of risk Score was used as the boundary to classify the glioma cases into high and low risk groups. Then risk curves were plotted, and survival curves and receiver operating characteristic (ROC) curves were drawn by “survminer” (version 0.4.8) and “survivalROC” package (version 1.0.3) respectively. Also, we plotted risk curves, survival curves, and ROC curves in the validation set.

To further investigate the relationship between clinicopathological characteristics and risk model, age, gender, survival status, Grade, IDH status, MGMT status, ATRX status, DAXX status, and BCR status were correlated with risk Scores in the training and validation cohorts. First, the number of patients with different clinical subtypes was compared between the high and low risk groups, and then the risk Scores were compared in the training set with different clinical information using rank sum tests. In addition, according to the abundance of model genes and their risk coefficient, the risk Score of normal samples was calculated, and the risk Score between tumor and normal samples was compared by the rank-sum test. Finally, the ROC curves were used to investigate the ability of the risk score to distinguish tumor tissue from normal tissue, and to distinguish between different Grades.

To discovery the independent prognosis of clinicopathological characteristics, risk scores and clinicopathological factors such as age, gender, Grade, IDH status, MGMT status, ATRX status, DAXX status, and BCR status were included in the risk model for Cox analyses in the training set. Next, a nomogram plot predicting survival of glioma patients was constructed. Finally, a calibration curve was plotted.

We performed GSVA enrichment analysis based on KEGG pathway for all genes using “GSVA” package to find the functional difference between high and low risk groups (53). First, we obtained the enrichment scores in each pathway for the samples, and then performed a differential analysis of the pathways by the “limma” package with the screening condition of adjust p < 0.05. Next, we executed a correlation analysis of the differential pathways with the risk scores by setting the screening condition of |cor| > 0.6 and p < 0.01.

The immune status difference between high and low risk groups were executed by ESTIMATE analysis. The Immune Score, Stromal Score, and the combined ESTIMATE score of TCGA glioma samples were obtained by “ESTIMATE” package (version 1.0.13). This study compared the three scores between the two groups by rank sum test. Then the proportion of 22 immune cells was calculated using the Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm (version 1.03) in the TCGA-glioma dataset. The proportion of each immune cell in each sample were gotten by CIBERSORT algorithm. The samples with p > 0.05 were excluded, and the proportions of immune cells in the remaining samples were displayed in box plots and bar charts. In addition, the immune cells between high and low risk groups were compared by the rank sum test and visualized by plotting box line plots using the “ggplot2” and “ggpubr” packages. Finally, the expression of 24 HLA related genes and nine immune checkpoints (LAG3, ICOS, TIGIT, CD274, PDCD1, IDO1, CD27, PDCD1LG2, and HAVCR2) were compared by rank sum test respectively.

In this study, 14 functions (angiogenesis, cell cycle, apoptosis, DNA damage, cell differentiation, DNA repair, EMT, cellular hypoxia, cancer cell invasion, metastasis, inflammation onset, proliferation, cell resting, and stem cell properties) of model genes were predicted using CancerSEA. Then the functional status of the model genes in glioma were queried and the significantly relevant functions were displayed.