Integrative Analysis of a Novel Eleven-Small Nucleolar RNA Prognostic Signature in Patients With Lower Grade Glioma

Objective The present study used the RNA sequencing (RNA-seq) dataset to identify prognostic snoRNAs and construct a prognostic signature of The Cancer Genome Atla (TCGA) lower grade glioma (LGG) cohort, and comprehensive analysis of this signature. Methods RNA-seq dataset of 488 patients from TCGA LGG cohort were included in this study. Comprehensive analysis including function enrichment, gene set enrichment analysis (GSEA), immune infiltration, cancer immune microenvironment, and connectivity map (CMap) were used to evaluate the snoRNAs prognostic signature. Results We identified 21 LGG prognostic snoRNAs and constructed a novel eleven-snoRNA prognostic signature for LGG patients. Survival analysis suggests that this signature is an independent prognostic risk factor for LGG, and the prognosis of LGG patients with a high-risk phenotype is poor (adjusted P = 0.003, adjusted hazard ratio = 2.076, 95% confidence interval = 1.290–3.340). GSEA and functional enrichment analysis suggest that this signature may be involved in the following biological processes and signaling pathways: such as cell cycle, Wnt, mitogen-activated protein kinase, janus kinase/signal transducer and activator of tran-ions, T cell receptor, nuclear factor-kappa B signaling pathway. CMap analysis screened out ten targeted therapy drugs for this signature: 15-delta prostaglandin J2, MG-262, vorinostat, 5155877, puromycin, anisomycin, withaferin A, ciclopirox, chloropyrazine and megestrol. We also found that high- and low-risk score phenotypes of LGG patients have significant differences in immune infiltration and cancer immune microenvironment. Conclusions The present study identified a novel eleven-snoRNA prognostic signature of LGG and performed a integrative analysis of its molecular mechanisms and relationship with tumor immunity.


INTRODUCTION
Tumors derived from the neuroepithelium are collectively called gliomas and are the most common malignant primary intracranial tumor (1). Lower-grade gliomas are welldifferentiated gliomas, and their preferred treatment is surgery, followed by radiotherapy, chemotherapy, immunotherapy, and targeted therapy. Just like any other cancers, glioma is also due to the interaction of innate genetic risk factors and environmental carcinogenic factors (e.g., radiation exposure). Some known genetic disease, such as neurofibromatosis (type I), is genetic susceptibility factor for glioma (2). With the development of high-throughput sequencing technology, we have found that more and more genetic variants are connected with the occurrence, development and clinical outcome of cancers. TCGA is a database resource for whole-genome multi-omics sequencing of 33 human cancers based on high-throughput sequencing technology (3). It can provide the whole genome multi-omics data set of glioma, so that we can further explore its genetic variation (4).
SnoRNA belong to a small non-coding RNA, mostly enriched in nucleolus (5). SnoRNA plays an important role in the splicing processing and post-transcriptional modification of ribosomal precursor RNA. Previous research reports suggested that snoRNAs were dysregulated between glioma tumor and adjacent non-tumor tissues, and involved in the regulation of cancer cell apoptosis, proliferation and invasion, including in glioma (6)(7)(8). However, there are few reports on the prognostic value of snoRNA in gliomas (8). We found that there is no research report on screening LGG prognostic snoRNAs based on genome-wide dataset. To fill this gap, we designed and implemented this study. The present study was used the RNA sequencing (RNA-seq) dataset to screen prognostic-related snoRNAs and construct a prognostic signature for The Cancer Genome Atla (TCGA) lower grade glioma (LGG) cohort, and conduct a comprehensive investigation of its molecular mechanisms and relationship with tumor immunity.

