This study integrated multiple independent datasets to construct and validate a prognostic model for low-grade glioma. Initially, we established a histone modification gene set based on Füllgrabe et al. (2011)’s research and the GeneCards database (https://www.genecards.org/) (relevance score >20). Subsequently, single-cell transcriptome data from 4 LGG samples (GSE182109) were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) for further single-cell analysis. TCGA-LGG (513 samples) and GTEx (2,642 normal control samples) datasets were merged after batch effect removal to calculate differentially expressed genes in LGG. Additionally, TCGA-LGG served as the training set, while two independent cohorts from the Chinese Glioma Genome Atlas (CGGA) database (https://www.cgga.org.cn/), containing 693 and 325 GBM samples respectively with clinical and survival information, were used as external validation sets.

WGCNA analysis was conducted based on the previously obtained histone-related differentially expressed gene set (Xu et al., 2023) to reveal gene co-expression relationships and their associations with phenotypes. During data preprocessing, gene expression data from the TCGA-LGG dataset were normalized, and quality was ensured by removing genes with zero standard deviation and samples with missing values. Subsequently, the goodSamplesGenes function was used for quality assessment, and outlier samples were detected using hierarchical clustering (Grabski et al., 2023). In network construction, the optimal soft threshold power = 5 was determined by analyzing the scale-free topology fit index and average connectivity under different soft thresholds. The blockwiseModules function was used to construct the co-expression network, setting the minimum module size to 50 genes, module merging similarity threshold to 0.15, and using unsigned network type. Through calculating Module Eigengenes (MEs) (Han et al., 2019), we analyzed the correlation between modules and histone scores, using Pearson correlation coefficients to evaluate module-trait relationships, and assessed correlation significance using Student’s t-test. The correlation intensity between modules and phenotypic features was visualized through heatmaps. Finally, in modules significantly correlated with histone scores, the biological significance was validated by analyzing the relationship between Module Membership and Gene Significance, with scatter plots demonstrating their positive correlation, further confirming these modules’ central role in the histone regulatory network.

Differential expression analysis was performed using the TCGA-LGG dataset and GTEx dataset through DESeq2, with selection criteria of absolute logFC >1 and p-value <0.05. Feature genes were selected by intersecting the differentially expressed genes with module genes obtained from WGCNA.

To construct a reliable prognostic prediction model, this study used the TCGA-LGG dataset as the training set and CGGA325 and CGGA693 datasets as independent validation sets. Initially, all datasets underwent standardization and feature space consistency was ensured, with missing values eliminated through strict data quality control. During model construction, we systematically evaluated multiple machine learning algorithms and their combinations, including Random Survival Forest (RSF), Elastic Net (Enet), Stepwise Cox regression (StepCox), CoxBoost, Partial Least Squares Regression (plsRcox), Principal Component Regression (SuperPC), Gradient Boosting Machine (GBM), Survival Support Vector Machine (survival-SVM), Ridge regression, and Lasso regression as base models. To enhance prediction performance, we explored various combinations of these base models, such as RSF with CoxBoost and Lasso with GBM combinations. For the Elastic Net model, performance was optimized by adjusting the α parameter (0.1–0.9); for stepwise Cox regression, forward, backward, and bidirectional feature selection methods were employed. Model evaluation used C-index (Song et al., 2022) as the primary evaluation metric, assessing predictive ability and generalization performance through comprehensive performance on training and two independent validation sets. Finally, model performance across different datasets was visualized through heatmaps, and models were ranked based on average C-index values from validation sets to select the final prognostic prediction tool with optimal predictive efficacy and stability.

Based on the comprehensive evaluation results of machine learning models, the optimal prediction model underwent thorough performance validation and risk score system construction. Specifically, through feature importance analysis, 20 key features contributing most significantly to prognosis prediction were identified from the original features. Based on the model’s predictions, a risk score system was established, dividing patients into high and low-risk groups using the median score as the threshold, and Kaplan-Meier survival analysis was performed to observe survival differences between risk groups. Additionally, to evaluate the model’s time-dependent predictive ability, time-dependent ROC curves (Obuchowski and Bullen, 2018) and corresponding AUC values were calculated for 1-year, 3-year, and 5-year predictions.

The prognostic value of clinical features and risk scores was evaluated through systematic survival analysis (Schober and Vetter, 2018). Initially, univariate Cox regression analysis assessed the impact of clinical features including age, gender, tumor type, grade, and risk score on patient prognosis. Subsequently, multivariate Cox regression analysis identified independent prognostic factors. Based on significant independent prognostic factors, an integrated nomogram prediction model (Gittleman et al., 2019) was constructed, and its prediction accuracy was verified through calibration curves for 1-year, 3-year, and 5-year survival predictions. Time-dependent C-index analysis was used to compare the predictive performance between the nomogram model and individual predictive factors, while decision curve analysis (DCA) evaluated the clinical net benefit of the model at different decision thresholds to validate this prediction tool’s value in clinical practice.

This study explored molecular pathway differences between high and low risk score groups through systematic functional enrichment analysis (Canzler and Hackermüller, 2020). First, Gene Set Enrichment Analysis (GSEA) evaluated significantly enriched Hallmark pathways in the high-risk group. Subsequently, Gene Set Variation Analysis (GSVA) scored all samples, and limma differential analysis identified pathways with significantly different activities between high and low-risk groups. To verify the clinical relevance of key pathways, patients were divided into high and low pathway activity groups based on GSVA scores, with Kaplan-Meier survival analysis evaluating prognostic differences, and Cox proportional hazards regression model calculating hazard ratios (HR) and their 95% confidence intervals. For significantly correlated pathways, forest plots were generated to visualize their prognostic value. Finally, survival curves were verified for the six most significant pathways, comprehensively assessing these pathways’ potential roles in glioma development and progression.

This study explored the association between tumor heterogeneity and risk scores through comprehensive analysis of somatic mutation data from the TCGA-LGG cohort. First, the Mutant-Allele Tumor Heterogeneity (MATH) score (Timar and Kashofer, 2020) was used to quantify tumor heterogeneity levels for each sample, comparing differences between high and low-risk groups. Distribution characteristics were visualized through violin plots, with statistical significance assessed using the Wilcoxon rank-sum test. Subsequently, patients were divided into high and low heterogeneity groups based on the median MATH score, with Kaplan-Meier survival analysis evaluating the relationship between tumor heterogeneity and prognosis. Further combining MATH scores with risk scores, patients were classified into four subgroups (high MATH/high risk, high MATH/low risk, low MATH/high risk, low MATH/low risk) to explore the joint predictive value of both indicators. Finally, maftools was used to analyze mutation characteristics of high and low-risk groups, with waterfall plots displaying distribution characteristics of top 20 mutated genes, while co-occurrence and mutual exclusivity analysis (Zang et al., 2023) revealed interaction patterns among key driver genes.

This study conducted systematic analysis of the LGG tumor immune microenvironment using multiple algorithms. Initially, the ESTIMATE algorithm (Tennant et al., 2022) calculated stromal scores, immune scores, and overall scores for each sample, comparing differences between high and low-risk groups. Subsequently, immune-related pathways were scored using ssGSEA, with heatmaps visualizing differential patterns of immune pathway activity between risk groups. Further, the CIBERSORT algorithm (Chen et al., 2018) was employed to deconvolute the infiltration proportions of 22 immune cell types, with violin plots showing immune cell composition differences between high and low-risk groups. Additionally, ssGSEA analysis was performed using 28 immune cell characteristic gene sets, with box plots clearly displaying abundance differences of various immune cell types between risk groups. Finally, correlation analysis explored relationships between key gene expression and immune cell infiltration, as well as associations between risk scores and immune cell infiltration levels. Significant correlation patterns were displayed through heatmaps and correlation scatter plots, revealing potential connections between the risk score model and tumor immune microenvironment.

Systematic drug sensitivity prediction analysis was performed on the TCGA-LGG cohort using the pRRophetic package (Yan et al., 2022). Initially, half maximal inhibitory concentration (IC50) values (Sebaugh, 2011) were predicted for all available drugs based on drug response data from the Cancer Genome Project (CGP) database. For each drug, drug sensitivity differences between high and low-risk groups were compared, with statistical significance assessed using the Wilcoxon rank-sum test. For drugs showing significant differences (P < 0.05), box plots were used to visually display drug sensitivity distribution characteristics across different risk groups.

Based on previous analysis results, we selected the top five key genes with the highest weights for expression validation in the HPA database (https://www.proteinatlas.org/). Through immunohistochemical staining images, protein expression levels and distribution patterns of these key genes were visually demonstrated in normal brain tissue and glioma tissues of different grades. The validation results from the HPA database not only confirmed the differential expression characteristics of these genes in glioma development at the protein level but also provided histological evidence for understanding their potential roles in tumor progression.

To comprehensively visualize the pan-cancer analysis results, we constructed a stratified forest plot integrating survival analysis outcomes across 29 cancer types, organized by eight major organ systems (Breast and Gynecologic System, Respiratory System, Digestive System, Urinary System, Nervous System, Endocrine System, Male Reproductive System, Hematologic System, and Others). In the forest plot, dot sizes represent -log10 (p-value), error bars indicate 95% confidence intervals, and color-coding displays risk levels (red for high-risk groups, cyan for low-risk groups, and gray for non-significant differences). A horizontal dashed line (HR = 1) serves as a reference, and hazard ratios are presented on a logarithmic scale to better illustrate the relative magnitude of risk differences.

Through dimensionality reduction clustering analysis and cell type annotation of LGG single-cell data (Figure 2A), we successfully identified five major cell subpopulations: SOX2 and OLIG2-expressing LGG tumor cells, CD68, AIF1, and FOLR2-expressing macrophages, MBP-expressing oligodendrocytes, PDGFRB-expressing pericytes, and CCL5, CD3E, CD8A, and CD44-expressing T cells. To validate the accuracy of cell type annotation, we constructed a dot plot displaying the expression patterns of marker genes for each cell type (Figure 2B).

To deeply explore the heterogeneity of histone modification levels among different cell types, we calculated and visualized histone modification scores for each cell subgroup (Figure 2C). Results showed significant differences in histone modification levels among different cell subgroups. Further statistical analysis (Figure 2D) revealed that LGG tumor cells and oligodendrocytes exhibited higher histone modification levels (scores >0.25). Based on the overall cellular histone modification levels (HMs), we separated LGG cells into high HMs and low HMs groups (Figure 2E). Differential expression analysis identified 5,278 differentially expressed genes (|logFC|>1, p < 0.05), with 3,638 genes upregulated in the high HMs group and 1,640 genes upregulated in the low HMs group.

Differential expression analysis revealed significant transcriptomic differences between high and low HMs groups. Volcano plot analysis showed numerous genes with significant differential expression (Figure 3A), suggesting these genes may participate in the histone modification regulatory network. To deeply analyze key regulatory genes, we visualized the top 50 up- and downregulated genes with the most significant differences in a circular plot (Figure 3B). To systematically identify co-expression modules related to histone modifications, we performed WGCNA analysis on differentially expressed genes. The hierarchical clustering dendrogram displayed gene co-expression relationships, while the bottom heatmap reflected HMs score variation patterns among samples (Figure 3C). Based on the dynamic tree-cutting algorithm, we ultimately identified 15 functional modules with significant co-expression characteristics (Figure 3D). Module-trait correlation analysis indicated that the magenta module (225 genes) showed the strongest positive correlation with histone modification scores (cor = 0.44, P < 0.005) (Figure 3E). Further module membership analysis revealed that genes in the magenta module showed significant positive correlation between Gene Significance (GS) and Module Membership (MM) (cor = 0.58, P < 1e-21) (Figure 3F), strongly suggesting this module plays a core role in the histone modification regulatory network.

Through integrative transcriptomic analysis of TCGA-LGG and GTEx datasets, we initially identified 6,672 LGG-related differentially expressed genes (DEGs), including 5,798 upregulated and 874 downregulated genes. Intersection analysis of these DEGs with previously determined magenta module genes yielded 129 LGG-specific histone modification-related genes (LGG-HMRgenes) (Figure 4A). Functional enrichment analysis revealed these LGG-HMRgenes were significantly enriched in pathways including heterocycle catabolic process, nucleobase-containing compound catabolic process, RNA catabolic process, and negative regulation of cellular macromolecule biosynthetic process (Figure 4B, P < 0.05). These results emphasize the regulatory role of histone modifications in LGG development and provide theoretical basis for developing targeted therapeutic strategies.

To construct a robust prognostic prediction model, we systematically evaluated the predictive performance of 105 machine learning algorithms using 129 HMR genes as feature inputs (Figure 4C). Through comprehensive comparison of C-index performance across validation sets, the Random Survival Forest (RSF) model demonstrated optimal predictive performance. Based on feature importance analysis of the RSF model, we further selected 20 core feature genes with the strongest predictive contributions (Figure 4D).

To systematically evaluate the time-dependent predictive performance of the HMRS model, we first conducted time-dependent receiver operating characteristic (ROC) analysis in the TCGA-LGG cohort (513 samples). Results showed that the model demonstrated excellent discriminative ability in 1-year, 3-year, and 5-year survival predictions, with corresponding areas under the curve (AUC) reaching 0.77, 0.73, and 0.71 respectively (Figure 5A). After stratifying patients into high and low-risk groups based on the optimal cutoff value, Kaplan-Meier survival analysis revealed significant prognostic differences between the groups (log-rank test, P < 0.001) (Figure 5B).

To validate the external applicability of the HMRS model, we conducted validation in two independent CGGA validation cohorts (CGGA325 and CGGA693). In the CGGA325 cohort (325 samples), the model demonstrated time-dependent prediction accuracy comparable to the training set (Figure 5C), with survival stratification differences showing statistical significance (P < 0.001) (Figure 5D). These results were further confirmed in the CGGA693 cohort (693 samples) (Figures 5E, F). Multi-center validation results confirmed that the HMRS model possesses robust prognostic prediction capability and broad clinical applicability.

To systematically evaluate the prognostic value of clinical features and HMRS scores, we first conducted Cox proportional hazards regression analysis. Univariate analysis showed that age, WHO grade, and HMRS score were significant prognostic factors (Figure 6A, all P < 0.001). Multivariate analysis further confirmed the independent prognostic value of these three factors (Figure 6B). Based on these independent prognostic factors, we constructed an integrated nomogram prediction model. Calibration curve analysis showed that the model demonstrated excellent calibration in 1-year, 3-year, and 5-year survival predictions (Figure 6C).

Clinical decision curve analysis (DCA) indicated that the integrated nomogram model demonstrated greater net benefit compared to single prognostic factors (Figure 6D). Dynamic analysis of time-dependent C-index showed that the nomogram model’s prediction accuracy (C-index >0.70) consistently outperformed single prognostic factors (Figure 6E). Notably, although age performed well in short-term (1–3 years) prediction, its long-term prediction stability was insufficient. Finally, we established a visual nomogram incorporating HMRS score, WHO grade, and age (Figure 6F), providing an intuitive quantitative tool for clinical prognostic assessment.

To elucidate the molecular biological basis of HMRS prognostic stratification, we conducted systematic functional enrichment analysis between high and low-risk groups. GSEA showed that the high-risk group was significantly enriched in multiple cancer-related Hallmark pathways, including Allograft Rejection, E2F Targets, Interferon Gamma Response, MYC Targets V1, and TNFα Signaling via NFκB (Figure 7A, FDR <0.05).

Gene Set Variation Analysis (GSVA) further revealed risk stratification-specific signaling pathway activity characteristics (Figure 7B). The high-risk group showed significant activation of TGF Beta Signaling, Mitotic Spindle, and IL6 JAK STAT3 Signaling; while the low-risk group was characterized by Pancreas Beta Cells, Oxidative Phosphorylation, and KRAS Signaling DN. Correlation analysis between HMRS scores and these pathway activity scores further validated these findings (Figure 7C, P < 0.05).

To assess the clinical prognostic significance of key pathways, we selected six most significant signaling pathways for survival analysis, including Epithelial Mesenchymal Transition, Angiogenesis, Glycolysis, Apoptosis, Coagulation, and IL2 STAT5 Signaling. Kaplan-Meier analysis showed that high activity in these pathways was significantly associated with poorer overall survival (Figure 7D, all P < 0.001). Hazard ratio (HR) analysis further confirmed these outcomes and found that the role of Hedgehog signaling pathway as the sole protect prognostic factor (Figure 7E).

To deeply understand the genomic characteristics of LGG patients, we conducted systematic analysis of histone modification gene mutation patterns and tumor heterogeneity. Using MATH (Mutant-Allele Tumor Heterogeneity) scores to quantify intratumoral heterogeneity levels, results showed significant difference between high score group and low score group (Figure 8A, p < 0.001). Survival analysis based on MATH scores indicated that low MATH scores were significantly associated with poorer prognosis (Figure 8B, p = 0.005). Further analysis integrating MATH scores with HMRS risk stratification showed that the “low risk + low MATH” subgroup had the most favorable prognosis (Figure 8C, p < 0.001).

Through systematic analysis of mutation landscapes in high and low-risk groups (Figures 8D, E), we found: 1) TP53, as a key tumor suppressor gene, had a mutation frequency of 63% in the high-risk group, significantly higher than 28% in the low-risk group; 2) The transcription repressor CIC had a mutation frequency of 4% in the high-risk group compared to 36% in the low-risk group; 3) IDH1 mutations, associated with specific cytogenetic abnormalities and 1p/19q codeletion, showed mutation frequencies of 68% and 85% in high and low-risk groups respectively; 4) ATRX gene, involved in transcriptional regulation and chromatin remodeling, had mutation frequencies of 46% and 20% in high and low-risk groups respectively; 5) Additionally, characteristic mutations in the high-risk group included TTN (16%), while the low-risk group included FUBP1 (15%).

