Molecular subtypes based on PANoptosis-related genes and tumor microenvironment infiltration characteristics in lower-grade glioma

The growth of cancer, the effectiveness of treatment, and prognosis are all closely related to PANoptosis (include pyroptosis, apoptosis, and necroptosis). It remains unclear whether PANoptosis genes (PANGs) may contribute to lower-grade glioma (LGG) tumor microenvironment (TME). In this study, we collected 1203 LGG samples from three public databases and reported that PANoptosis involves TME interaction and prognosis. Firstly, we provided a comprehensive review of the pan-cancer landscape of PANGs in terms of expression characteristics, prognostic value, mutational profile, and pathway regulation. Then, we identified two distinct PANclusters, each with its own molecular, clinical, and immunological profile. We then developed a scoring system for LGG patients called PANscore. As well as investigating immune characteristics, tumor mutational characteristics, and drug sensitivity, we examined the differences between groups with high PANscores and those with low PANscores. Based on this PANscore and clinicopathological variables, an instant nomogram for predicting clinical survival in LGG patients was developed. Our thorough examination of PANGs in LGG revealed their probable function in TME, as well as their clinicopathological characteristics and prognosis. These discoveries could deepen our comprehension of PANGs in LGG and provide doctors fresh perspectives on how to forecast prognosis and create more efficient, individualized treatment plans.


Introduction
Brain gliomas develop from glial cells and are the most common primary malignant tumors in the central nervous system. Today, lower-grade glioma (LGG), which refers to grade II and III oligoastrocytomas, oligodendrogliomas, and astrocytomas (Zeng et al. 2018), is used to define better the features of slow infiltration and relatively slow progression.
There is currently no satisfactory standard of care for gliomas, which consists of surgical resection followed by radiotherapy and temozolomide chemotherapy. With the development of tumor genomics and the introduction of molecular pathology diagnosis of glioma, the treatment of glioma has also entered the era of molecularly targeted therapy (Louis et al. 2021). Molecular targeted therapy has achieved good clinical results in hematological and some solid tumors. However, the current treatment effect in glioma is still poor because the pathogenesis of glioma is still unclear. Patients with gliomas are highly heterogeneous, and conventional pathological histological classification is difficult to predict the prognosis of LGG (Ostrom et al. 2018). As genetic alterations are drivers of glioma formation, new genomics-based subtypes are urgently needed.
In healthy tissues, cell death occurs frequently and plays an important role in maintaining tissue function and morphology, and includes active cell death, programmed death, apoptosis, and passive cell death-necrosis (Barth et al. 2020). There are three main genetic pathways to programmed cell death (PCD): pyroptosis, apoptosis, and necroptosis (Yan et al. 2023). They are involved in homeostasis and disease development in various complex ways. The PANoptosome complex regulates PANoptosis, an inflammatory PCD pathway activated by no specific trigger, which highlights the interplay between cytokinesis, apoptosis, and necroptosis, but whose mechanisms cannot be explained by any of these three PCD pathways alone . Although identifying the key initiators, effectors, and actuators of these three PCD pathways varies, a growing body of evidence highlights their extensive interactions . The inhibitory effects of PANoptosis on tumor growth across the cancer spectrum have been previously revealed, providing an additional clinical rationale for patient-targeted therapy studies ). This suggests a broader role for PANoptosis in cancer, where it promotes tumor cell death and thus improves patient prognosis (Karki et al. 2020). Therefore, understanding the underlying mechanisms of PANoptosis could open up further avenues for developing promising new strategies for cancer treatment.
However, the role of PANoptosis in patient prognosis and the tumor microenvironment (TME) of LGG is unknown. Based on clinical data from patients with LGG obtained from the TCGA and GEO database and expression data from 66 PANoptosis genes identified in prior studies, we conducted a bioinformatic analysis of the present study. We showed that molecular clustering and prognostic characteristics based on PANoptosis could be used to predict survival and the tumor microenvironment in LGG patients. This study provides a further understanding of the role of PANoptosis in LGG and enables the development of a clinical basis for more effective treatment strategies for patients with glioma.

