Development of prognostic indicator based on NAD+ metabolism related genes in glioma

Background: Studies have shown that Nicotinamide adenine dinucleotide (NAD+) metabolism can promote the occurrence and development of glioma. However, the specific effects and mechanisms of NAD+ metabolism in glioma are unclear and there were no systematic researches about NAD+ metabolism related genes to predict the survival of patients with glioma.

Methods: The research was performed based on expression data of glioma cases in the Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases. Firstly, TCGA-glioma cases were classified into different subtypes based on 49 NAD+ metabolism-related genes (NMRGs) by consensus clustering. NAD+ metabolism-related differentially expressed genes (NMR-DEGs) were gotten by intersecting the 49 NMRGs and differentially expressed genes (DEGs) between normal and glioma samples. Then a risk model was built by Cox analysis and the least shrinkage and selection operator (LASSO) regression analysis. The validity of the model was verified by survival curves and receiver operating characteristic (ROC) curves. In addition, independent prognostic analysis of the risk model was performed by Cox analysis. Then, we also identified different immune cells, HLA family genes and immune checkpoints between high and low risk groups. Finally, the functions of model genes at single-cell level were also explored.

Results: Consensus clustering classified glioma patients into two subtypes, and the overall survival (OS) of the two subtypes differed. A total of 11 NAD+ metabolism-related differentially expressed genes (NMR-DEGs) were screened by overlapping 5,995 differentially expressed genes (DEGs) and 49 NAD+ metabolism-related genes (NMRGs). Next, four model genes, PARP9, BST1, NMNAT2, and CD38, were obtained by Cox regression and least absolute shrinkage and selection operator (Lasso) regression analyses and to construct a risk model. The OS of high-risk group was lower. And the area under curves (AUCs) of Receiver operating characteristic (ROC) curves were >0.7 at 1, 3, and 5 years. Cox analysis showed that age, grade G3, grade G4, IDH status, ATRX status, BCR status, and risk Scores were reliable independent prognostic factors. In addition, three different immune cells, Mast cells activated, NK cells activated and B cells naive, 24 different HLA family genes, such as HLA-DPA1 and HLA-H, and 8 different immune checkpoints, such as ICOS, LAG3, and CD274, were found between the high and low risk groups. The model genes were significantly relevant with proliferation, cell differentiation, and apoptosis.

Conclusion: The four genes, PARP9, BST1, NMNAT2, and CD38, might be important molecular biomarkers and therapeutic targets for glioma patients.

Introduction

Glioma originates mainly from glial cells and is the most common primary brain tumor. It involves a broad category of central nervous system tumors, including astrocytoma, oligodendroglioma, and glioblastoma. Currently, the standard clinical strategy for glioma includes maximum surgical resection followed by radiotherapy, and temozolomide (TMZ) chemotherapy (1). However, the overall treatment effect was not ideal. Even with standardized treatment, the median survival time of glioblastoma multiforme (GBM) was only about 15 months, and the 5-year survival rate was less than 10% (2). Thus, considering the limited treatment strategies for glioma, there is an urgent need to develop reliable prognostic biomarkers and therapeutic targets.

Nicotinamide adenine dinucleotide (NAD+) is one of the most important co-enzymes (or co-factors) in oxidation-redox reactions, and NAD+ is also the core of energy metabolism (3, 4). Playing a key role in energy transduction and cell signal transduction, NAD+ can be transformed into NADP, NAADP and cADPR (4, 5). Moreover, NAD+ degradation products, such as nicotinamide and N-methyl nicotinamide, have also been considered as key regulators of energy metabolism, epigenetics and disease status (6–8). NAD+ pathway metabolites can also serve as substrates for a group of diverse enzymes (9, 10), including SARM1, ARTs, PARPs, CD38, sirtuins, and RNA polymerases, which were involved in several aspects of cellular homeostasis. Previous studies mainly focused on sugar metabolism, lipid metabolism, amino acid metabolism, nucleotide metabolism and energy metabolism. Interestingly, an increasing number of studies have revealed that NAD+ metabolism was closely related to the pathogenesis of many tumors (11, 12). For example, in GBM patients, the expression of NAMPT was associated with poor prognosis (13). More than that, NAD+ metabolism was reported to play an important role in tumor immunity (11, 14–16). But the mechanisms of NAD+ metabolism in glioma were unclear and there were no systematic researches about NAD+ metabolism related genes to predict the survival of patients with glioma. Different from previous researches, our study aimed to explore the prognostic value of NAD metabolism related genes in glioma by bioinformatics analysis.