Mutation co-occurrence analysis (Figures 8F, G) revealed that in the high-risk group, TP53 mutations significantly co-occurred with IDH1 and ATRX. In the low-risk group, besides observing co-occurrence patterns of TP53, IDH1, and ATRX, significant mutation co-occurrence characteristics were also found between COL6A3 and PTEN.

We conducted multi-dimensional analysis of the immune microenvironment in high and low-risk groups. ESTIMATE algorithm assessment results showed (Figures 9A–C) that the high-risk group had significantly higher stromal scores, immune scores, and overall scores than the low-risk group (p < 0.001), suggesting more active immune responses and stromal components in the high-risk group.

Furthermore, ssGSEA algorithm analysis revealed all 15 significantly different immune-related pathways between high and low-risk groups (Figure 9D).

CIBERSORT algorithm analysis of immune cell infiltration characteristics showed (Figures 9E, F): 1) Memory CD4⁺ T cells were significantly higher in the high-risk group (P < 0.001); 2) Memory B cells, CD8⁺ T cells, and follicular helper T cells were more abundant in the low-risk group (P = 0.001); 3) Except for CD56dim NK cells, other immune cells generally showed higher expression levels in the high-risk group.

Correlation analysis between key gene expression and immune cell infiltration (Figure 9G) revealed that histone modification genes might participate in LGG progression by regulating immune cell infiltration. Correlation analysis between risk scores and immune cell infiltration (Figure 9H) indicated: 1) Significant positive correlations with memory CD4⁺ T cells, regulatory T cells, dendritic cells, and neutrophils; 2) Significant negative correlations with naive T cells and memory B cells. These findings suggest that HMRS can effectively quantify the immune status of LGG patients, reflecting significant immune landscape differences among patients with different risk levels.