Acquiring and processing data
We downloaded RNA-seq data and corresponding clinical information for LGG patients from the TCGA database (https:// portal. gdc. cancer. gov/). Transcriptome profiles of normal brain tissue were obtained from TCGA and the GTEx project (https:// gtexp ortal. org/ home/). We used FPKM data for subsequent data analysis. For batch effects between the TCG-LGG and GTEx datasets, we used the "limma" method to remove batch effects and used it to draw heat maps and box plots. In addition, we downloaded copy number variation (CNV) data and somatic mutation data from the TCGA and UCSC Xena browsers (http:// xena. ucsc. edu/). Then, we downloaded the sequence matrix file for GSE16011 from the Gene Expression Omnibus (GEO). The TCGA cohort contains 505 LGG samples, and the GEO cohort (GSE16011) contains 106 LGG samples. We used the "Combat" algorithm to eliminate the batch effect and combined the GEO dataset with the TCGA-LGG dataset. In addition, we screened LGG samples from two RNA-seq cohorts in the Chinese Glioma Genome Atlas (CGGA) database (http:// www. cgga. org. cn/) for external validation in this study; these included CGGA325 (172 cases) and CGGA693 (420 cases). According to previous studies, 66 PANoptosis genes were combined with scorched, apoptotic, and necrotic genes, removing overlapping genes to create the PANoptosis gene list (Pan et al. 2022).

Pan-cancer analysis
First, in our study, we convert FPKM values to TPM values and then perform subsequent difference analysis based on "limma" package (|log2(FC)| > 1.5, FDR 0.05). We derived the survival picture for PANoptosis from the TCGA analysis of the association between gene expression and patient survival, with the criterion of hazard ratio HR < 1 for protective genes and HR > 1 for risk genes. The TCGA database also contained data on 10,234 single nucleotide variants (SNV) from 33 cancer types. Seven distinct types of mutation were taken into account for this study. Mutation frequencies for pan-cancer were summarized using percentage heat maps. To identify cases of significantly altered CNV amplification or deletion in the patient group, we analyzed CNV data from 11,495 samples obtained from the TCGA database. Copy number alterations were improved by heterozygosity and purity of amplifications and deletions for each gene, of which more than 5% were considered high-frequency CNV (Feng et al. 2022a;Feng et al. 2022b). On the other hand, we annotated the methylation probes for each gene promoter. Using Wilcoxon signed rank tests, we determined whether tumors and normal samples had different levels of methylation, and genes that were significantly hypermethylated or hypomethylated were determined using a p-value cutoff of 0.05. Gene set variation analysis (GSVA) is an analysis method that enriches gene sets from microarray and RNA-seq data under parameterless and unsupervised conditions. Finally, we used "GSVA" package to calculate pathway activity scores (7876 cancer samples from the TCGA and TCPA databases), and the relevant gene sets were obtained from the previous references (Ye et al. 2018).

A consistent clustering analysis of differentially expressed genes (DEGs) and PANGs
To decrease the dimensionality of consensus unsupervised clustering, univariate cox regression analysis was used to exclude out genes having prognostic significance. Only prognostically relevant genes closely associated with PANoptosis were retained as clustering factors. Using the R software's "ConcensusClusterPlus" program, PANoptosis molecular subtypes were determined using consensus clustering. After selecting the ideal number of clusters between k = 2 and 10, the process was repeated 1000 times to guarantee that the findings were reliable and repeatable. According to PANG expression, LGG patients were classified into different PANclusters. Consequently, it was important to import the gene names together with pertinent sample data and expression levels in order to examine the transcriptional patterns of the various PANoptosis subtypes. This is done by using the "limma" package of the princomp function for PCA and the "ggplot2" package of the R software to present the results. In addition, DEGs from different molecular subtypes (p-value < 0.05 and log2 |fold change (FC) | > 2) were classified into different genomic subtypes.

Prognosis and clinicopathological features of molecular subtypes
For the purpose of discovering the clinical significance of different PANclusters, we compared the clinicopathological characteristics of LGG patients with their genetic subgroups. Age, sex, grade, and tumor histological type were among the clinicopathological variables. With the R survivorship and survminer packages, we performed a survival analysis, and we looked at the relationship between PANclusters and overall survival (OS) in LGG patients using Kaplan-Meier plots and log-rank.

Relationship between different molecular isoforms and TME
To confirm the TME's characterization in many molecular PANclusters, we performed gene variance analysis (GSVA) on the set of marker genes for TME. The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) were implemented by the MSigDB database. We used Gene Set Enrichment Analysis (GSEA) based on two gene sets to understand the specific functional profiles of molecular isoforms. Normalized enrichment scores (NES) with absolute values >1 and nominal p < 0.05 were considered statistically significant. Additionally, tumorinfiltrating immune cells (TICs) were calculated for each LGG sample using the immune cell infiltration algorithm. Based on the single sample GSEA algorithm, the level of TICs in each LGG sample was calculated. Of course, to get a more comprehensive picture of the TIME characteristics of the different molecular subtypes, we simultaneously looked at immune landscape differences using multiple algorithms such as CIBERSORT, Abs, EPIC, IPS, MCPcounter, and Xcell. In addition, the differences between the two groups were verified by estimate algorithms for the StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity. Finally, we also cared for immune checkpoint genes, immunoinhibitors, immunostimulators, and MHC gene expression levels between subtype groups to validate the TIME differences.

Identification of different molecular isoforms of DEGs and functional enrichment analysis
DEGs were identified between different molecular subtypes of copper rot using R's "limma" package. The fold change was 1.5, adjusted for p < 0.05, for screening DEGs. We used the "ClusterProfiler", "org.Hs.eg.db", "enrichplot" and "ggplot2" packages in the R software for GO and KEGG enrichment analysis of DEGs.

Construction of a PANoptosis risk scoring system
We developed a PANoptosis risk scoring system to identify PANoptosis patterns in individual tumors. In order to identify DEGs related with survival in LGG patients, DEGs from various PANoptosis molecular subtypes were initially screened. We next did univariate Cox regression analysis. Secondly, patients were classified into different genetic subtypes of PANoptosis, including PANcluster A and PANcluster B, based on the DEGs affecting prognosis by consensus clustering analysis. Third, using the R "caret" package, the cohort's total number of LGG patients was randomly split into a training group and a test group in a 1:1 ratio. Finally, in the training group, we developed a PANoptosis risk grading system. We also tested the system on trial and with both groups. To reduce the possibility of overfitting, we specifically used the "glmnet" R package for LASSO regression analysis. To create the model, the trajectories of the relevant variables were examined and cross-validated. In order to identify potential PANoptosis risk genes, we once more used multivariate Cox analysis. Using the training data, we developed a predictive PANoptosis risk score system. The PANoptosis risk score was calculated as follows: PANscore = Σ (Expi * coef). Expi denotes the expression of PANoptosis essential risk genes, and coef represents the risk factor. Molecular or genetic subtypes were assessed through correlation analysis. A total of LGG patients in the training group were divided into a high-PANscore group and a low-PANscore group based on the median PANscore. Furthermore, Kaplan-Meier survival analysis was used to compare survival rates between two PANscore groups. Likewise, validation and combined groups were divided into two PANscore groups. ROC curves were plotted for each group based on survival analysis.

Evaluation of TME and different PANscore groups
By using box line plots, we looked at the expression of PANoptosis in two PANscore subgroups. We also calculated TICs in the TME of each LGG sample using different immune cell infiltration algorithms. TICs and prognostic risk genes were studied using correlation analysis. We analyzed the correlation between CSC and TMB in two PANscore groups. We used the "maftools" R package to analyze the frequency of mutations in the two PANscore risk groups.

Drug sensitivity analysis
Furthermore, as previously reported, we carried out a connection map (CMap) study to forecast certain small molecule compounds ). The connection map database (CMap, https:// clue. io/) had 1309 drug signatures in total, and the top 150 upregulated and 76 downregulated expression profiles were chosen as input data. The algorithm limit sum (XSum) was used to obtain CMap scores, and the top four small molecule compounds with the lowest CMap scores were chosen for display. The free binding energy for binding between putative macromolecular proteins and small molecule compounds was then determined using AutoDock Vina software using the small molecule compounds with the lowest CMap scores. Finally, we used PyMOL software to show the molecular docking data. Pan-cancer analysis of PANGs. A Graph showing log2 FC and FDR for PANGs in each cancer. Red and blue indicate genes significantly up-and downregulated. B Based on COX regression, PANG expression and patient OS are correlated in each cancer. Red represents higher expression associated with poorer survival (HR > 1) and blue represents associated with better survival (HR < 1); p < 0.05 is considered statistically significant. C For a given malignancy, PANG mutation frequency represents the number of samples with a mutant gene. D SNV oncoplot. E Bar graph showing the frequency of CNV change for each PANG in each cancer type. F Heat map depicts the differential methylation of PANGs in malignancies; genes that are hypermethylated and hypomethylated are indicated, respectively, in red and blue (Wilcoxon rank-sum test). G Effect of PANGs on cancer-related pathways (FDR≤0.05); the number in each cell indicates the percentage ◂

Development of nomogram scoring system
A nomogram is created using the "rms" package, combining clinical and pathological features and PANoptosis scores based on the patient's survival. By summing the scores of the variables, we were able to obtain the total score. In order to verify the projected values between the anticipated 1-, 3-, and 5-year survival rates and the actual results, follow-up calibration plots of the nomogram scoring system were carried out.

Statistical analyses
A statistically significant difference was defined as p < 0.05, and R version 4.2.2 was used for all statistical calculations. Figure 1 depicts the flowchart of this study. Our next study analyzed 66 PANoptosis such as TNF, TLR3, RIPK3, RIPK1, and MLKL (Table S1). First, we compared tumor samples with normal samples from the TCGA database for expression of PANGs. We found that 42 PANGs were generally upregulated in tumor tissues, including RIPK3, RIPK1, MLKL, and FASLG ( Fig. 2A). Following this, we constructed a survival picture of these genes based on the association between gene expression levels recorded in TCGA and patient survival, summarizing the correlation between PANoptosis expression and patient survival. With LGG. H Two PANclusters were defined using consensus clustering analysis. I CDF curves with consensus scores between k = 2-10 and the area under the curve HR < 1 for protective genes and HR > 1 for risk genes at p < 0.05, we found the most protective genes in KIRC and SARC and the riskiest genes in LGG and UVM (Fig. 2B). Meanwhile, SNV, CNV, and methylation alterations were evident for PANoptosis-related genes, including LGG. TP63, CASP8, UASA, and APAF1 showed significant SNV alterations in most tumor types (Fig. 2C, D). rIPK1, MLKL, FASLG, and FAS showed significant CNV amplification; however, TNF, TLR3, RIPK3, and CHMP4A showed significant CNV deletion in most tumor types (Fig. 2E). PANoptosis genes were differentially methylated relative to normal tissue in most cancers, notably, TP73 and RIPK1 as representative of isogenic mean level hypermethylation, and CASP8 and TNF as representative of isogenic average hypomethylation (Fig. 2F). We then presented detailed information on the potential effects of PANoptosis genes on various biological activity pathways. In brief, it can be seen that PANoptosis genes may generally have a potential activating effect on pathway activity such as apoptosis, which is unquestionable, and conversely, a potentially inhibiting effect on DNA damage pathway activity, and indeed, PANoptosis genes may be associated with many biologically relevant pathways; more information is shown in Fig. 2G.

Identification of PANoptosis molecular subtypes
We first learned the detailed location of CNV alterations in the PANoptosis gene on the chromosome (Fig. 3A). Next, we performed a general analysis of its somatic mutation frequency, which showed a low frequency of mutations in the PANoptosis gene in TCGA-LGG samples, with the highest mutation frequency of only 1% in CASP8 (Fig. 3B). Further CNV frequency analysis showed a prevalence of  Comparisons of immune scores, stromal scores, estimate scores, and tumor purity between PANclusters A and B (M). * indicated p < 0.05, ** indicated p < 0.01, *** indicated p < 0.001 CNV alterations in the PANoptosis gene (Fig. 3C). Our next objective was to compare the expression levels of PANoptosis genes between LGG tumors and standard brain samples. We found that 34 PANoptosis expressions were upregulated in tumor tissue (Fig. 3D, E). There was a significant association between 46 PANGs and overall survival (OS) among LGG patients in a survival analysis (Fig. 3F), and we visualized an integrated landscape of complex relationships between PANoptosis genes and LGG predictive value using a network diagram that suggested a general and complex interaction between PANoptosis genes (Fig. 3G). Regarding the integrated association between PANoptosis genes, we classified LGG patients into two groups based on the expression profile of PANoptosis genes and used the ConsensusClusterPlus algorithm for classification. The findings imply that k=2 might be the most appropriate option for classifying patients into two categories, including PANcluster A (n=307) and B (n=304) (Fig. 3H, I). Additionally, in accordance with the various LGG molecular subtypes, we examined the association between PANoptosis expression and clinical traits such as grade, gender, age, and tissue type (Fig. 4A). Significant variations in PANoptosis-related transcriptional patterns between subtypes A and B were validated by PCA (Fig. 4B). Survival analysis revealed that patients with LGG of subtype A had a greater likelihood of surviving than those with subtype B (Fig. 4C).

TME differences between molecular subtypes
To further characterize the different subtypes of TME, GO GSVA showed that the PANcluster B was significantly enriched in pathways associated with the regulatory function of the immune response, the favorable regulatory mechanism of γ-interferon production, and the signaling pathway mediated by tumor necrosis factor (Fig. 4D). GSEA showed that genes of the PANcluster B were also enriched in the pathways mentioned above in the gene ontology set C5 collection (Fig. 4F). KEGG GSVA also revealed that the PANcluster B was mainly enriched in pathways associated with antigen processing and presentation and interaction between cytokines and cytokine receptors (Fig. 4E). GSEA showed that in the KEGG gene set C2 collection, genes of the PANcluster B were also mainly enriched in the pathways cited earlier (Fig. 4G). The relationship between immune cell infiltration and immune modulators was examined between the two PANclusters in order to determine whether the PANGs may be involved in TIME. Human immune cell subpopulations were calculated for each LGG sample based on seven immune cell infiltration algorithms, including CIBERSORT, CIBERSORT-Abs, EPIC, IPS, MCPcounter, ssGSEA, and Xcell (Fig. 5H). Between PANclusters A and B, there were noticeable changes in the infiltration of the majority of immune cells (Fig. 5A-G). These results suggest that the PANcluster A correlates with the cold immune phenotype and the PANcluster B correlates with the immune hot phenotype. We also observed significant differences in the expression levels of immune checkpoints (Fig. 5I), immunostimulatory molecules (Fig. 5J), immunosuppressive molecules (Fig. 5K), and MHC molecules (Fig. 5L) between PANclusters A and B. These expression levels were generally higher in PANcluster B than in A, confirming subtype B as the immune hot phenotype. Subsequently, based on the ESTIMATE algorithm, subtype B had higher StromalScore, ImmuneScore, and ESTIMATEScore, while subtype A had higher TumorPurity (Fig. 5M).

Identification of PANoptosis gene isoforms based on DEGs
We discovered 3251 DEGs between PANclusters A and B to further examine the probable biological activity of the various molecular isoforms of PANoptosis. Following these analyses, we identified relevant biological pathways based on KEGG and GO enrichment analysis. The results showed that KEGG was enriched in Fc gamma R-mediated phagocytosis, cell adhesion molecules, and endocytosis (Fig. 6A, B). BP regulates chemical synaptic transmission, trans-synaptic signaling, and neuronal death; CC is involved in the synaptic membrane, cell-substrate junction, and neuron-toneuron synapse; MF is involved in tubulin binding, GTP binding, and ubiquitin protein ligase binding (Fig. 6C, D). To assess the predictive value of 2353 DEGs associated with subtypes (p < 0.05), we performed univariate Cox regression analysis (Supplementary File 1). We classified patients into two genomic subtypes, gene clusters A and B, based on a consensus cluster analysis of 2353 DEGs associated with prognosis (Fig. 6E). Patients with genotype B had a greater chance of survival than those with genotype A (Fig. 6F). We also used heat maps to demonstrate the clinical characteristics between genotypes with differences in gene expression used for consensus clustering (Fig. 6G). The two PANoptosis genotypes differed significantly in the face of PANGs; this is consistent with the results of the PANoptosis model (Fig. 6H).

Development and validation of the PANscore
Our study began with dividing LGG patients at random into two groups, one for training (n=306) and one for testing (n=305). We used LASSO and multivariate Cox analysis of 2353 prognostic DEGs to seek the best predictive features (Fig. 7A, B). Next, we conducted a multivariate Cox regression analysis, and five genes were identified (PTGFRN, SULF2, LRRC4, CCNY, ARL3), of which four genes (SULF2, LRRC4, CCNY, ARL3) were negatively associated, and PTGFRN was positively associated. We constructed the PANscoring system in the training set: PANscore = (−0.034386217*LRRC4 expression) + (−0.056261861*CCNY expression) + (0.040297237*PTG-FRN expression) + (−0.057615884*ARL3 expression) + (−0.016362642* SULF2 expression) (Fig. 7C). The distribution of molecular subtypes, genetic subtypes, and risk subtypes in the LGG samples is shown in Fig. 7D. In addition, we calculated risk scores for molecular and genetic subtypes for each LGG sample. We found that molecular subtype B and genetic subtype A had significantly higher risk scores than molecular subtype A and genetic subtype B (Fig 7E, F). The expression of PANs differed significantly between high-and low-PANscore groups (Fig. 7G). Next, based on PANscore, Kaplan-Meier survival curves showed that LGG patients with low PANscores had better survival ( Fig. 8A-E). Scatter plots of PANscore in training and individual validation groups showed that patients' survival time decreased as PANscores increased and that the expression of the five essential risk score genes differed more between the high-and low-PANscore groups (Fig. 8F-J). In addition, the area under the time condensed curve (AUC) values for 1-, 3-, and 5-year survival for PANscore were 0.873, 0.876, and 0.799 in the training set, respectively; good diagnostic efficacy was present in each validation set ( Fig. 8K-O).

Relationship between TME and different groups of PANscores
Next, by applying seven immune cell infiltration algorithms, including CIBERSORT, CIBERSORT-Abs, EPIC, IPS, MCPcounter, ssGSEA, and Xcell, TICs and PANscore were correlated using correlation analysis, and to examining the relationship between five key. (Fig. 9A). A significant association was found between these five genes and the majority of immune cells. Of these, PANscore showed a strong positive correlation with type 1 T helper cell, plasmacytoid, dendritic cell, natural killer cell, natural killer T cell, and activated dendritic cell (Fig. 9B). Therefore, it is possible that the PANscore is related to LGG's TIME. Subsequently, based on the ESTIMATE algorithm, the high PANscore subtype had higher StromalScore, ImmuneScore, and ESTIMATEScore, while the low PANscore subtype had higher TumorPurity (Fig. 9C). Based on their self-renewal capacity and differentiation potential, tumor stem cells have been considered promising targets for cancer therapy. Therefore, we investigated the correlation between PANscore and tumor stem cell index values and showed a linear negative correlation result between PANscores and tumor stem cell index, suggesting that LGG cells with low PANscores have more diverse stem cell properties and a lower degree of cell differentiation (Fig. 9D).
Our analysis of the TCGA-LGG cohort mutation data showed that the high-PANscore group had a higher TMB than the low-PANscore group. A subsequent Spearman correlation analysis found a positive correlation between TMB and PANscore (Fig. 9E). We also discussed variations in the distribution of somatic mutations in the TCGA-LGG cohort between the two sets of PANscore groups using the maftools package. The top three mutated genes in the two PANscore sets were IDH1, TP53, and ATRX, respectively. The frequency of these mutations was significantly higher in patients in the low-PANscore set compared to the high-PANscore set (Fig. 9F).

Analyses of chemotherapy drug sensitivity and screening of small molecules
Chemotherapy remains an indispensable treatment for LGG, and one of the main reasons for poor prognosis in LGG patients is chemoresistance. In this study, we used the "oncoPredict" R package to investigate the potential sensitivity of clinical agents in two PANscore groups. We screened 10 chemotherapeutic agents for the treatment of LGG, such as vorinostat, osimertinib, linsitinib, and lapatinib, among others (Fig. 10A). These drugs showed high sensitivity in high-grade patients, suggesting that high-PANscore patients may benefit from these chemotherapeutic agents. In addition, we used the eXtreme Sum (XSum) approach to perform CMap analysis to screen potential small molecule drugs that we could use to treat LGG. The volcano plot shows DEGs in both PANscore groups, including 178 upregulated and 76 downregulated genes  (Fig. 10B). We screened five potential small molecule drugs based on these DEGs: butein, imatinib, MK-886, arachidonyltrifluoromethane, and TTNPB (Fig. 10C). Finally, to assess the affinity of the drug candidates for their targets, we performed a molecular docking analysis. The binding pose and interaction of the candidate with the highest CMap score, namely butein, with its core molecular target TNF, were obtained using AutoDock Vina v.1.2.2. The results showed that the butein drug bound to its protein target through visible hydrogen bonding and strong electrostatic interactions with low binding energy of −7.91 kcal/mol, indicating a highly stable binding (Fig. 10D). Overall, these results suggest that genetic markers of PANoptosis are associated with drug sensitivity.

Establishment of a nomogram for predicting survival in LGG patients
A nomogram survivor was developed and evaluated based on the PANscore in univariate and multivariate Cox regression studies. For LGG patients, the PANscore was an independent prognostic factor (Fig. 11A, B). Afterwards, multivariate Cox was used to construct a nomogram model for estimating OS at 1, 3, and 5 years. Model variables included age, tumor tissue type (astrocytoma, oligoastrocytoma, and oligodendroglioma), tumor grade (GII, GIII), and PANscore (Fig. 11C). For the analysis of the model's AUC, the nomogram exhibited a high diagnostic value, according to ROC curves (Fig. 11E). Calibration curves demonstrated that the nomogram accurately predicted survival in LGG patients for 1, 3, and 5 years (Fig. 11F). Based on the DCA results, the nomogram model yielded substantial net gains under a wide range of risks (Fig. 11G, H, I). All these findings reveal the strong performance of the nomogram model.

Single-cell analysis and immunohistochemistry of five risk model genes
Using the TISCH glioma database, GSE141460, we explored gene expression in TME by using five risk model genes. This dataset had 15 cell groups annotated and 6 intermediate cell types (Fig. 12A). LRRC4, CCNY, PTGFRN, ARL3, and SULF2 were all expressed in malignant tumor cells (Fig. 12B, C). We obtained IHC staining pictures of five risk model gene-associated proteins in LGG and healthy brain tissue using the HPA database. We utilized them to examine if the levels of protein expression for these five genes varied in LGG. LRRC4, CCNY, PTGFRN, ARL3, and SULF2 protein expression levels were significantly higher and gene expression panels. K-O PANscore's predictive efficiency for 1, 3, and 5 years is illustrated by the ROC curves in LGG samples than in standard samples, supporting these findings (Fig. 13A-D).

Discussion
The mode of cell death usually occurs in a mixed form, as cells can undergo extensive crosstalk under pathological conditions (Zheng and Kanneganti 2020;Karki et al. 2021). PANoptosis has been identified as a complex form of regulated cell death that includes apoptosis, pyroptosis, and necroptosis and cannot be characterized by any single form of death, often coinciding with the pathophysiological conditions of infectious diseases (Malireddi et al. 2019). PANoptosis can favor tumor suppression by inducing a mechanism whereby the host activates alternative cell death defense mechanisms ). Many PANoptosis-associated genetic markers have been identified in various tumor types, including bladder, ovarian, and colorectal cancers (Chen et al. 2021;Song et al. 2021;Ye et al. 2021), although it has been documented that PANoptosis can be a target for regulating various central nervous system diseases (Yan et al. 2022). However, there are fewer studies on the effect of PANoptosis on LGG invasive mobility and drug resistance and its role in predicting LGG prognosis. PANoptosis is a novel concept, and it is essential to explore in depth the role of these genes in LGG in terms of expression signature, predictive value, cancer signaling, CNV, and SNV.
Fortunately, large-scale public databases, represented by databases such as TCGA, CGGA, and GEO, allow us to freely access and bioinformatically analyze data from multiple groups of cancers, allowing us to gain a comprehensive understanding of the genetic landscape of cancer, examine possible biomarkers, establish treatment plans, and assess the prognosis of patients, including LGG (Tomczak et al. 2015;Barrett et al. 2011;Zhao et al. 2021). In our study, we offer a comprehensive view. We first analyzed these 66 PANGs from a pan-cancer perspective and obtained many interesting results. For example, in terms of expression characteristics, PANGs were generally upregulated in most tumors; in terms of predictive value, we found that in KIRC and SARC, PANGs mostly acted as protective genes, whereas in LGG and UVM, PANGs primarily acted as risky genes. Importantly, PANGs were significantly altered at the SNV, CNV, and methylation levels, suggesting that these genetic alterations crosstalk with each other to influence cancer signaling and thus indirectly or directly affect patient prognosis.
LGG is highly heterogeneous, which explains the low OS of patients with LGG (Lemasson et al. 2013). More than a decade has passed since the great advances in chemotherapy for LGG, namely the creation of temozolomide. That is because of the frequent occurrence of tumor resistance when traditional histological classifications guide antitumor therapy. These have highlighted the critical importance of identifying accurate molecular subtypes of LGG and developing effective biomarkers for early stratification and prophylactic treatment of high-PANscore patients with poor prognoses (Huang et al. 2020). Although previous studies have also identified several molecular subtypes of LGG, there remains considerable heterogeneity between subtypes. Therefore, there is an urgent need for a more accurate classification of LGG to improve patient survival.
We used an unsupervised consistent clustering algorithm based on the PANGs to classify LGG patients into two potential subgroups that differed significantly in terms of prognosis, immune cell infiltration, immune-related targets, and functional biological pathways. The TME plays a crucial role in distant tumor metastasis, immune evasion, and local drug resistance, and the components of TME influence the cancer immune phenotype and affect patient prognosis (Chen et al. 2015). Compared to PANcluster B, patients with PANcluster A had poorer OS; PANcluster B not only has a high level of immune infiltration but is also significantly associated with the interstitial activity. GSVA results also showed that PANcluster B was mainly enriched for regulating immune responses, positive gamma-interferon production regulatory mechanisms, and tumor necrosis factor-mediated signaling pathways, among other related pathways.
In addition, the DEGs between the two clusters allowed us to identify two genetic subtypes. The A pattern of PANcluster almost overlapped with gene cluster B. These results suggest that the two gene clusters represent different patterns of PANoptosis, with glioma cells being more prone to PANoptosis in gene cluster A. Subsequent prognostic survival analyses were performed to look for signature genes (PTGFRN, SULF2, LRRC4, CCNY, ARL3) associated with survival in LGG patients and to calculate PANscores as an indicator for stratification of LGG patients. PTGFRN is a Ig superfamily CAM that is upregulated in various tumor cells, including gliomas (Goenaga et al. 2007;Guilmain et al. 2011;Kolesnikova et al. 2009), and is an indicator of poor prognosis (Mala et al. 2022). Similar results have recently been shown for PTGFRN overexpression in gliomas, promoting cell growth and radiation resistance Fig. 9. Characterization between PAN risk score subtypes. A Correlation of immune cells with 5 signature genes. B Correlation of immune cells with PAN risk score. C Comparisons of immune scores, stromal scores, estimate scores, and tumor purity between PAN risk score subtypes. D Correlation analysis of PANscore and CSC index. E Comparisons of TMB between PANscore groups and examination of the relationship between PANscore and TMB. F The features of somatic mutations in groups with high and low PANscores are shown in a waterfall plot ◂ through PI3K-AKT signaling (Aguila et al. 2019). Another study has identified the oncogenic role of PTGFRN and its regulation by miR-137 in gliomas, thus suggesting it as a potential therapeutic target (Mala et al. 2022). SULF2 expression is increased in various human cancers, including glioma, lung, breast, head and neck, and pancreatic cancers (Wade et al. 2013;Phillips et al. 2012;Morimoto-Tomita et al. 2005;Nawroth et al. 2007), and promotes the development of gliomas (Phillips et al. 2012). LRRC4 plays an important role in early nervous system development and Fig. 10. Predicting new chemotherapy regimens. A Drug sensitivity analysis of chemotherapeutic agents between PANscore subgroups based on "oncoPredict." B Volcano map analysis based on differentially expressed genes between PAN risk score subgroups. C Top 5 potential small molecule chemotherapeutic agents for the treatment of LGG screened by CMap analysis based on the "XSum" approach. D Based on the best selected conformation of TNF and butein, the binding site and interaction results are shown differentiation, particularly during synapse formation (Kim et al. 2006;Woo et al. 2009;Xu et al. 2015). LRRC4 has also been identified as an essential component of glioma's immune microenvironment (Feng et al. 2020), which binds to PDK1 and HSP90, promotes NF-κB translocation and cytokine production in glioma cells, and influences Treg cell infiltration in the glioma microenvironment (Li et al. 2017). In glioma cells, CCNY stimulates cell proliferation, Fig. 11. Nomogram construction and validation. A, B Univariate and multifactorial analyses showed the prognostic value of the PAN risk score. C A nomogram to predict 1-, 3-, and 5-year survival in LGG patients. D Nomogram-based survival curves. E ROC curves based on nomogram predicted survival and sensitivity at 1, 3, and 5 years. F Predicted nomogram calibration curves for 1-, 3-, and 5-year OS across the cohort. G, H, I DCA curves were compared at 1, 3, and 5 years for LGG patients colony formation, and cell cycle advancement, whereas its knockdown suppresses these processes (Xu et al. 2010). As a member of the ARF family, ARL3 regulates cell morphology and cell division through microtubule-based processes (Zhou et al. 2006). It has also been shown that ARL3 expression is reduced in gliomas and correlates with tumor grade and that low ARL3 expression is associated with poor outcomes in gliomas and resistance to radiotherapy (Wang et al. 2019). Previous studies have shown that our signature gene plays an essential role in the development of glioma.
TMB for the high-PANscore group were higher than those for the low-PANscore group, according to our analysis of the mutation data for the TCGA-LGG cohort. A subsequent Spearman correlation analysis found a positive correlation between PANscore and TMB. We also discussed differences in the distribution of somatic mutations in the TCGA-LGG cohort between the two PANscore sets using maftools. The top three mutated genes in the two PANscore sets were IDH1, TPR3, and ATRX, respectively. IDH mutations in tumor cells result in a loss of function in their normal tissues, accumulation in mutated tumor cells, and ultimately, DNA or histone hypermethylation. The frequency of these mutations was significantly higher in patients in the low-PANscore group compared to the high-PANscore group. The PANscore was therefore developed as a prognostic predictor and validated. We found significant differences in prognosis, mutations, TME, clinical features, and drug sensitivity between the two PANscore subgroups of patients.
Although several genetic markers have been established for predicting prognosis in glioma patients in recent years, they are not associated with PANoptosis. Furthermore, simply classifying samples into low-or high-PANscore groups with different prognoses based on specific genetic characteristics lacks adequate biological explanation. We did this by combining the PANscore with clinical traits and demonstrating that PANscores are independent prognostic factors. To facilitate clinical practice, we created a quantitative nomogram with high AUC values (0.886, 0.855, and 0.789 at 1, 3, and 5 years, respectively) to predict gliomas' prognosis. The results of the C-index and DCA confirm the excellent performance of the columnar line graph model. By stratifying LGG based on the predictive model, we can better understand its molecular mechanism and propose new therapeutic approaches.

Conclusion
It is important to note, however, that we mostly rely on bioinformatic analyses when studying the relationship between PANGs and TME in LGG. The specific mechanisms by which PANGs affect TME await further ex vivo experimental studies, which may have important implications for the treatment of LGG. In conclusion, this study identified two PANoptosis-associated molecular subtypes in LGG, indicating different prognoses. PANGs are essential contributors to TME heterogeneity in LGG. The PANscore is a promising biomarker for identifying prognosis, molecular subtypes, TME cell infiltration characteristics, and targeted therapeutic efficacy in LGG patients.
Author contributions AA, AM: the study's design, data analysis, and article writing were all completed. NY, QF, SL, YL: checked the text and assisted with data analysis and collation. YW, QZ: reviewed and revised the manuscript. The article's submission was reviewed and approved by all authors.
Funding This work was supported by the Xinjiang Key Laboratory of Neurological Disorder Research (grant number: XJDX1711-2202).

Data availability
The original contributions presented in the study are included in the article and further enquiries can be directed to the corresponding author.

Declarations
Competing interests The authors declare no competing interests.