Therefore, we used the data downloaded from the Cancer Genome Atlas (TCGA) in this study to perform bioinformatics analysis based on NAD+ metabolism related genes and then build a new risk model to predict the prognosis of glioma patients. Then we validate the model using Chinese Glioma Genome Atlas (CGGA). This model could find prognostic molecular markers and potential therapeutic targets for glioma patients, which provided an important reference value on selection of treatment strategies for clinicians, and also provided new ideas for the future basic research of glioma.

Results

Classification of glioma subtypes

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).

FIGURE 1

Figure 1. NMRGs classify glioma patients into different two subtypes. (A1,A2) Cumulative Distribution Function (CDF) plot for consistent clustering (B) consistency clustering heat map. The colours indicate the ease of clustering. (C) Survival curves for cluster1 and cluster2.

Identification of NAD+ metabolism-related differentially expressed genes (NMR-DEGs)

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).

FIGURE 2

Figure 2. NAD+ metabolism-related differentially expressed genes were selected from glioma database. (A) Volcano map of differentially expressed gene analysis (DEGs) between glioma and normal samples. The horizontal coordinate log2FC indicates the difference multiplicity and the vertical coordinate indicates the confidence level - log10(adj. Pvalue). Each point in the graph represents a gene, with blue and red points representing significantly differentially expressed genes. Red dots indicate that their gene expression was up-regulated in the glioma samples and blue dots indicate that the gene was down-regulated in the glioma samples. (B) Venn diagram of differentially expressed genes with NMRGs. (C,D) Bar graphs of Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis for NMR-DEGs.

TABLE 1

Table 1. Differentially expressed genes between glioma and normal samples.

Building and validating of risk model

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).

FIGURE 3

Figure 3. Building of the risk model (A) forest plot for univariate cox regression analysis. (B1,B2) Plot of gene coefficients for Lasso regressions. Error plots for 10-fold cross-validation.

FIGURE 4

Figure 4. Building and validating of a risk model (A) the risk curve, scatter diagram, and model gene expression heat map of the training cohort (B) survival curves of the high- and low- risk groups of the training cohort. (C) The ROC curves of glioma patients in training cohort for 1, 3, 5 years (D) the risk curve, scatter diagram and model gene expression heat map of the CGGA 693 cohort. (E) Survival curves demonstrating that glioma cases in the high- and low- risk groups of the CGGA 693 cohort. (F) The ROC curves of glioma patients in the CGGA 693 cohort. (G) The risk curve, scatter diagram, and model gene expression heat map of the CGGA 325 cohort. (H) Survival curves for glioma cases of the CGGA 325 cohort in the two risk groups. (I) The ROC curves of glioma patients in CGGA 325 cohort.

TABLE 2

Table 2. Coef of 4 model genes in lasso analysis.

Correlation analysis of risk model and clinical factors

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).

FIGURE 5

Figure 5. Correlation analysis of risk model and clinical factors (A–D) comparison of risk scores of different clinical subgroups of TCGA glioma patients.

FIGURE 6

Figure 6. Correlation analysis of risk model and clinical factors (A) comparison of differences in risk scores between tumor and normal samples. (B) The ROC curves for risk score differentiation between Tumor and Normal samples. (C) The ROC curves for Grade4 and Grade2 + Grade3 samples with risk score differentiation. (D) The ROC curves for GBM and non-GBM samples with risk score differentiation.

TABLE 3

Table 3. Risk and clinical data in TCGA.

TABLE 4

Table 4. Risk and clinical data in CGGA.

Independent prognostic analysis of risk model

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.

FIGURE 7

Figure 7. Independent prognostic analysis of risk model (A,B) univariate Cox and multivariate Cox independent prognostic analysis (C) the nomogram in TCGA glioma sample to predict 1, 3 and 5-year survival rate (the C-index was 0.8617) (D) the correction curve of nomogram.

Correlation of risk score with KEGG pathway

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.

FIGURE 8

Figure 8. Correlation of risk scores with KEGG pathway.