Data Source
The original level 3 raw RNA sequencing (RNA-seq) dataset of lower grade glioma used in our study comes from TCGA data portal (4). Clinical parameters of TCGA LGG cohort were obtained from UCSC Xena (http://xena.ucsc.edu) (9). We obtained 940 snoRNAs from TCGA RNA-seq dataset into the current analysis. The original HTSeq-Counts RNA-seq data set was normalized in the R platform through the edgeR package (10). We eliminated snoRNAs with a mean value less than 1, and obtained a total of 137 snoRNAs for subsequent survival analysis. We obtained a total of 529 RNA-seq data sets of 511 patients from the TCGA data portal, and we excluded 18 recurrent cancer tissues. Of the 511 patients, we excluded five patients whose survival time was 0 or who had no survival time. We also excluded 18 patients who had a history of other malignancies or who received radiotherapy or chemotherapy before surgery. Finally, 488 patients were recruited for the subsequent comprehensive analysis, which included survival analysis and functional mechanism analysis. IDH1 mutation data of TCGA LGG cohort were obtained from cBioPortal for Cancer Genomics (https://www.cbioportal.org/index.do).

Survival Analysis of snoRNAs Signature
Prognostic snoRNAs screening was performed on the R platform using survival packet and univariate Cox proportional hazards regression model. The prognostic snoRNAs signature is constructed in the R platform using the step function. We use prognostic-related snoRNAs as variables into the multivariate Cox proportional hazard regression model, and the Cox regression coefficient (b) as the weight of each snoRNAs included in the signature to compose the risk score. The snoRNAs signature calculation formula was as follows: risk score = Exp of snoRNA 1 × b 1 snoRNA 1 + Exp of snoRNA 2 × b 2 snoRNA 2 + Exp of snoRNA n × b n snoRNA n (Exp: expression) (11)(12)(13). Subsequently, we also used the time-dependent ROC curve, nomogram, and combined effect survival analysis to conduct a comprehensive analysis of the prognostic value of this risk score signature in LGG. The nomogram was drawn in the R platform using the rms package. The time-dependent ROC curve was performed in the R platform using the survival ROC package. The combined effect survival analysis was performed in SPSS version 22.0 using the multivariate Cox proportional hazards regression model.

Functional Enrichment Analysis of snoRNAs Signature Based on Whole-Gene RNA-seq Dataset
SnoRNA is a small non-coding RNA encoded by introns, and its function is exerted through the regulation of protein-coding genes (PCGs). We use the Cor function and whole genome encoded RNA-seq dataset of LGG tumor tissues samples to screen the PCGs that related to snoRNAs, were performed on the R platform. We identified that the absolute value of Pearson correlation coefficient (r) between snoRNAs and genes were greater than 0.4 and P <0.05 considered to be the coexpression PCGs of snoRNA. The normalization of PCGs' RNA-seq dataset was also performed by edgeR. Then, we use the Database for Annotation, Visualization and Integrated Discovery v6.8 (DAVID v6.8, https://david.ncifcrf.gov/home. jsp) online tool to perform functional annotations on the snoRNAs we have screened through the snoRNA coexpression PCGs (14). In addition, we also used the gene set enrichment analysis (GSEA, http://software.broadinstitute.org/ gsea/index.jsp) desktop installer to perform functional enrichment analysis between high-risk and low-risk phenotypes for purpose of further investigate the molecular mechanisms of prognostic differences between this two phenotypes (15). For further investigate the molecular mechanism, we further used edgeR to screen differentially expressed genes (DEGs) between high-risk and low-risk phenotypes. Then we use DAVID v6.8 to annotate the functions of these DEGs. We subsequently used these DEGs and connectivity map (CMap, https://portals.broadinstitute.org/ cmap/) online tool to identify potential drugs for LGG patients with a high-risk phenotype (16,17). PubChem (https://pubchem. ncbi.nlm.nih.gov) and STITCH (http://stitch.embl.de/cgi/) online tools were applied to obtain the chemical structure of drugs and the drug-gene interaction networks, respectively (18,19).

Tumor Immune Infiltration and Microenvironment Analysis
CIBERSORT package was used to compute the abundance and scale of 22 immune cells in LGG tumor tissues to assess the degree of immune infiltration of LGG tumor samples (20,21). We use the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) package and LGG's whole genome RNA-seq dataset to calculate the tumor immune microenvironment on the R platform, including the values of immune and stromal cells (22).

Statistical Analysis
The calculation of FDR adopts Benjamini-Hochberg method (23). Survival analysis uses log-rank test of Kaplan-Meier algorithm, as well as univariate and multivariate Cox proportional hazards regression models. Clinical parameters with P <0.05 were included in multivariate Cox proportional hazards regression model for adjustment. Survival curves, heat maps and volcano maps were drawn on the R platform using the ggplot2 package. P <0.05 identified the difference to reach statistical significance. All statistics are computed by SPSS version 22.0.