Based on the risk score model, we predicted sensitivity differences to common drugs between high and low-risk groups. Through comparison of IC50 values (Figure 10), significant response differences were found for the following drugs: Temozolomide (A), Pazopanib (B), Paclitaxel (C), Rapamycin (D), Sorafenib (E), and Gefitinib (F). These findings provide important references for risk stratification-based individualized medication.

To further validate the expression characteristics of key genes in the risk model, we selected five genes with the highest weight coefficients (ADK, UGCG, RPN2, CAPZA1, KDELR2) for HPA database immunohistochemical validation (Figure 11). Results showed that all these genes exhibited significantly high expression in LGG tissue.

Our pan-cancer analysis revealed distinct survival patterns across different cancer types and organ systems (Figure 12). Among the 29 cancer types analyzed, four cancer types demonstrated statistically significant associations with survival outcomes (p < 0.05). In the Urinary System, both Kidney Renal Clear Cell Carcinoma (KIRC, HR = 1.72, p = 0.00016) and Kidney Renal Papillary Cell Carcinoma (KIRP, HR = 2.87, p = 0.00077) showed significantly higher risk in the high-risk group. Within the Respiratory System, Lung Adenocarcinoma (LUAD) exhibited significantly poorer outcomes (HR = 1.57, p = 0.0014). Additionally, in the others category, Sarcoma (SARC) demonstrated significantly worse survival (HR = 1.59, p = 0.024).