Impact of risk model on immune heterogeneity

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).

FIGURE 9

Figure 9. Impact of risk model on immune heterogeneity (A) violin diagram of immune scores, stromal scores and ESTIMATE scores for high and low risk groups. The horizontal coordinates represent the groupings (high-/low- risk groups) and the vertical coordinates represent the different scores, blue for the low risk group and pink for the high risk group. (B,C) The box line plots and stacked bar graphs show the abundance of immune cells.

FIGURE 10

Figure 10. Impact of risk model on immune heterogeneity (A) the box line plots comparing the difference of immune cells between high and low risk groups. (B) The box line plots comparing the difference of HLA between high and low risk groups.

Single-cell functional analysis of model genes

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).

FIGURE 11

Figure 11. Single-cell functional analysis of model genes (A–D) functional status of (A) NMNAT2, (B) PARP9, (C) CD38, (D) BST1.

Discussion

As the most common neuroepithelial tumor of the central nervous system, gliomas are a heterogeneous group of central nervous system tumors (17). A unique feature of glioma cells is the high demand for energetic molecules namely adenosine triphosphate (ATP), to sustain their rapid growth rate and the biosynthesis of DNA and proteins (18, 19). Products of NAD+ degradation, NAD+ pathway metabolites and many molecules converted from NAD+ were closely associated with the development and progression of gliomas (12, 14, 20–22). Nevertheless, there were no systematic researches about NAD+ metabolism related genes to predict the survival of patients with glioma. Therefore, it was necessary to establish a new glioma prognostic model based on NAD+ metabolism related genes to verify the prognostic molecular markers and potential therapeutic targets for glioma patients.

Currently, numerous glioma prognostic models have been reported. These glioma prognosis models can be divided into the following types (23): RNA models, methylation models, other multimolecular models. However, due to limitations, such as the lack of external validation, lack of gold markers, not providing quantitative results or absolute risk stratification, difficulties of data collection, the complexity of analysis, and low adherence to complete and transparent reporting, these models are still not recommended for clinical application (23, 24). The NAD+-related model constructed in this study contained only four model genes, it passed internal validation and external validation, was an accurate glioma model.

In this study, based on consensus cluster analysis of 42 NAD related genes (NMRGs), we first divided glioma patients into cluster 1 and cluster 2, with significant survival difference between the two clusters. Then, based on the differential expression analysis of genes between the two groups, Cox analysis and LASSO regression analysis, we further screened out four genes significantly related to the survival of glioma patients, and then established a risk scoring model and nomogram to evaluate their prognostic prediction and clinical application value. Through the correlation analysis between risk score and immune invasion, immunotherapy, KEGG pathway and single cell function enrichment analysis, we further explored the potential functions and signal pathways of the four model genes involved in glioma, with a view to providing more candidate targets for the treatment of glioma patients.

As one of the PARP family, PARP9 was highly expressed in gliomas and high PARP9 expression was associated with poor prognosis and clinicopathological features (15), which supported our findings. In addition, it was shown that PARP9 played an important role in the immune microenvironment of glioma and that PARP9 may be a prepotential immunotherapeutic target for glioma (15). CD38 mediated intracellular ATP levels and glioma cell survival (25). Moreover, CD38 deficiency regulated microglia activation through microglia-associated mechanisms to attenuate glioma progression, and it could modulate tumor-associated microglia/macrophage characteristics (16), which was also consistent with our study. However, to the best of our knowledge, BST1 and NMNAT2 were found to be closely associated with glioma for the first time. BST1 regulated nicotinamide riboside metabolism through glycohydrolase and base exchange activity, which has beneficial effects on aging and aging-related diseases (26). NMNAT2 was a target for the treatment of many malignancies. For example, upregulation of NMNAT2 was associated with the presence, depth and stage of colorectal cancer (27). Furthermore, specific deletion of NMNAT2 in mouse oocytes interfered with meiotic apparatus assembly and metabolic activity, where supplementation with niacin or forced expression of NMNAT2 in aging oocytes during in vitro culture reduced reactive oxygen species (ROS) production and the incidence of spindle/chromosome defects (28). Therefore, it is speculated that NMNAT2 may also act on glioma through NAD+ metabolism.