Survival Analysis of snoRNAs Signature
Clinical parameters of TCGA LGG cohort were summarized in Table 1. Through survival analysis, we obtained 21 snoRNAs that are markedly connected with the clinical outcome of LGG ( Figure 1 and Table S1). Then we use the step function to screen a risk score signature from these 21 prognostic-related snoRNAs. Through the screening of the step function, we finally constructed a signature composed of 11 prognostic snoRNAs. The expression of risk score is below: risk score = (-0. LGG patients and risk score, we found that patients with high risk score had poorer clinical outcome than those with low risk score ( Figures 3A, B, log rank P <0.0001). We used multivariate Cox proportional hazard regression model to analyze risk score and found that LGG patients with high risk score had a higher risk of death (adjusted P = 0.003, adjusted HR = 2.076, 95%CI = 1.290-3.340). Then we used survivalROC to assess the accuracy of this risk score for predicting the clinical outcome of LGG patients. We found that this risk score had the highest accuracy in assess the one year prognosis of LGG patients, and the area under curve (AUC) was 0.850 ( Figure 3C). The accuracy of the prognosis signature exceeded 0.7 in the clinical outcome prediction of LGG patients from one to five years ( Figure 3C). We also use this prognostic signature and clinical parameters with P <0.01 to perform a joint effect survival analysis, these clinical parameters including age, grade, IDH1 mutation, radiation therapy and tumor location. The joint effect survival analysis can significantly more accurately classify LGG patients and identify subtypes of LGG patients with different prognostic outcomes ( Table 2 and Figures 4A-E). Through the combined effect survival analysis, we found that the combination of age and risk score can significantly divide LGG patients into four subtypes with significant different prognosis ( Table 2 and Figure 4A, all adjusted P <0.05). For understand the contribution of this prognostic signature to the clinical outcome of LGG patients, we constructed a nomogram base on eight LGG prognostic-related clinical parameters and risk score. Through the nomogram, we found that the risk score has the highest contribution to LGG survival, with scores ranging from 0 to 100 ( Figure 5).

Functional Enrichment Analysis of snoRNAs Signature
Through genome-wide co-expression analysis, 990 co-expression pairs of snoRNAs-PCGs were identified, and the snoRNAs-PCGs co-expression interaction network is shown in Figure  S1, since we did not screen the co-expressed PCGs of SNORD91A in this cohort. Finally, we will screen the coexpressed genes of 10 snoRNAs into the subsequent functional enrichment analysis. Enrichment analysis reveals that these snoRNA co-expressed PCGs may work by participating in the following biological functions: DNA repair, neuromuscular process controlling balance, covalent chromatin modification, neurotransmitter secretion, voltage-gated calcium channel activity, NIK/NF-kappaB signaling, nervous system development, neuromuscular junction, neuronal action potential, protein polyubiquitination, chemical synaptic transmission, brain morphogenesis, regulation of long-term neuronal synaptic plasticity, positive regulation of GTPase activity, MAPK cascade, positive regulation of cyclindependent protein serine/threonine kinase activity involved in G1/S transition of mitotic cell cycle, protein kinase binding, regulation of cell cycle, synaptic vesicle cycle and glutamatergic synapse (Table S3). We also performed a multivariate survival analysis on these snoRNA co-expressed PCGs, and obtained 180 prognostic PCGs of LGG ( Figure 6A). The top three significant PCGs are zinc finger protein 821 (ZNF821, Figure 6B), solute carrier family 8 member A3 (SLC8A3, Figure 6C) and cellular communication network factor 4 (CCN4, also named as WISP1, Figure 6D). In order to reveal the molecular mechanisms that cause the prognosis difference of LGG patients with high-and low-risk score phenotypes, we also compared patients with two phenotypes using the GSEA approach. GSEA using the c5 reference gene set suggested that LGG patients with high-risk score phenotypes are significantly different from low-risk score phenotypes in the following biological processes: integrin mediated signaling pathway, natural killer cell mediated immunity, cell adhesion mediated by integrin, kappa B kinase NF kappa B signaling, toll like receptor signaling pathway, T cell mediated immunity, regulation of tumor necrosis factor mediated signaling pathway, T cell receptor signaling pathway, Wnt activated receptor activity, hippo signaling, Wnt protein binding, receptor signaling pathway via STAT, regulation of cell cycle G1/S phase transition, NIK/NF kappa B signaling, tumor necrosis factor mediated signaling pathway, stem cell proliferation ( Figures 7A-P). GSEA using the c2 reference gene set suggested that LGG patients with high-risk score phenotypes are significantly different from low-risk score phenotypes in the following pathways: CD8/TCR pathway, natural killer cell mediated cytotoxicity, B cell receptor signaling pathway, PI3KCI pathway, caspase pathway, TH1/ TH2 pathway, P53 pathway, metastasis up, tumor vasculature up, JAK/STAT signaling pathway, cell adhesion molecules CAMs, P38/MAPK pathway, Wnt pathway require Myc, signaling by Rho gtpases, tumorigenesis up and NF-kB pathway ( Figures S2A-P).
We screened 1319 DEGs between high-and low-risk score phenotypes by edgeR, of them, 122 DEGs were down-regulation and 1197 DEGs were up-regulation ( Figure 8, Figure S3 and Table S7). Functional enrichment revealed that these DEGs may function by participating in the following biological processes: positive regulation of cell proliferation, cell-cell signaling, positive regulation of ERK1 and ERK2 cascade, interferon-gammamediated signaling pathway, regulation of immune response, neuron migration, angiogenesis, positive regulation of cell migration, cytokine-mediated signaling pathway, cell adhesion, cellular response to interleukin-1, cell differentiation, positive regulation of T cell proliferation, positive regulation of cytosolic calcium ion concentration involved in phospholipase C-activating G-protein coupled signaling pathway, cell surface receptor signaling pathway, neuropeptide signaling pathway, JAK-STAT cascade, mitotic sister chromatid segregation, cell division, mitotic nuclear division, integrin binding, cell-matrix adhesion, T cell receptor signaling pathway, response to hypoxia, positive regulation of protein kinase B signaling, epithelial to mesenchymal transition, regulation of vascular endothelial growth factor production, positive regulation of MAPK cascade, cellular response to fibroblast growth factor stimulus, phosphatidylinositol-4,5-bisphosphate 3-kinase activity, Wnt-protein binding, cytokine-cytokine receptor interaction, ECM-receptor interaction, cell adhesion molecules (CAMs), focal adhesion, PI3K-Akt signaling pathway, Proteoglycans in cancer, cell cycle, chemokine signaling pathway, Jak-STAT signaling pathway, and transcriptional misregulation in  (Table S8). We also performed a multivariate survival analysis on these DEGs, and obtained 320 prognostic DEGs of LGG (Table S9 and Figure 9A). The top three significant DEGs were calcium binding protein 4 (CABP4, Figure 9B), elastin microfibril interfacer 3 (EMILIN3, Figure 9C) and ISL LIM homeobox 2 (ISL2, Figure 9D). We use these DEGs through the CMAP to identify the potential targeted therapy drugs of this snoRNAs prognostic signature in LGG. We have obtained ten targeted therapy drugs for this snoRNAs prognostic signature in LGG. These drugs are 15-delta prostaglandin J2, MG-262, vorinostat, 5155877, puromycin, anisomycin, withaferin A, ciclopirox, chloropyrazine and megestrol ( Table 3). The chemical structures of nine drugs (except 15-delta prostaglandin J2) can be obtained through PubChem, which are shown in Figures 10A-I. Subsequently, we also constructed the drug-genes interaction network by STITCH, In addition to 15-delta prostaglandin J2 and 5155877, we can obtain drug-genes interaction network for the remaining 8 drugs ( Figure S4). We found that these drugs can work by regulating genes that are differentially expressed between high and low-risk score phenotypes. MG-262 may function in LGG by

Tumor Immune Infiltration and Microenvironment Analysis
By using CIBERSORT to assess the degree of infiltration of 22 immune cells in LGG tumor tissues, we finally obtained distribution histogram and abundance heatmap of immune cell infiltration of 116 LGG patients ( Figures S5A-B), of which 52 patients were the low-risk score phenotype and 64 patients were high-risk score phenotype. By comparing the distribution of 22 immune cells between high-and low-risk score phenotypes of LGG patients, we found that eight immune cells, including T cell CD8, T cells follicular helper, T cells regulatory (Tregs), NK cells activated, Monocytes, Macrophages M1, Macrophages M2 and Eosinophils, showed significant differences between high-and low-risk score phenotypes ( Figure 11).

Among the eight types of immune cells, four types of immune cells [T cell CD8, T cells regulatory (Tregs), Macrophages M1and
Macrophages M2] in the high-risk phenotype have a higher fraction than in the low-risk phenotype, whereas, four types of immune cells (T cells follicular helper, NK cells activated, Monocytes and Eosinophils) in the low-risk phenotype have a higher fraction than these in the high-risk phenotype ( Figure 11). Subsequently, we also scored the immune cells and stromal cells in tumor tissues of 488 LGG patients. We found that all the stromal, immune, and ESTIMATE score for the tumor microenvironment of LGG patients were higher in the high-risk score phenotype ( Figures  12A-C).

DISCUSSION
There has been a small amount of literature that used the TCGA genome-wide data set to screen snoRNA prognostic biomarkers.   (25). Both of the above studies confirmed that U3 can be used as a prognostic marker for corresponding tumors (24,25). Similarly, Zhao et al. also constructed a renal clear cell carcinoma with six prognosisrelated snoRNA signatures, and conducted functional enrichment analysis on these snoRNAs (26). Similar to previous studies, our current study also used TCGA RNA-seq dataset to investigate prognostic snoRNAs, and also found that U3 can be used as a prognostic biomarker for LGG. Therefore, we inferred that U3 could be used as prognostic biomarker in multiple cancers. Lafaille et al. found that the genetic variation of SNORA31 was closely related to encephalitis in the forebrain. After alleles at corresponding sites of SNORA31 were knocked out, in vitro experiments proved that SNORA31 could change the sensitivity of neurons to herpes simplex virus-1 (27). The expression of SNORA31 was markedly down-regulated in CD19+ b-cells of chronic myeloid leukemia (CLL) patients compared with normal B-cells (28). Studies have shown that relapsing-remitting multiple sclerosis (RRMS) patients ncRNA-mRNA network is seriously affected by the disease, resulting in a large number of ncRNAs and mRNAs imbalance. By means of RNA sequencing, Irizar et al. found that SNORA40 was disorder in RRMS and might be a potential therapeutic target (29). We found that the rest of snoRNAs had not been reported before. We found for the first time that that these snoRNAs are related to LGG prognosis through whole-genome RNA sequencing data. Through the snoRNAs co-expression genes, GSEA and differentially expressed genes, we have a further understanding the function of this eleven-snoRNA prognostic signature. A large number of pathways and biological function gene sets enriched by this eleven-snoRNA prognostic signature have been reported as classic or novel cancer-related signaling pathways in previous studies, such as JAK/STAT, p38//MAPK, and Wnt signaling pathways. Xiong et al. screened the differentially expressed genes between ovarian serous cystadencinoma (OSC) and ovarian specimens, and used CMap method to determine that three drugs, namely MG-132, puromycin and 15-delta prostaglandin J2, could be potential target drugs for OSC (30). Shi et al. also identified 15-delta Prostaglandin J2 as a potential therapeutic agent by comparing gene expression profiles in patients with diabetic nephropathy using bioinformatics analysis (31). , and its potential target in HNSCC is proliferating cell nuclear antigen (32). Vorinostat in the treatment of glioma has been widely reported in previous studies. Clinical trials have shown that vorinostat combined with Bevacizumab or Temozolomide can be used in the treatment of glioma (33)(34)(35), as well recurrent malignant gliomas (36)(37)(38)(39).
Local Delivery of vorinostat can alter the tumor microenvironment, thereby inhibiting glioma recurrence (40). Vorinostat can also mediate the regulation of cell cycle regulatory proteins in glioma (41). Bortezomib can enhance the apoptosis of glioma cells induced by vorinostat (42). Puromycin exerts antitumor function in multiple cancers and can induce cancer cell apoptosis. The aminopeptidase inhibitors based on puromycin show high anti-tumor effect in vitro in hematologic  malignancies and are expected to be potential therapeutic drugs for hematologic diseases (43). Puromycin can induce apoptosis of breast cancer cells MCF-7. Its apoptosis mechanism is not exerted by insulin-like growth factor-1 (IGF-1), but is affected by the level of IGF-1 receptor (44). Study have shown that puromycin can induce apoptosis of glioma cells, but its apoptotic mechanism is not exerted through the Fas/Fas ligand pathway and Bcl-2 has only a small protective effect on puromycin-induced apoptosis of glioma cells (45). Previous studies have shown that anisomycin has anticancer effects in a variety of cancers, but its molecular mechanisms are different in different cancers. Anisomycin can inhibit angiogenesis in ovarian cancer (OV) through take part in the ceRNA regulatory network (lncRNA−Meg3/miR−421/platelet derived growth factor receptor a axis), thereby inhibiting tumor cell growth (46). LncRNA b-site APP cleaving enzyme 1 antisense strand (BACE1-AS) is also a potential target of anisomycin in OV (47). Anisomycin plays an anti-tumor role in primary hepatocellular carcinoma by mediating natural killer cells, and its main molecular mechanism has been found by genome sequencing to be mediated by immunomodulatory genes (48). Anisomycin can be involved in inhibiting the proliferation of colorectal cancer cells by mediating GATA binding protein 6 (49). Anisomycin can improve the efficacy of BCR-ABL tyrosine kinase inhibitors in CLL by mediating Wnt/b-catenin signaling pathway (50). Anisomycin can also mediate phosphorylated mitogen-activated protein kinases (MAPKs) p38 and Jun Nterminal kinase (JNK) to increase the apoptosis of glucocorticoids-resistant acute lymphoblastic leukemia cell lines (51). Anisomycin can also increase the sensitivity of melanoma anti-tumor drugs (52). Anisomycin exerting an anti-tumor effect in osteosarcoma through induces cell cycle arrest and apoptosis by mediating mitochondrial biogenesis (53). Anisomycin in renal cell carcinoma cells can mediate Bcl-2/c-FLIP(L)/Mcl-1/death receptor 4 (DR4) to promote cell apoptosis and exert anti-tumor effects (54,55). Anisomycin induces apoptosis of glioma cell lines (U251 and U87 cell lines) by down-regulating PP2A catalytic subunit in in vitro cell experiments (56). Studies found that withaferin A has great potential in anti-tumor therapy (57,58), as well as in brain tumor (59,60). The combination of Withaferin A and tumor treating fields can improve its anti-tumor effect (61). The potential molecular mechanism of Withaferin A in glioma may be exerted by NF-KB nuclear translocation, activation of caspase   (64,65). Ciclopirox is also a widely reported anti-tumor drug (66)(67)(68)(69). However, after reviewing the literature, we have not found an anti-tumor study of ciclopirox in glioma. Megestrol exert anti-tumor effects in multiple cancers, but has not been reported in the treatment of glioma (70)(71)(72). Among the ten drugs we screened, 5155877 and chloropyrazine have not been reported in antitumor studies. This study still has some limitations that need to be clarified. First, this study is a single-center study and lacks a validation cohort. Therefore, our results need further validation in a multi-center, large sample cohort. Second, the molecular mechanism analysis and drug prediction of snoRNA in this study are based on bioinformatics analysis methods, so our results require further verification by in vivo and in vitro experiments. Despite the above limitations, our research is still the first study to screen LGG prognosis snoRNAs based on whole-genome RNA-seq dataset. Our research is also the first to conduct a comprehensive analysis of LGG prognosis snoRNAs, including molecular mechanism, targeted drug prediction, tumor immune infiltration and tumor microenvironment. The above results can provide a basis for future study on clinical application and molecular mechanisms of snoRNAs in LGG.

CONCLUSION
In conclusion, our current study have identified 21 prognostic snoRNAs and constructed a novel eleven-snoRNA prognostic signature for LGG patients. GSEA and functional enrichment analysis suggest that this signature may be involved in the following biological processes and signaling pathways: such as cell cycle, Wnt, mitogen-activated protein kinase, janus kinase/ signal transducer and activator of tran-ions, T cell receptor, nuclear factor-kappa B signaling pathway. CMap analysis screened out ten targeted therapy drugs for this signature: 15-delta prostaglandin J2, MG-262, vorinostat, 5155877, puromycin, anisomycin, withaferin A , ciclopirox, chloropyrazine and megestrol. We also found that high-and low-risk score phenotypes of LGG patients have significant differences in immune infiltration and cancer immune microenvironment. The novel findings of our study are the first comprehensive genome-wide investigation of the prognostic related snoRNAs in LGG and the preliminary identification of a novel prognostic snoRNAs signature. We used the bioinformatics analysis to identify the prognostic value and biological mechanisms of this signature. Since this study is based on a single-center study, our results still need to be further verified.

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.

AUTHOR CONTRIBUTIONS
TD: Study design and preparation of manuscript (equal contributors). YG: Study planning and literature search (equal contributors). XL: Data collection and analysis. XW: Literature analysis. XZ: Literature analysis. GZ: Data interpretation. LM: Study planning. All authors contributed to the article and approved the submitted version.