Although not reaching statistical significance, several other cancer types showed notable trends:

1. Adrenocortical Carcinoma (ACC) displayed a strong tendency toward higher risk (HR = 2.06, p = 0.065).

Notably, most cancer types (24/29) maintained hazard ratios above 1.0, suggesting a consistent trend toward worse outcomes in high-risk groups across different cancer types, although not all reached statistical significance. However, some cancer types, including Cholangiocarcinoma (CHOL), Head and Neck Squamous Cell Carcinoma (HNSC), Acute Myeloid Leukemia (LAML), Ovarian Serous Cystadenocarcinoma (OV), and Testicular Germ Cell Tumors (TGCT), showed hazard ratios below 1.0, indicating potentially better outcomes in the high-risk group, though these associations were not statistically significant.

In this study, we first revealed the cellular heterogeneity characteristics of histone modifications in LGG through multi-omics analysis and developed a prognostic prediction model with clinical application prospects. Single-cell level analysis showed that tumor cells and oligodendrocytes exhibited higher levels of histone modifications, which is consistent with previous reports on the key role of histone modifications in glioma stem cell maintenance (Liu et al., 2024; Zhu et al., 2024; Liu et al., 2023). High levels of histone modifications may promote tumor cell proliferation and stemness maintenance through precise regulation of gene expression networks (Sharma et al., 2023; Shen et al., 2020), and this epigenetic level cellular heterogeneity may be one of the important reasons leading to LGG treatment resistance and recurrence. Based on this finding, we constructed the HMRS prognostic model, which demonstrated excellent predictive performance across multiple independent cohorts. The model’s prediction accuracy was further improved after integration with clinical features. The advantages of the HMRS model lie in its robustness, practicality, and individualization characteristics, providing a new tool for LGG patients’ prognostic assessment and treatment decision-making. These findings not only deepen our understanding of LGG epigenetic heterogeneity (Chaligne et al., 2021) but also provide an operational prognostic assessment method for clinical practice.