Some scholars reported that glioma patients with high immune scores and stromal scores were correlated with increased malignancy and reduced survival (29, 30). Our result, which showed that immune scores, stromal scores and ESTIMATE scores were higher in high-risk group than in low-risk group, was consistent with these reports. The majority of the non-neoplastic cells in gliomas are glioma-associated microglia/macrophages (GAMs), comprised macrophages of peripheral origin and brain-intrinsic microglia, which supported tumor progression (31). GAMs are usually divided into two categories, macrophages M1 and macrophages M2. Macrophages M2 are the most abundant immune cells within the glioma stroma (32, 33). Our study showed that Mast cells activated, NK cells activated, and B cells naive were less in the high-risk group, which supported this view. In addition, Tanzhu et al. obtained a similar result (34).

HLA family genes are essential for immunosurveillance and cancer immunotherapy (35, 36). Many HLA molecules, for example, HLA-E, HLA-F, HLA-G, HLA-H and HLA-DR were reported to be overexpressed on cancer cells (35–38). In our model, all 24 HLA family genes were higher in the high-risk group. More data is needed to verify it in the future. In general, most immune checkpoints (including PD-1, TIM-3, CD96, PDCD1, IDO1, PDCD1LG2, and CTLA-4) were highly expressed in glioma cells (39–42). Higher expressions of immune checkpoints were observed in more severe grades of glioma, and this indicated that it was linked to a worse prognosis (39–42). Our result showed that all the immune checkpoints were high in the high-risk group except for TIGIT, suggesting that they synergistically regulated the immune response in the tumor microenvironment. These findings opened up new possibilities for the treatment of gliomas.

An enquiry into the functional status of biomarkers in glioma showed that PARP9 was positively associated with EMT, cell metastasis and cell differentiation and negatively associated with DNA damage. It has been shown that PARP9 was essential for B cell development (43). In addition, the PARP inhibitor veliparib could act on DNA damage repair in prostate cancer cells (44), which supported our study. CD38 showed a significant negative correlation with both cell invasion and apoptosis. Previous studies suggested that CD38 was involved in cell differentiation and inflammatory processes and played a key role in the inflammatory process of autoimmunity (45, 46). NMNAT2 was significantly and positively associated with stemness and proliferation. Wu et al. found that NAD+ deficiency due to reduced NMNAT2 expression affected cell metabolism and meiosis (28), which provided evidence for our study.

There are several limitations in our study. First, the amount of data used in the analysis is not large, so our results may have certain deviation. More data is needed to validate this model in the future. Second, it is the result of bioinformatics analysis without experimental verification. More basic experiments are required to verify the specific mechanism of these genes in glioma. Third, more prospective studies are needed to prove the prognostic function of the four genes. In the future, we will also collect more clinical samples for in-depth in vitro and in vivo experimental research to further verify the relevant mechanisms of the four genes in glioma.

Conclusion

For the first time we built a NAD+ metabolism related risk model, The model could predict the prognosis of glioma patients from another new perspective. Most importantly, it proved that CD38, NMNAT2, PARP9, and BST1 might be important prognostic molecular markers and potential therapeutic targets for glioma patients.

Methods

Data source

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).

TABLE 5

Table 5. Data source.

Classification of glioma subtypes

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).

Screening of NMR-DEGs

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).

Building and validating of a risk model

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.

$\mathrm{risk}\mathrm{score}=\sum _{n=1}^n\left(\mathrm{coe}{\mathrm{f}}_i\times x_i\right)$

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.

Correlation analysis of risk model and clinical factors

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.

Independent prognostic analysis of risk model

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.

Correlation of risk scores with KEGG pathway

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.

Impact of risk model on immune heterogeneity

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.

Single-cell functional analysis of model genes

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.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author contributions

JW: designed this study and contributed substantially to the design of the search strategy. XC: searched and selected the trials and extracted data. MW and WW: performed the analysis and interpreted the data. YW and JX: wrote the manuscript. BZ and XL: critically reviewed the manuscript. HY, HZ and ML: participated in the data extraction and critically revised it. WX and XB: proofread the final version. All authors contributed to the article and approved the submitted version.

Funding

This work was supported in part by the National Natural Science Foundation of China (grant no. 82173285); the Shaanxi International Cooperation Project (grant no. S2017-ZDYF-GJHZ-MSXM-0115).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsurg.2023.1071259/full#supplementary-material.