Compared with currently widely used prognostic assessment models, the HMRS model shows unique advantages. Traditional prognostic assessments mainly rely on WHO grading, IDH mutation status, and 1p/19q codeletion as molecular markers. Although these indicators have important prognostic implications, they often cannot fully reflect the molecular heterogeneity and dynamic evolution characteristics of tumors (Eckel-Passow et al., 2015). Recently developed radiomics-based prediction models, such as those combining MRI imaging features with machine learning (Li et al., 2021b), although advantageous in non-invasive assessment, still need improvement in prediction accuracy and stability. In contrast, our HMRS model not only integrates histone modifications as an important epigenetic feature but also ensures its predictive reliability through multi-center validation. Notably, the HMRS model can reflect the epigenetic heterogeneity of tumor cells, giving it potential advantages in predicting treatment response and guiding individualized treatment. Furthermore, the inclusion of clinical features makes it more easily applicable in actual clinical work, an advantage not possessed by other single molecular markers or complex models.

Through functional enrichment analysis of high and low-risk groups, we deeply revealed the molecular mechanisms behind HMRS stratification. The study found multiple key signaling pathways significantly activated in the high-risk group, most notably the TGF-β and IL6-JAK-STAT3 signaling pathways. TGF-β signaling pathway activation may enhance tumor cell invasion and metastasis capabilities through inducing epithelial-mesenchymal transition (EMT) (Meng et al., 2016), while IL6-JAK-STAT3 pathway activation suggests the important role of inflammatory microenvironment in LGG progression, with this chronic inflammatory state potentially promoting tumor malignant progression through multiple mechanisms (Zhang et al., 2023). Meanwhile, we observed significant differential expression of metabolism-related pathways, reflecting the metabolic heterogeneity developed by tumor cells to adapt to malignant proliferation. These findings not only explain the intrinsic mechanisms of prognostic differences between different risk stratification patients at the molecular level but also provide theoretical basis for developing new therapeutic strategies, particularly targeted therapy against these abnormally activated pathways as potential treatment options for high-risk patients. Importantly, the interactions between these pathways form a complex regulatory network, suggesting the need to consider multi-target combined intervention when developing treatment strategies.

Compared to previous studies, the uniqueness of our research lies in being the first to systematically reveal the role of histone modification-mediated signaling pathway networks in LGG progression. Previous studies mainly focused on single pathways, such as IDH mutation-mediated metabolic reprogramming reported by Wang et al. (2023), or Guo et al. (2023)’s discovery of Notch signaling regulation in glial cell development and tumorigenesis. Through integrative analysis, our study not only confirmed the importance of these known pathways but also discovered complex regulatory relationships between them. Notably, we found that the TGF-β signaling pathway may influence global epigenetic states by regulating histone-modifying enzyme expression, a finding that echoes Mao et al. (2020)’s recent discovery in glioma but is the first report in LGG. Furthermore, our study first revealed potential connections between inflammatory pathways and metabolic reprogramming, providing new perspectives for understanding LGG heterogeneity and theoretical basis for developing multi-target combination therapy strategies.

Further pathway analysis revealed comprehensive activation of both innate and adaptive immune response networks in high-risk patients. The enrichment of pattern recognition receptor pathways (RIG-I-like, NOD-like, and Toll-like receptor signaling) indicates heightened innate immune surveillance, potentially triggered by tumor-derived danger signals. The concurrent activation of B cell receptor signaling pathway, T cell receptor signaling pathway, and natural killer cell-mediated cytotoxicity suggests broad engagement of adaptive immune responses. The upregulation of leukocyte transendothelial migration and chemokine signaling pathways points to active immune cell trafficking within the tumor microenvironment. Additionally, the enrichment of antigen processing and presentation pathways, along with the complement and coagulation cascades, indicates robust immune recognition and response mechanisms. However, despite this extensive immune activation, the apparent ineffectiveness in tumor control suggests potential immune dysfunction or suppression. The enhanced FC gamma R-mediated phagocytosis pathway might reflect increased clearance of antibody-coated tumor cells, yet the overall immune response appears insufficient to prevent disease progression in high-risk patients. These pathway alterations, combined with the observed immune cell composition changes, paint a picture of a complex but potentially dysfunctional immune response that may contribute to tumor progression (Giannone et al., 2020; Kaur et al., 2022).