Data sheet 1
Highly expressed DEGs of glioma in TCGA database.

Data sheet 2
Low expressed DEGs of glioma in TCGA database.

Data sheet 3
Correlation of risk Score with KEGG pathway.

References

1. Velazquez VJ, Brat DJ, Ryken TC, Olson JJ. The role of neuropathology in the management of newly diagnosed glioblastoma: a systematic review and evidence-based clinical practice guideline. J Neurooncol. (2020) 150:143–64. doi: 10.1007/s11060-020-03616-3

PubMed Abstract | CrossRef Full Text | Google Scholar

2. McKinnon C, Nandhabalan M, Murray SA, Plaha P. Glioblastoma: clinical presentation, diagnosis, and management. Br Med J. (2021) 374:n1560. doi: 10.1136/bmj.n1560

CrossRef Full Text | Google Scholar

3. Chini CCS, Zeidler JD, Kashyap S, Warner G, Chini EN. Evolving concepts in NAD+ metabolism. Cell Metab. (2021) 33:1076–87. doi: 10.1016/j.cmet.2021.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Covarrubias AJ, Perrone R, Grozio A, Verdin E. NAD(+) metabolism and its roles in cellular processes during ageing. Nat Rev Mol Cell Biol. (2021) 22:119–41. doi: 10.1038/s41580-020-00313-x

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lee HC, Deng QW, Zhao YJ. The calcium signaling enzyme CD38 - a paradigm for membrane topology defining distinct protein functions. Cell Calcium. (2022) 101:102514. doi: 10.1016/j.ceca.2021.102514

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lundt S, Ding S. NAD(+) metabolism and diseases with motor dysfunction. Genes. (2021) 12:1776. doi: 10.3390/genes12111776

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Eckert MA, Coscia F, Chryplewicz A, Chang JW, Hernandez KM, Pan S, et al. Proteomics reveals NNMT as a master metabolic regulator of cancer-associated fibroblasts. Nature. (2019) 569:723–8. doi: 10.1038/s41586-019-1173-8

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Reiten OK, Wilvang MA, Mitchell SJ, Hu Z, Fang EF. Preclinical and clinical evidence of NAD(+) precursors in health, disease, and ageing. Mech Ageing Dev. (2021) 199:111567. doi: 10.1016/j.mad.2021.111567

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Rodriguez KM, Cohen MS. Chemical genetic methodologies for identifying protein substrates of PARPs. Trends Biochem Sci. (2022) 47:390–402. doi: 10.1016/j.tibs.2021.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Bhatt V, Tiwari AK. Sirtuins, a key regulator of ageing and age-related neurodegenerative diseases. Int J Neurosci. (2022):1–26. doi: 10.1080/00207454.2022.2057849. [Online ahead of print]35549800

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zapata-Perez R, Wanders R, van Karnebeek C, Houtkooper RH. NAD(+) homeostasis in human health and disease. Embo Mol Med. (2021) 13:e13943. doi: 10.15252/emmm.202113943

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Garcia-Peterson LM, Li X. Trending topics of SIRT1 in tumorigenicity. Biochim Biophys Acta Gen Subj. (2021) 1865:129952. doi: 10.1016/j.bbagen.2021.129952

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Lucena-Cacace A, Umeda M, Navas LE, Carnero A. NAMPT As a dedifferentiation-inducer gene: NAD+ as core axis for glioma cancer stem-like cells maintenance. Front Oncol. (2019) 9:292. doi: 10.3389/fonc.2019.00292

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Gasparrini M, Audrito V. NAMPT: a critical driver and therapeutic target for cancer. Int J Biochem Cell Biol. (2022) 145:106189. doi: 10.1016/j.biocel.2022.106189

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu H, Chai S, Wang Y, Wang J, Xiao D, Li J, et al. Molecular and clinical characterization of PARP9 in gliomas: a potential immunotherapeutic target. CNS Neurosci Ther. (2020) 26:804–14. doi: 10.1111/cns.13380

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Levy A, Blacher E, Vaknine H, Lund FE, Stein R, Mayo L. CD38 Deficiency in the tumor microenvironment attenuates glioma progression and modulates features of tumor-associated microglia/macrophages. Neuro Oncol. (2012) 14:1037–49. doi: 10.1093/neuonc/nos121

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bush NA, Chang SM, Berger MS. Current and future strategies for treatment of glioma. Neurosurg Rev. (2017) 40:1–14. doi: 10.1007/s10143-016-0709-8

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Pavlova NN, Zhu J, Thompson CB. The hallmarks of cancer metabolism: still emerging. Cell Metab. (2022) 34:355–77. doi: 10.1016/j.cmet.2022.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lemasters JJ. Metabolic implications of non-electrogenic ATP/ADP exchange in cancer cells: a mechanistic basis for the Warburg effect. Biochim Biophys Acta Bioenerg. (2021) 1862:148410. doi: 10.1016/j.bbabio.2021.148410

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Parsons RB, Facey PD. Nicotinamide N-methyltransferase: an emerging protagonist in cancer macro(r)evolution. Biomolecules. (2021) 11:1418. doi: 10.3390/biom11101418

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Makhov P, Uzzo RG, Tulin AV, Kolenko VM. Histone-dependent PARP-1 inhibitors: a novel therapeutic modality for the treatment of prostate and renal cancers. Urol Oncol. (2021) 39:312–15. doi: 10.1016/j.urolonc.2020.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Colombo G, Gelardi E, Balestrero FC, Moro M, Travelli C, Genazzani AA. Insight into nicotinamide adenine dinucleotide homeostasis as a targetable metabolic pathway in colorectal cancer. Front Pharmacol. (2021) 12:758320. doi: 10.3389/fphar.2021.758320

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Liang X, Wang Z, Dai Z, Zhang H, Cheng Q, Liu Z. Promoting prognostic model application: a review based on gliomas. J Oncol. (2021) 2021:7840007. doi: 10.1155/2021/7840007

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Molinaro AM, Wrensch MR, Jenkins RB, Eckel-Passow JE. Statistical considerations on prognostic models for glioma. Neuro-Oncology. (2016) 18:609–23. doi: 10.1093/neuonc/nov255

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ma Y, Jiang J, Ying W. CD38 Mediates the intracellular ATP levels and cell survival of C6 glioma cells. Neuroreport. (2014) 25:569–73. doi: 10.1097/WNR.0000000000000139

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yaku K, Palikhe S, Izumi H, Yoshida T, Hikosaka K, Hayat F, et al. BST1 Regulates nicotinamide riboside metabolism via its glycohydrolase and base-exchange activities. Nat Commun. (2021) 12:6767. doi: 10.1038/s41467-021-27080-3

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Qi J, Cui C, Deng Q, Wang L, Chen R, Zhai D, et al. Downregulated SIRT6 and upregulated NMNAT2 are associated with the presence, depth and stage of colorectal cancer. Oncol Lett. (2018) 16:5829–37. doi: 10.3892/ol.2018.9400

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Wu X, Hu F, Zeng J, Han L, Qiu D, Wang H, et al. NMNAT2-mediated NAD(+) generation is essential for quality control of aged oocytes. Aging Cell. (2019) 18:e12955. doi: 10.1111/acel.12955

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang C, Cheng W, Ren X, Wang Z, Liu X, Li G, et al. Tumor purity as an underlying key factor in glioma. Clin Cancer Res. (2017) 23:6279–91. doi: 10.1158/1078-0432.CCR-16-2598

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Chaudhary R, Morris RJ, Steinson E. The multifactorial roles of microglia and macrophages in the maintenance and progression of glioblastoma. J Neuroimmunol. (2021) 357:577633. doi: 10.1016/j.jneuroim.2021.577633

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zeiner PS, Preusse C, Golebiewska A, Zinke J, Iriondo A, Muller A, et al. Distribution and prognostic impact of microglia/macrophage subpopulations in gliomas. Brain Pathol. (2019) 29:513–29. doi: 10.1111/bpa.12690

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Andersen JK, Miletic H, Hossain JA. Tumor-associated macrophages in gliomas-basic insights and treatment opportunities. Cancers. (2022) 14:1319. doi: 10.3390/cancers14051319

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Tanzhu G, Li N, Li Z, Zhou R, Shen L. Molecular subtypes and prognostic signature of pyroptosis-related lncRNAs in glioma patients. Front Oncol. (2022) 12:779168. doi: 10.3389/fonc.2022.779168

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang Y, Fan X, Li H, Lin Z, Bao H, Li S, et al. Tumor border sharpness correlates with HLA-G expression in low-grade gliomas. J Neuroimmunol. (2015) 282:1–6. doi: 10.1016/j.jneuroim.2015.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Benevolo M, Mottolese M, Tremante E, Rollo F, Diodoro MG, Ercolani C, et al. High expression of HLA-E in colorectal carcinoma is associated with a favorable prognosis. J Transl Med. (2011) 9:184. doi: 10.1186/1479-5876-9-184

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Roychowdhury A, Samadder S, Das P, Mazumder DI, Chatterjee A, Addya S, et al. Deregulation of H19 is associated with cervical carcinoma. Genomics. (2020) 112:961–70. doi: 10.1016/j.ygeno.2019.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Hazini A, Fisher K, Seymour L. Deregulation of HLA-I in cancer and its central importance for immunotherapy. J Immunother Cancer. (2021) 9:e2899. doi: 10.1136/jitc-2021-002899