Then, our study revealed close associations between HMRS risk stratification and tumor immune microenvironment through multi-dimensional analysis. Through ESTIMATE algorithm analysis, we found that high-risk group patients showed significantly elevated immune scores and stromal scores, suggesting the existence of a more complex immune regulatory network. Notably, we observed significantly increased memory CD4⁺ T cells in the high-risk group, while other immune cells such as memory B cells, CD8⁺ T cells, and follicular helper T cells showed relative deficiency. This immune cell component remodeling may reflect the establishment of tumor immune escape mechanisms. Meanwhile, the activation of multiple immune-related signaling pathways, including RIG-I-like receptor signaling pathway, B cell receptor signaling pathway, and Toll-like receptor signaling pathway, further supports the key role of immune microenvironment in LGG progression. These findings not only deepen our understanding of the LGG immune microenvironment but more importantly provide new perspectives for developing immunotherapy strategies. Particularly for high-risk patients, rebuilding effective anti-tumor immune responses, such as enhancing memory CD8⁺ T cell function or regulating specific immune pathway activity, may become important strategies for improving treatment efficacy. Additionally, these immune characteristic differences suggest the need to consider individualized immune microenvironment differences when designing treatment plans, potentially requiring treatment strategy adjustments based on patients’ immune status, including whether to combine immune checkpoint inhibitors and other immunotherapy approaches.

In terms of drug applications, temozolomide, the standard first-line treatment for LGG, showed significant therapeutic differences between high-risk and low-risk groups (Tomar et al., 2021), providing direct guidance for clinical medication decisions. Meanwhile, we observed that the high-risk group demonstrated good sensitivity to certain multi-target tyrosine kinase inhibitors such as pazopanib (Miyamoto et al., 2018), which is consistent with the abnormal pathway activation patterns we previously identified. Interestingly, high-risk group patients also showed increased sensitivity to traditional chemotherapy drugs like paclitaxel (Alqahtani and Aleanizy, 2019), suggesting that cell cycle regulation may be a crucial factor affecting drug response. Individualized treatment strategies not only hold promise for improving therapeutic outcomes but may also reduce unnecessary drug toxicity, thereby enhancing patients’ quality of life.

Our pan-cancer analysis demonstrates that the model initially developed for LGG exhibits predictive value across multiple cancer types, achieving statistical significance particularly in kidney cancers (KIRC, KIRP) and lung adenocarcinoma (LUAD). This cross-organ system predictive capacity can be understood through the lens of evolutionary conservation in developmental biology. Although these tissues originate from different germ layers (neuroectoderm vs. mesoderm), they share several key regulatory networks during embryonic development.

Notably, kidney development involves complex mesenchymal-epithelial transition (MET) and neuroectoderm-mesoderm interactions, which share many similarities with cell fate determination mechanisms in neural system development. For instance, WT1 and the PAX gene family play crucial roles in both kidney and neural system development (Discenza et al., 2003; Ochi et al., 2022). Similarly, lung development requires precise EMT and complex cell lineage determination, sharing multiple regulatory pathways with neural crest cell migration and differentiation mechanisms (Jolly et al., 2018). These common developmental features may explain why epigenetic markers derived from neural system tumors demonstrate significant predictive value in these cancer types.

From an evolutionary developmental perspective, this conservation of predictive features reflects the shared pluripotent state and fundamental regulatory mechanisms of different tissues during early embryonic development. Despite subsequent divergence in differentiation pathways, the basic epigenetic regulatory networks remain highly conserved throughout evolution.

Although our study has made several important discoveries in LGG histone modifications and prognostic prediction, there are still some limitations that need to be addressed in future research. The primary limitation is the lack of prospective clinical validation; although our model has demonstrated good predictive performance across multiple independent cohorts, its performance in real clinical settings still needs to be verified through prospective studies. Secondly, the functional mechanisms of key genes still need more experimental data support, particularly in vivo and in vitro functional experiments will help to deeply understand these genes’ specific roles in LGG progression. Additionally, while drug sensitivity analysis provides important clues for individualized treatment, these predicted results still need to be validated through standardized clinical trials for accuracy and reliability. Finally, the dynamic nature of epigenetic modifications suggests that longitudinal sampling might provide additional insights not captured in our current cross-sectional analysis. Based on these limitations, we have planned several important future research directions: first, we will conduct multi-center prospective clinical studies to systematically evaluate the clinical application value of the HMRS model; second, through in-depth molecular biology experiments, we will elucidate the regulatory mechanisms of key genes and their roles in tumor progression; third, we will explore individualized combination treatment strategies based on risk stratification, particularly in optimizing combinations of immunotherapy and targeted therapy; finally, we plan to integrate multi-modal data such as radiomics and metabolomics to develop more comprehensive and precise prediction models. These studies will help to further improve the diagnosis and treatment level of LGG patients, ultimately achieving better clinical outcomes.