CrossRef Full Text | Google Scholar

39. Ghouzlani A, Kandoussi S, Tall M, Reddy KP, Rafii S, Badou A. Immune checkpoint inhibitors in human glioma microenvironment. Front Immunol. (2021) 12:679425. doi: 10.3389/fimmu.2021.679425

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Qi Y, Liu B, Sun Q, Xiong X, Chen Q. Immune checkpoint targeted therapy in glioma: status and hopes. Front Immunol. (2020) 11:578877. doi: 10.3389/fimmu.2020.578877

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Romani M, Pistillo MP, Carosio R, Morabito A, Banelli B. Immune checkpoints and innovative therapies in glioblastoma. Front Oncol. (2018) 8:464. doi: 10.3389/fonc.2018.00464

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Bacolod MD, Barany F, Fisher PB. Can CpG methylation serve as surrogate markers for immune infiltration in cancer? Adv Cancer Res. (2019) 143:351–84. doi: 10.1016/bs.acr.2019.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Robert I, Gaudot L, Yelamos J, Noll A, Wong HK, Dantzer F, et al. Robust immunoglobulin class switch recombination and end joining in Parp9-deficient mice. Eur J Immunol. (2017) 47:665–76. doi: 10.1002/eji.201646757

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yin L, Liu Y, Peng Y, Peng Y, Yu X, Gao Y, et al. PARP Inhibitor veliparib and HDAC inhibitor SAHA synergistically co-target the UHRF1/BRCA1 DNA damage repair complex in prostate cancer cells. J Exp Clin Cancer Res. (2018) 37:153. doi: 10.1186/s13046-018-0810-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Piedra-Quintero ZL, Wilson Z, Nava P, Guerau-de-Arellano M. CD38: an immunomodulatory molecule in inflammation and autoimmunity. Front Immunol. (2020) 11:597959. doi: 10.3389/fimmu.2020.597959

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Shubinsky G, Schlesinger M. The CD38 lymphocyte differentiation marker: new insight into its ectoenzymatic activity and its role as a signal transducer. Immunity. (1997) 7:315–24. doi: 10.1016/S1074-7613(00)80353-2

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Li C, Zhu Y, Chen W, Li M, Yang M, Shen Z, et al. Circulating NAD+ metabolism-derived genes unveils prognostic and peripheral immune infiltration in amyotrophic lateral sclerosis. Front Cell Dev Biol. (2022) 10:831273. doi: 10.3389/fcell.2022.831273

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Briere G, Darbo E, Thebault P, Uricaru R. Consensus clustering applied to multi-omics disease subtyping. BMC Bioinformatics. (2021) 22:361. doi: 10.1186/s12859-021-04279-1

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wilkerson MD, Hayes DN. Consensusclusterplus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–73. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Chen H, Boutros PC. Venndiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. (2011) 12:35. doi: 10.1186/1471-2105-12-35

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Shurson GC, Ku PK, Waxler GL, Yokoyama MT, Miller ER. Physiological relationships between microbiological status and dietary copper levels in the pig. J Anim Sci. (1990) 68:1061–71. doi: 10.2527/1990.6841061x

PubMed Abstract | CrossRef Full Text | Google Scholar