Machine learning-based investigation of regulated cell death for predicting prognosis and immunotherapy response in glioma patients

Glioblastoma is a highly aggressive and malignant type of brain cancer that originates from glial cells in the brain, with a median survival time of 15 months and a 5-year survival rate of less than 5%. Regulated cell death (RCD) is the autonomous and orderly cell death under genetic control, controlled by precise signaling pathways and molecularly defined effector mechanisms, modulated by pharmacological or genetic interventions, and plays a key role in maintaining homeostasis of the internal environment. The comprehensive and systemic landscape of the RCD in glioma is not fully investigated and explored. After collecting 18 RCD-related signatures from the opening literature, we comprehensively explored the RCD landscape, integrating the multi-omics data, including large-scale bulk data, single-cell level data, glioma cell lines, and proteome level data. We also provided a machine learning framework for screening the potentially therapeutic candidates. Here, based on bulk and single-cell sequencing samples, we explored RCD-related phenotypes, investigated the profile of the RCD, and developed an RCD gene pair scoring system, named RCD.GP signature, showing a reliable and robust performance in predicting the prognosis of glioblastoma. Using the machine learning framework consisting of Lasso, RSF, XgBoost, Enet, CoxBoost and Boruta, we identified seven RCD genes as potential therapeutic targets in glioma and verified that the SLC43A3 highly expressed in glioma grades and glioma cell lines through qRT-PCR. Our study provided comprehensive insights into the RCD roles in glioma, developed a robust RCD gene pair signature for predicting the prognosis of glioma patients, constructed a machine learning framework for screening the core candidates and identified the SLC43A3 as an oncogenic role and a prediction biomarker in glioblastoma.

Gliomas are the most common primary malignant tumors of the central nervous system, with approximately 10,000 new cases diagnosed each year in the United States 1 .According to the 2016 World Health Organization (WHO) classification of CNS tumors, gliomas are categorized into four grades (I-IV).Grades I and II are considered as low-grade gliomas (LGG), while grades III and IV fall under high-grade gliomas (HGG).Among these, grade IV gliomas, also known as glioblastoma multiforme (GBM), exhibit the highest degree of growth aggressiveness 2 .Unfortunately, despite advances in diagnosis and drug therapy, glioblastoma remains incurable,

Single cell RNA-seq process
The processing flow of single-cell data was referenced to the source literature of the data (https:// github.com/ TheJa ckson Labor atory/ singl ecell glioma-verha aklab, 24 ).The parameters of quality control, normalization and batch effect correction were totally same as these in the study.Out of the 55,284 single cells, a total of 12 cell clusters were obtained through cell clustering and dimension reduction.Each cluster was well-annotated using the R package "Seurat", following the guidelines from the reference literature.
To gain a deeper understanding of the underlying biological processes and cellular heterogeneity within the complex biological system, RCD pathway enrichment analysis at the single-cell level was conducted using the R package "AUCell".Furthermore, we utilized the "UMAP" function, which downscaled the data, to efficiently visualize the highdimensional profiles of gene expression or pathway activity.

Identification of RCD subtypes in glioma
The Single Sample Gene Set Enrichment Analysis (ssGSEA) was conducted using the R package "GSVA" to investigate the activated level of the 18 RCD in glioma cohorts.
Additionally, compared to k-means, pam's clustering method is insensitive to outlier data, receives fewer images of biased data, and yields classes with higher confidence.So, we performed consensus clustering analysis with the "pam" clustering method, Euclidean correlation distance, and 1000 repetitions using the R packages "ConsensusClusterPlus" and "limma", just like the previous studies 28,29 .Subsequently, we conducted Principal Component Analysis (PCA) to assess the distribution of the transcriptomic data for the identified RCD subtypes.

Functional pathway enrichment analysis
In this section, several methods were employed to investigate potential activated or suppressed signals in glioma.Based on a study of selecting an optimal method for differential expression gene analysis, for population-level RNA-seq studies with large sample sizes, the Wilcoxon rank-sum test was recommended 30 .Moreover, it also recommended the significance threshold for this analysis.So, here, to identify Differential Expression Genes (DEGs) we applied the Wilcoxon rank-sum test with a significance threshold of p < 0.05.
Gene Ontology (GO) and KEGG pathway enrichment analyses of the DEGs were conducted using the R package "ClusterProfile" 31 .The KEGG pathways were obtained from the website (www.kegg.jp/ kegg/ kegg1.html).
Additionally, we performed Gene Set Variation Analysis (GSVA) to assess the enrichment or activity of gene sets in biological samples 32 .We used the R package "GSVA" along with downloaded gene lists that included cancer hallmark gene sets, KEGG gene sets, GO biological process gene sets, and REACTOME gene sets.A p value < 0.05 was considered statistically significant for this analysis.

Multidimensional analysis of immune
The absolute percentage of immune cells was computed using the CIBORTSORT algorithm 33 .The infiltration of immune cells and stromal cells were represented by the immune score and stromal score, which were calculated by the ESTIMATE algorithm 34 .Furthermore, we compiled a set of immune checkpoint genes, including CD274, PDCD1, CD247, PDCD1LG2, CTLA4, TNFRSF9, TNFRSF4 and TLR9, based on previous literature 35 .In addition, we gathered several other signatures that describe specific features of immunity.These signatures included the anti-tumor immune cycle, tumor micro-environment and various immune traits [36][37][38] .More detailed information about these signatures is provided in supplementary table 5 and 6.

Construction of the RCD-related gene pair signature in glioma
Gene pair signatures have the potential to enhance the predictive power of computational models when compared to single gene-based signatures.The consideration of interactions between pairs of genes allows for the capture of additional information about biological processes, regulatory networks, and molecular interactions, ultimately leading to more accurate predictions.Here, we have established an applicable framework for the development of a robust RCD-related gene pair signature in glioma.The framework consists of three steps.Screening out the RCD related genes in glioma.To achieve this goal, a three-step strategy was employed.Firstly, the consensus clustering analysis was used to identify different RCD clusters.Genes with a Wilcoxon rank sum test p value < 0.05 and an absolute log2 fold change |log2FC|> 1 between distinct RCD patterns in TCGA glioma were considered as RCD-related DEGs.Secondly, the pathway activity of the 18 RCD signals was calculated.For each RCD, differently expressed gene analysis and Spearman's correlation analysis were conducted.For instance, genes related to alkaliptosis (Spearman's correlation p value < 0.05) were identified as Ga.DEGs between the high and low alkaliptosis subgroups, determined based on the median of the ssGSEA score of alkaliptosis were identified as Gb (Wilcoxon rank sum test p value < 0.05).Genes shared between Ga and Gb of alkaliptosis in the TCGA-GBM Hiseq dataset were considered as G1.The same process was followed for the rest of the RCDs, resulting in G1-G18.The geometric mean of Spearman's R was then computed for each gene across G1-G18 using the formula: where x1, x2, x3, …, xn represent the Spearman's correlations.Genes with a geometric mean of Spearman's R > 0.3 were considered as C.genes.Finally, the final RCD-related genes were obtained by taking the intersection of RCD-related DEGs from the first step and C.genes from the second step.
Identifying the prognostic RCD-related gene pairs.After obtaining the final RCD related genes, the gene pair matrix was constructed using the previously identified RCD-related genes.In this matrix, Gene x (Gx) and Gene y (Gy) represented the RCD-related genes in glioma.The expressions of Gx and Gy were denoted as Ex and Ey, respectively.The score of the gene pair (Gx, Gy) was calculated as follows: In both the TCGA GBM cohort and the TCGA glioma cohort, univariate Cox regression analysis was conducted for the gene pairs that were constructed.Gene pairs with p-values less than 0.05 in both cohorts, as determined by the univariate Cox regression analysis, were retained to form the final scoring system.Developing the RCD related gene pair scoring system.For each optimal RCD gene pair, RCD.GP (i) was assigned a value of either 1 or 0. The risk score of the sample was then determined based on all the optimal gene pairs.Specifically, the scoring system involved calculating the sum of the scores of the gene pairs.The RCD gene pair risk score for each patient or cell was computed as follows:

Survival analysis
After scoring all patients in the cohorts with the RCD.GPscore system, we divided the patients into high and low RCD.GP score subgroups based on the maximally selected rank statistics through the function "surv_cutpoint" in the R package "survminer".The log-rank sum test was conducted to compare the statistical difference between the distinct subgroups.For the univariate and multivariate Cox regression analysis, we used the R package "survival".

Immunotherapy related analysis
Several methods were used to estimate the value of the RCD.GP score for immunotherapy.We collected a number of signatures related to immunotherapy response from previous literatures, including T cell-inflamed GEP, IFNG, CD274, CD8, TLS-melanoma, TLS, T cell dysfunction, MDSC, TAM M2, T cell exclusion and CAF [39][40][41] .Some researchers have shown that these signatures are closely related to the response to immunotherapy and that they can predict the response to immunotherapy 42,43 .So, the relationship between the RCD.GP score and the immunotherapy response could be evaluated by computing the link between the RCD.GP score and these signatures, just like the previous studies did 44,45 .We calculated the Spearman's correlation between the GSVA score of these signatures and the RCD.GP score.Additionally, we compared the difference between the high and low RCD.GP score subgroups for each of these signatures.

Screening out the potential therapeutic targets in glioma
Machine learning is well-known and has been used in previous bioinformatics studies 59,60 .Although different machine learning algorithms have different optimal usage scenarios and different usage characteristics 60 .XgBoost is a gradient boosting algorithm used to solve regression and classification problems in machine learning 60 .It also performs feature selection by assigning a score to each feature based on its importance in the model.Various machine learning algorithms have been used for feature selection to screen for more central genes 61,62 .Here, with some machine learning algorithms, we present a computational framework for identifying potential therapeutic targets for glioma treatment.The framework involves four main steps.
-Identification of prognostic genes in glioma.We performed univariate Cox regression analysis in three cohorts: TCGA LGG, TCGA GBM and TCGA LGG-GBM cohort.Genes with a p value less than 0.05 in all three cohorts were considered as prognostic genes in glioma.-Selection of RCD-related genes.RCD-related genes were filtered out by identifying the intersection of RCDrelated DEGs between the distinct RCD patterns and the previously extracted C.genes.-Screening important genes in glioma using six machine learning methods.Six machine learning methods including CoxBoost, Boruta, random survival forest (RSF), least absolute shrinkage and selection operator (Lasso), eXtreme Gradient Boosting (XgBoost), and Elastic net (Enet) were employed.For the elastic net method, we used different mixing parameters (α) values ranging from 0.1 to 0.9.Genes identified by at least ten methods with different parameters in both TCGA GBM and TCGA GBMLGG cohorts were considered potential therapeutic genes in glioma.

Cell culture
The human astrocytes (HA1800) and three GBM cell lines, including A172, LN229 and U87, were obtained from the Cancer Research Institute of Central South University.The GBM cell lines were cultured in Dulbecco's Modified Eagle Medium (DMEM)(Gibco) supplemented with 10% fetal bovine serum (FBS).The HA1800 was cultured with the astrocyte medium which was purchased from ScienCell Research Laboratories (Carlsbad, CA, USA).The FBS was from the Thermo Fisher Scientific.All cells were maintained at 37 °C with 5% CO2.

Collection of glioma tissue sample
Clinical glioma samples and normal brain tissues used in this study were obtained from patients surgically treated at the Neurosurgery Department of Xiangya Hospital of Central South University from September 2019 to August 2023, of which WHO Grade II (n = 12), Grade III (n = 12), Grade IV (n = 12), and normal brain tissues (n = 12).Normal brain tissues were normal tissues that had to be removed during glioma surgery due to the need for surgical access as well as book exposure of the lesion.Detailed information about the tissue samples was summarized in supplementary table 8. Informed agreement was obtained from all patients.This study was performed with ethical approval of the ethics committee of Xiangya Hospital in accordance with the Declaration of Helsinki and written informed consents were obtained from all of the enrolled subjects.

qRT-PCR
Based on the TRIzol (Accurate Biology, China, AG21101) and RNA extraction kit (Thermo Scientific, K0731), we extracted RNA from clinical glioma samples, which was then stored at -80 °C.The extracted RNA was utilized to generate cDNA through amplification, which was subsequently used for qRT-PCR.To perform qRT-PCR, we queried the gene sequence in GenBank and designed the primer sequences for the target genes.The relative levels of indicated genes were analyzed through the 2 − ΔΔCt method.The primer sequences of BMP2, IGFBP2, SPP1, SLC43A3, P2RY6, PTX3, and STEAP1, GAPDH, were exhibited in supplementary table 8. Two-tailed Student's t test was performed for comparison between two groups of samples.

Statistical analysis
All analysis was performed in R (version 4.1.3,http:// www.rproj ect.org/).For the continuous variables and the categorical variables, the Wilcoxon rank sum test and the chi-squared test were adopted, respectively.Correlation analysis between the different variables was performed using Spearman's coefficients.The Benjamini and Hochberg method was used to adjust the P value.P value less than 0.05 was considered to be statistically significant in this study.The graphic abstract was provided in Fig. 1.

Ethics approval and consent to participate
The study was approved by the ethics committee of Xiangya Hospital, and the written informed consent was obtained from all patients.

Dysfunction of regulated cell death in glioma
Dysfunction of regulated cell death in tumors, specifically the failure of programmed cell death mechanisms, plays a significant role in tumor development and progression 11 .In normal physiological conditions, cells undergo programmed cell death, including apoptosis, autophagy, and necroptosis, to maintain tissue homeostasis and eliminate damaged or abnormal cells.Here, we comprehensively investigated the RCD level between glioma and normal brain cortex with collected 18 RCD signatures.It might seem counterintuitive that the level of RCD was dysregulated and higher in gliomas compared to normal tissue (Fig. 2A).With the Spearman's correlation analysis, we also estimated the inner regulated network among these RCD in glioma and the normal tissue (Fig. 2B).The alteration of the correlation also provided solid evidence of the dysregulated RCD in glioma.This abnormal increase and dysregulation of RCD are often detrimental to the survival of patients with many types of tumors, as reported in various publications [63][64][65] .
To confirm whether this also applied to gliomas, we performed univariate Cox regression in several cohorts of LGG, GBM and glioma.The results revealed that the high level of the RCD was a risky factor for glioma patients (Fig. 2C).The transcriptomic profiles from the CCLE project indicated that the level of most RCDs was higher in most high grade glioma (HGG) cells than in LGG cells (H4) (Fig. 2D).From previous publications, we selected some core genes in individual RCDs, and immunohistochemically stained tissue sections from the HPA of these genes showed that the level of RCD was higher in glioma than in normal brain tissue (Fig. 2E) www.nature.com/scientificreports/Collectively, at both the mRNA and protein levels, we found a significant increase in these 18 already reported regulated modes of cell death and functional abnormalities.It is important to note that while the level of regulated cell death might be higher in gliomas, the effectiveness of these cell death mechanisms could be impaired or evaded by tumor cells, contributing to the survival and progression of gliomas.

RCD-based patterns show distinct micro-environments
According to the consensus clustering analysis based on the profile of the RCD, the patients in the TCGA glioma dataset could be divided into two distinct RCD patterns, named RCD cluster A (n = 220) and RCD cluster B (n = 469) (Fig. 3A, supplementary table 9).The PCA results indicated different characteristics of the RCD profiles between the two RCD clusters (Fig. 3B).Next, we described the differences and connections between these two RCD clusters at multiple levels, including clinical features, immune microenvironment, signaling pathways and so on.
The patients in RCD cluster A had significantly better overall survival than patients in RCD cluster B (Fig. 3C, log-rank test p < 0.001, HR: 3.195, 95%CL: 2.46-4.13).Compared with RCD cluster A, RCD cluster B had more patients with un-methylated MGMT promoter, chr 19/20 co-gain, wild type IDH, higher-level pathological tissue types and older ages (Fig. 3D).The levels of most RCDs were higher in RCD cluster A than in RCD cluster  The ssGSEA score of the RCD signatures in glioma cell lines.The red represented that the score in this GBM cell line was higher than in LGG cell line (H4), while the white represented the opposite.The right panel showing the number of the GBM cell lines in which the score was higher than in the H4.(E) The immunohistochemically stained tissue sections images from the HPA of the core RCD genes indicated that the level of the RCD was different between in glioma and in normal brain tissues.The detailed clinical information of the images was provided in supplementary table 3. www.nature.com/scientificreports/B (Fig. 3D), which was consistent with previous results showing that higher RCD levels were detrimental to patients' prognosis.The expression of some classical immune checkpoint genes, including CD274, PDCD1, CD247, PDCD1LG2, CTLA4, TNFRSF9, TNFRSF4 and TLR9, was investigated in the RCD clusters, and the results indicated that these genes, except for TLR9, had higher expression in RCD cluster A. Both innate and adaptive immune cells infiltrated at significantly higher levels in RCD cluster A than in RCD cluster B (Fig. 3D).The TMB increased in RCD cluster A (Supplementary Fig. 1A).The pathway enrichment analysis based on the DEGs between the two RCD clusters revealed that the dysfunction and abnormality of the RCD were significantly related to neutrophil chemotaxis, monocyte chemotaxis, cytokine-cytokine receptor interaction, and collagen degradation (Fig. 3E).
We also performed the same analysis method in the TCGA GBM cohort and obtained similar results to those in the TCGA glioma cohort (Supplementary Fig. 1B, C, supplementary table 10).Generally speaking, in the context of tumor development, dysregulation of regulated cell death pathways could impact both tumor cells and immune cells, leading to alterations in the immune microenvironment 69 .Some forms of regulated cell death, such as necroptosis and pyroptosis, could induce inflammation and release danger signals (damage-associated molecular patterns-DAMPs) 70,71 .DAMPs can activate innate immune cells and promote an inflammatory response 72 .This inflammation can influence the recruitment and activation of various immune cells within the tumor microenvironment, ultimately affecting the prognosis of glioma patients.

Integrated single-cell level analysis of the RCD in glioblastoma
Based on the previous results indicating dysfunction and an abnormal increase in RCD in glioblastoma with bulk-level data, we subsequently investigated RCD using single-cell data to enable a finer scale assessment of the RCD landscape in gliomas.After clustering, dimensionality reduction, and annotation, we present the 12 types of cells and their corresponding numbers (Fig. 4A).For each single cell, we estimated the signal pathway activity of RCD using the AUCell method (Supplementary Fig. 2).The pathway activity of a cell type was represented by the average of all cells of the same type.Disulfidptosis and immunogenic cell death had higher AUCell scores compared to other RCD types (Fig. 4B-D).Notably, disulfidptosis exhibited higher signal levels in endothelial and pericyte cells than in other cell types.The dominant RCD type of each single cell was defined as the type of RCD with the highest level, and we observed a significant proportion of cells with disulfidptosis and immunogenic cell death as the dominant RCD (Fig. 4E).
Given the previous results showing an absolute predominance of disulfidptosis in glioblastoma, we investigated the core disulfidptosis genes in different cell types.SLC7A11 was highly expressed in fibroblast cells, while MYH10 showed high expression levels in malignant cells, especially in Prolif.stem-likecells (Fig. 4F).We found that mitotic cell death was the dominant RCD type in almost all the Prolif.stem-likecells.The UMAP view of cells with the dominant RCD (top) and cell density (bottom) displayed the RCD profile distribution across the different cell types (Fig. 4G).Furthermore, we investigated the substantial changes in the different RCD cell landscape when exposed to external intervening factors including hypoxia and irradiation (Fig. 4H-K, supplementary Figs. 3, 4).
As the duration of hypoxia increased, the proportion of cells in which disulfidptosis predominated gradually increased, reaching levels similar to controls, while immunogenic cell death gradually decreased to baseline levels (Fig. 4J).Overall, we explored the unique landscape of RCD in glioblastomas from a single-cell perspective and found that disulfidptosis and immunogenic cell death appeared to play more important roles in glioblastomas than other RCD types.Additionally, it is worth noting that in the presence of external environmental stimuli, the transformation of these two modes of cell death may be a potentially possible mechanism used by cells to adapt to external stresses.

Development and validation of a novel RCD-related gene pair signature
The prognosis of glioblastoma is poor, with an overall survival rate remaining relatively low.Even with aggressive treatment, the median survival is typically around 12-16 months 73 .The results presented above suggested that RCD had a significant impact on patients with glioblastoma.Therefore, the purpose of this study is to establish a reliable RCD-related model that can be used in the clinic to accurately predict a patient's prognosis, quality of life, and response to treatment.This model aims to assist doctors in achieving accurate treatment for their patients.The framework for constructing the scoring system is displayed in Fig. 5A.
With calculated geometric spearman's correlation, we identified 6162 C.Genes related RCD.Detailed descriptions of this process are provided in the methods section (Fig. 5B).Additionally, 8084 genes with a p-value less than 0.05 and |log2FC| greater than 1 were identified as DEGs between distinct RCD patterns (Fig. 5C).The intersection of C.Genes and DEGs resulted in the final set of RCD-related genes (n = 725).After removing missing genes in most validation datasets, we constructed 639 RCD-related gene pairs.Among the 203,841 gene pairs based on these 639 genes, 377 gene pairs in TCGA-GBM and 40,788 gene pairs were identified as prognostic RCD-related gene pairs using univariate Cox regression (Fig. 5D, supplementary Fig. 5A, p < 0.05).Finally, 118 gene pairs, including 102 genes, were selected as the optimal gene pairs for constructing the scoring system, named RCD.GP score (Supplementary table 11).

Glioma with high RCD.GP score possesses strongly malevolent biology and activated immune characteristics
The tumour microenvironment plays a crucial role in tumor growth, immune surveillance, immune escape, and response to therapy 74 .Immune checkpoints are molecules on immune cells that regulate the immune response, preventing excessive activation and tissue damage 30 .However, tumors can hijack these checkpoints to suppress immune responses and avoid immune destruction.Inhibitory immune checkpoint molecules, such as PD-1, PD-L1, and CTLA-4, are frequently expressed within the tumor microenvironment and can dampen immune responses against cancer cells.Our investigation revealed that glioblastomas with high RCD.GP scores highly expressed immune checkpoint genes, including CD274, CD247 and so on (Fig. 6A).As the immune microenvironment consists of various immune cells, including lymphocytes, macrophages, dendritic cells, and myeloidderived suppressor cells (MDSCs) which can infiltrate the tumor site and interact with cancer cells, influencing tumor progression, we assessed the infiltration of the immune cells, and found that higher RCD.GP risk scores correlated with increased infiltration of cells such as M0 macrophages, neutrophils, and Tregs (Fig. 6A).Similar results were observed in the CGGA dataset (Supplementary Fig. 6A).
To further confirm the association of RCD.GP scores with malignant biology and tumor-associated immunity modulation in glioblastoma, we performed GSVA with the cancer hallmark signature database, KEGG database, REACTOME database and GO BP database in TCGA-Glioma.The high RCD.GP scores were associated with activated malignant pathways and immune-related signals, such as epithelial mesenchymal transition, complement and coagulation cascades, cytokine-cytokine receptor interaction, and complement cascade (Fig. 6B, C, supplementary Fig. B-D).The enrichment analysis with GO, KEGG and REACTOME pathways further validated the activation of immune-related signals and malignant pathways in the high RCD.GP score subgroup (Fig. 6D, supplementary Fig. 6E, F).The anti-tumor immune cycle, representing the steps and interactions involved in mounting an effective immune response against tumors 75 , exhibited higher signals in the high risk score subgroup, indicating an activated anti-tumor immune response (Fig. 6E).
We also estimated the correlation between a number of immune-related scores 76 (n = 92) and RCD.GP scores, and the results indicated a positive relationship between most immune-related scores and the RCD.GP score (Supplementary Fig. 6G).Although patients with higher risk scores had worse prognoses, we investigated the features of the tumor microenvironment using signatures from the cancer single-cell states atlas 36 (n = 14) and the previous literatures (n = 29) 77 .The GSVA results revealed that factors promoting intense malignant progression of tumors, such as EMT, inflammation, metastasis, immune suppression by myeloid cells, pro-tumor cytokines, and tumor proliferation, were significantly higher in the high-risk group than in the low-risk group (Fig. 6F-G, supplementary Fig. 6H).One possible explanation for the apparent contradiction is that the tumor microenvironment, characterized by factors such as inflammation, immune cell infiltration, hypoxia, and nutrient deprivation, can influence the levels of cell death in tumors.Certain aspects of the tumor microenvironment could promote cell death, while others could protect tumor cells from undergoing programmed cell death.In glioblastoma, both might be at a high level, leading to an imbalance that favors increased cell death but does not necessarily result in tumour shrinkage.This imbalance could be due to the rapid and uncontrolled proliferation of tumor cells, leading to an increased number of cells requiring elimination.

Assessment of the capability of the RCD.GP score for the immunotherapy response
Considering the impressive association of the RCD.GP score with immune-related characteristics, we hypothesized that the RCD.GP score could be strongly linked to the response to immunotherapy.To investigate this, we collected some signatures related to immunotherapy response and found that the RCD.GP score positively correlated with the GSVA scores of these signatures (Fig. 7A).
Previous studies have demonstrated that tumors with high TMB are more likely to respond to immunotherapy due to enhanced recognition and targeting by the immune system 42,78 .We observed a positive correlation between the RCD.GP score and the TMB (Fig. 7B, R = 0.39, p < 0.001).Furthermore, the TMB was higher in the high-risk score group compared to the low-risk score group (Fig. 7C, Wilcoxon rank sum test p < 0.001).Other immunotherapy response-related scores, such as TCR richness, TGF-beta response, and lymphocyte infiltration signature score, also showed positive associations with the RCD.GP score (Fig. 7D, supplementary table 13).In the highrisk score group, most of these scores were higher than those in the low-risk score group (Supplementary Fig. 7).
To further verify the predictive efficacy of the RCD.GP score in immunotherapy response, we analyzed multiple cohorts with immunotherapy treatment.Strikingly, the higher RCD.GP score group exhibited better prognosis and immunotherapy response in patients with UC, NSCLC and SKCM, while intriguingly, the K-M curves exhibited the opposite trend in patients with glioblastoma (Fig. 7E-J, supplementary table 14).Overall, our study suggests that patients with a high RCD.GP score may have a higher potential to benefit from immunotherapy treatment.

Identification of SLC43A3 as novel and potential RCD-related oncology target
In this study, we present a novel framework for screening important features (Fig. 8A).Differentially expressed gene analysis was performed among GBM, LGG and normal brain cortex (Supplementary Fig. 8A-C), leading  www.nature.com/scientificreports/ to the identification of 3260 genes as DEGs across all groups (Fig. 8B).Among these, 725 genes were identified as RCD-related genes in glioma, and 1660 genes were associated with the prognosis of the glioma patients (Fig. 8A, B).From this analysis, 60 genes were identified as RCD-related prognostic DEGs (Fig. 8B, supplementary Fig. 8D).Using a machine learning framework on TCGA-Glioma cohort and TCGA-GBM cohort, we generated some candidate genes (Fig. 8C, D).The genes that were repeatedly screened out more than 10 times in both TCGA-GBM and TCGA-Glioma cohorts were identified as the most important genes for glioma patients (Fig. 8E, supplementary table 15).These genes included BMP2, IGFBP2, SPP1, SLC43A3, P2RY6, PTX3, and STEAP1.Among these, the role of SLC43A3 in tumors, especially gliomas, had hardly been reported, while the functions of the other genes in tumors were widely documented [79][80][81][82][83][84] .Interestingly, the expression of SLC43A3 increased with the pathological grade and was most highly expressed in the ME transcriptional subgroup, which has been associated with a malignant process and shorter median survival time in glioblastoma patients (Fig. 8F-I).We also comprehensively investigated the distribution of SLC43A3 in different clinical subgroups, revealing a potential relationship between SLC43A3 expression and poor clinical prognosis (Supplementary Fig. 9A-M).
To validate the findings, we assessed SLC43A3 expression in high-grade gliomas compared to normal brain cortex using qRT-PCR (Fig. 8J).Transcriptomic data from the CCLE project showed higher mRNA levels of SLC43A3 in most glioblastoma cells compared to LGG cells (H4) (Supplementary Fig. 9N).The qRT-PCR experiment with glioblastoma cell lines further confirmed higher expression of SLC43A3 in glioblastoma cell lines compared to astrocyte cell lines (Fig. 8K).Prognostic analysis using multiple glioma datasets revealed that high expression of SLC43A3 was associated with worse overall survival (Supplementary Fig. 10A-C).Meta-analysis and multiple variate Cox regression analysis provided strong evidence that SLC43A3 was a risky and independent prognostic factor for predicting glioma patients' outcomes (Fig. 8L, M).
In conclusion, we have provided a novel framework for screening important features, identified seven potential important therapeutic targets, and presented solid evidence supporting an oncogenic role of RCD-related SLC43A3 in glioblastoma.

Discussion
Glioblastoma is the most common and deadliest type of brain tumor in adults, accounting for approximately 15% of all primary brain tumors.It originates from glial cells in the brain, and is characterized by uncontrolled cell growth, invasion into surrounding tissues, and resistance to cell death.Despite advancements in treatment, the prognosis for glioblastoma remains poor, with a median survival time of about 12-16 months.Its highly invasive nature and resistance to therapies make it challenging to treat, emphasizing the urgent need for more accurate predictive tools to improve overall survival rates, symptom management, and patients' quality of life.
Previous predictive signatures for glioblastoma patient prognosis mainly relied on individual gene expression, disregarding the interactions or synergistic relationships between genes.In contrast, predictive models based on gene pairs offer several advantages including increased predictive power, enhanced interpretability and transferability to different datasets or experimental conditions 85 .Here, we comprehensively investigated the RCD landscape of glioblastoma from both bulk and single-cell perspectives.We developed a robust and accurate RCD-related gene pair signature that holds potential for transferability across different datasets.Additionally, we provided a framework for screening out relatively core genes and identified SCL43A3 as a potential therapeutic oncology target.
Our results revealed a complex and multifaceted relationship between RCD and glioblastoma.The elevated levels of almost all RCD in glioblastoma were associated with worse prognosis, strongly correlated with malignant biology processes and immune microenvironment dysfunction.Cancer cells often exploit immune checkpoint pathways to evade immune surveillance and promote tumor growth 86 .They upregulate immune checkpoint molecules on their surface or within the tumor microenvironment, leading to the suppression of anti-tumor immune responses.Glioma with higher levels of RCD exhibited increased expression of immune checkpoint genes, indicating potential suppression of anti-tumour immunity.Combining therapies that induce regulated cell death with immune checkpoint inhibitors have shown promising results in preclinical and clinical studies [87][88][89][90][91] .Inducing immunogenic forms of cell death can enhance the immunogenicity of tumors and improve the efficacy  www.nature.com/scientificreports/ of immune checkpoint blockade by promoting the activation of anti-tumor immune responses 92 .Signaling pathways promoting malignant biological processes such as JAK-STAT signaling pathway, TNF signaling pathway and NF-kappa B signaling pathway, were demonstrated to be extraordinarily related with the high level of RCD.However, we also observed enhanced anti-tumor immune responses, such as the increased infiltration of antitumor immune cells and the activated anti-tumor immune signals like IL-17 signaling pathway, cytokine-cytokine receptor interaction and neutrophil chemotaxis which might be a result of the dysregulated tumor microenvironment promoting both cell death and tumor protection mechanisms.Further investigation is required to understand the complexity of molecular mechanisms underlying this association and its therapeutic implications.Single-cell analysis suggested a significant role for disulfidptosis in glioblastoma.Disulfidoptosis, also known as disulfide bond deficiency, is a condition characterized by an impaired formation or maintenance of disulfide bonds 68 .Disulfide bonds are important for the stability and proper folding of proteins within cells.Disulfide bond deficiency can potentially impact various cellular processes that are relevant to tumorigenesis.One possible link is the role of disulfide bonds in protein folding and cell signaling.Proteins involved in cell growth regulation, regulated cell death, and DNA repair often require proper disulfide bond formation to function correctly.Disruptions in these processes can contribute to the development and progression of tumors.The expression of SLC7A11, a core gene of disulfidoptosis, was highly expressed in fibroblast cells, while MYH10, another core gene, showed high expression in malignant cells, particularly Prolif.stem-likecells.This suggests MYH10 as a potential therapeutic target for glioma.Nonetheless, the relationship between disulfidoptosis and glioblastoma is complex, requiring further investigation.
The RCD.GP score proved to be a robust and promising biomarker for predicting clinical outcomes and immunotherapy response in glioma patients.Traditional prediction models based on absolute gene expression values may be influenced by noise and biological variability, leading to less reliable predictions.In contrast, our model, incorporating gene pairs, demonstrates improved dtability and robustness by accounting for variability in gene expression data.In 16 cohorts of glioma patients, the RCD.GP score consistently exhibited superior reliability and generalizability after validation and optimization.The validation metrics including, K-M curves, C-index, 1-year AUC, 3-year AUC, 5-year AUC, meta-analysis results, and multivariate Cox regression analysis, unequivocally support the enhanced predictive power of the RCD.GP score system.
Patients with glioma and high RCD.GP scores experienced shorter overall survival times, which were strongly associated with malignant biological processes, including coagulation, epithelial-mesenchymal transition, and inflammatory response.Additionally, the subgroups with high RCD.GP risk scores displayed an increased profile of immune checkpoint genes, enhanced infiltration of the immune cells, and elevated levels of pro-tumor signaling pathways.Immunotherapy has demonstrated promising results in various cancers, such as melanoma, lung cancer, bladder cancer, and some types of lymphomas and leukemia [93][94][95][96][97] .Our study also revealed a significantly positive correlation between the RCD.GP score and immunotherapy response indexes, including CD8, MDSC, TLS, TMB and TGF-beta response.
Interestingly, while patients with high RCD.GP scores showed improved OS and enhanced immunotherapy response in various cancer types like UC, NSCLC and SKCM, a different trend was observed in GBM patients treated with anti-PD-1 therapy.The group with higher RCD.GP scores exhibited a worse prognosis and impaired immunotherapy response.Despite the promising results of immunotherapy in treating other cancers, immune checkpoint inhibitors have shown limited efficacy in glioblastoma compared to traditional treatment modalities such as radiotherapy, chemotherapy and surgery 98 .This may be attributed to the intricate interplay between multiple immune-related mechanisms affecting glioblastoma progression and the influence of RCD in promoting malignant progression.As such, exploring combination approaches, such as using different immunotherapies in conjunction with RCD inhibition, chemotherapy, radiation therapy, or targeted therapies, holds potential for further improving treatment outcomes in glioblastoma patients.
Seven genes including BMP2, IGFBP2, SPP1, SLC43A3, P2RY6, PTX3 and STEAP1 were identified as potential therapeutic targets through a machine learning framework designed to identify the most important features.Each of these genes plays a unique role in tumorigenesis and tumor progression.BMP2 is a multifaceted gene that can exhibit both tumor-promoting and tumor-suppressive effects in different types of cancers.It can stimulate cell proliferation, angiogenesis, and metastasis in certain contexts, while also inducing cell cycle arrest, apoptosis (programmed cell death), and differentiation of cancer cells, leading to tumor growth inhibition 79 .IGFBP2 has been implicated in promoting tumor growth and progression in several types of cancer 80 .High expression of IGFBP2 can help tumor macrophages to form an immunosuppressive microenvironment and exert a substantial inhibitory effect on T cell proliferation and activation, and this may be related to immunogenic death 99 .SPP1, a glycoprotein with diverse functions, is involved in various biological processes such as cell adhesion, migration, immune regulation, and tissue remodeling.Its ability to stimulate cell proliferation, survival, and angiogenesis, facilitating tumor formation 81,100 .Previous studies have shown that SPP1 controls UPR and ER stress-induced autophagy by regulating intracellular sphingosine-1-phosphate homeostasis 101 .P2RY6 and STEAP1 are considered tumor suppressor genes in certain cancer types 84,102 .Previous findings suggest that physiological P2RY6 ligands and specific P2RY6 agonists can restore normal monocyte differentiation by restoring autophagy in some primary myeloid cells of patients with chronic granulocytic leukemia, demonstrating a potential link between P2RY6 and programmed cell death 103 .There are some evidence suggesting that silencing the STEAP1 can induce the apoptosis of the LNCaP, which indicates the potential relationship between the STEAP1 and the regulated cell death 104 .PTX3 is known to enhance tumor cell proliferation, survival, invasiveness, and angiogenesis by modulating various signaling pathways in different tumor types, including breast cancer, lung cancer, ovarian cancer, and glioblastoma 105 .The key role of the PTX3 in regulating the ferritinophagy in glioma has been recently reported, and PTX3-deficient IDH1 mutant gliomas shown enhanced autophagic signature 106 .SLC43A3, a protein-coding gene that encodes a transporter involved in amino acid transport, has been associated with fatty acid flux, nucleotide metabolism and DNA repair [107][108][109] .DNA repair is one of the most important factors for the cell survival and cell death, including apoptosis, necrosis and autophagy 110 .This suggests in part that SLC43A3 may be essential for the cell death.Moreover, its role in cancer development and progression is not yet fully understood.Our focused investigation on SLC43A3 revealed a RCD-related oncogenic potential in glioma through systematic analysis of its expression profile and prognostic value in multiple datasets and qRT-PCR experiments.As for the other six genes, with qRT-PCR, we also validated the expression in GBM cell (Supplementary Fig. 11).Overall, the identification of these seven genes as potential therapeutic targets opens up new avenues for future research and the development of targeted therapies for glioma treatment.Further studies are warranted to elucidate the precise molecular mechanisms and potential.
We acknowledge that this study has some limitations.Firstly, the cohorts with immunotherapy included only four cancer types, and to draw a more accurate and broadly generalized conclusion, data from a wider variety of cancer cohorts with immunotherapy should be included in the analysis.Secondly, the single-cell analysis revealed the core role of disulfidoptosis in glioblastoma, but only a rudimentary analysis was conducted.To comprehensively elaborate its role, further investigations involving more systematic, comprehensive, and advanced analysis and experiments are necessary.Thirdly, although we identified seven candidate genes for therapy, we only validated the expression profile of SLC43A3 in tissues and cancer cell lines using qRT-PCR.To a deeper understanding of the potential and insightful mechanisms of SLC43A3, additional experiments need to be conducted.Fourthly, the specific mechanism and links between the identified genes and the regulated cell death patterns in gliomas should be validated and elaborated with more advanced and complex net-experiments.Finally, this study utilized retrospective data without prospective clinical trials data for validating the superior reliability and generalizability of the RCD.GP score.To further validate the performance of the RCD.GP score, prospective clinical trials data should be incorporated in future studies.Moreover, the study utilized data from multiple sources, including TCGA, CGGA, and GEO.The potential variations, biases, and quality issues among these datasets might exist.Differences in data collection methods, experimental conditions, and patient demographics could introduce confounding factors.

Conclusions
In conclusion, we have conducted a comprehensive investigation of the RCD landscape in glioblastoma, exploring both bulk and single-cell aspects.Our study resulted in the development of a robust and accurate RCDrelated gene pair signature, which holds promise for potential clinical applications in predicting the prognosis of glioblastoma.The glioma patients with low RCD.GP score had better prognosis.Moreover, we established a framework consisting of Lasso, RSF, XgBoost, Enet, CoxBoost and Boruta, for identifying relatively core genes and successfully identified SCL43A3 as a potential therapeutic target in oncology based on the bioinformatics and the qRT-PCR.Furthermore, our findings highlight the significance of the RCD.GP score as a predictor for adverse clinical outcomes and impaired immunotherapy response in glioblastoma patients.These insights have the potential to improve patient management and treatment decisions.

Figure 1 .
Figure 1.The graphic abstract of this study.

Figure 2 .
Figure 2. Dysfunction of regulated cell death in glioma.(A) The comparison of the ssGSEA score of the 18 RCD signatures between gliomas and normal brain cortex, with the large-scale bulk transcriptomic data from the TCGA and GTEx project.The Wilcoxon rank sum test was performed.The two sided p value < 0.001 was represented by "***".(B) The spearman's correlation of the 18 RCD ssGSEA scores in TCGA-GBM dataset and GTEx normal brain cortex dataset.(C) The univariate Cox regression result of the 18 RCD signatures in datasets with LGG, datasets with GBM and datasets with glioma.The p value < 0.05 was considered as significance.(D)The ssGSEA score of the RCD signatures in glioma cell lines.The red represented that the score in this GBM cell line was higher than in LGG cell line (H4), while the white represented the opposite.The right panel showing the number of the GBM cell lines in which the score was higher than in the H4.(E) The immunohistochemically stained tissue sections images from the HPA of the core RCD genes indicated that the level of the RCD was different between in glioma and in normal brain tissues.The detailed clinical information of the images was provided in supplementary table 3.

Figure 3 .
Figure 3. RCD-based patterns show distinct micro-environments.(A) The glioma patients in TCGA were classified into two RCD clusters, named RCD cluster A and RCD cluster B, based on the consensus clustering analysis.(B) The PCA analysis showing the different distribution of the RCD profile.(C) The K-M curves showing that the patients in RCD cluster A had worse OS than patients in RCD cluster B. (D) A heat map showing the clinical index, RCD profile, expression of immune checkpoint genes, immune score, stromal score and tumor microenvironment in the TCGA glioma cohort.The Wilcoxon rank sum test or the chi-square test was performed the assess the difference between the RCD cluster A and RCD cluster B. "*", "**", "***", and "****" represented that the p value < 0.05, 0.01, 0.001, and 0.0001.(E) The enrichment analysis of the DEGs between the two RCD clusters.KEGG, GO BP and REACTOME databases were included for the functional pathway analysis.

Figure 4 .Figure 5 .
Figure 4. Integrated single-cell level analysis of the RCD in glioblastoma.(A) UMAP plot showing that the 12 types of cells and their corresponding numbers were obtained after clustering, dimensionality reduction, and annotation.(B) The AUCell of the 18 RCD signatures in the 12 cell types.The pathway activity of a cell type was represented by the average of all cells of the same type.(C) The ridge plot showing the distribution of the RCD profile in the single-cell level.(D) Percentage of the AUCell scores of the different RCDs in 12 cell types.(E) For a single cell, we defined the dominant RCD type of this cell as the type of RCD that has the highest level of its RCD.This plot showing the percentage of the dominant RCD types in 12 cell types.(F) Core gene expression of disulfidptosis across defined cell clusters.Bubble size is proportional to the percentage of cells expressing a gene and color intensity is proportional to average scaled gene expression.(G) UMAP view of dominant RCD (top) and cell density (bottom) displaying the RCD profile distribution across the different cell types.High relative cell density is shown as bright magma.(H) UMAP view of cell types with different stress interventions.(I) The percentage of the RCD in different cell types and in different stress interventions.(J) The percentage of the dominant RCD in different cell types and in different stress interventions.(K) The profile of the 18 RCD in different conditions and cell types.

Figure 6 .
Figure 6.Glioma with high RCD.GP score possesses strongly malevolent biology and activated immune characteristics.(A) The heat map of the relationship of the RCD.GP score and clinical indexes, immune checkpoint genes, infiltration of the immune cells and immune microenvironment function in TCGA glioma cohort.(B, C) The spearman's correlation between the RCD.GP score and the GSVA score of the tumor hallmark signatures (B) and functional pathway signatures in KEGG database (C).(D) The pathway enrichment of the GO biological process.(E, F) The plot showing the relationship between the RCD.GP score and the antitumor cycle (E), and single cell state from the cancerSEA (F).The spearman's correlation and the Wilcoxon rank sum test were performed."*", "**", "***", and "****" represented that the p value < 0.05, 0.01, 0.001, and 0.0001.(G) The spearman's correlation between the RCD.GP score and the 29 signatures related with the tumor microenvironment.

Figure 7 .
Figure 7. Assessment of the capability of the RCD.GP score for the immunotherapy response.(A) The relationship of the RCD.GP score and the GSVA score of the signatures related with the immunotherapy response.The spearman's correlation analysis and the Wilcoxon rank sum test were performed."*", "**", "***", and "****" represented that the p value < 0.05, 0.01, 0.001, and 0.0001.(B) The spearman's correlation between the RCD.GP score and the tumor mutation burden.(C) The difference of the TMB between the high RCD.GP score subgroup and low RCD.GP score subgroup.The Wilcoxon rank sum test was used.(D) The spearman's correlation of the RCD.GP score and the immune scores related with the immunotherapy response.(E-I) The K-M curves showing the difference of the prognosis between the high RCD.GP score subgroup and the low RCD.GP score subgroup in UC, SKCM, NSCLC and GBM.(J) The difference of the RCD.GP score between the NR patients and R patients.NR: not response to immunotherapy response.R: response to immunotherapy response.

Figure 8 .
Figure 8. Identification of SLC43A3 as novel and potential RCD-related oncology target.(A) The framework of screening the important features with machine learning methods.(B) The numbers of the DEGs, RCD related genes, and prognostic candidates.(C, D)The number of the important genes screened by different machine learning methods with different parameters in TCGA-Glioma cohort (C) and TCGA-GBM cohort (D).(E) The number of times each genes was screened out by the machine learning framework in TCGA-Glioma cohort and TCGA-GBM cohort.(F) The profile of the SLC43A3 in different grade glioma and normal brain cortex tissues in TCGA and GTEx datasets.Wilcoxon rank sum test was used for assessing the difference between the two subgroups."*", "**", "***", and "****" represented that the p value < 0.05, 0.01, 0.001, and 0.0001.(G) The difference of the SLC43A3 expression among the different transcriptomic subtypes in TCGA-GBM cohort.(H, I) The expression of the SLC43A3 in different glioma grades and different transcriptomic subtypes in CGGA cohort.J.The qRT-PCR result showing the SLC43A3 highly expressed in GBM compared with LGG and normal brain cortex tissues.NT: normal brain cortex tissue.The t test was utilized.(K) The qRT-PCR of the SLC43A3 relative expression in different glioma cell lines.The human astrocytes (HA1800) and three GBM cell lines including A172, LN229 and U87 were used for evaluating the expression of the score genes.The t test was utilized for confirming the difference.(L) The meta analysis provided a comprehensive HR of the RCD.GP score.M. Multivariate Cox regression analysis of the SLC43A3 was performed in TCGA-Glioma, CGGA-325 and CGGA-693 cohort. https://doi.org/10.1038/s41598-024-54643-3 Confirmation of valuable DEGs in glioma.The genes with the Wilcoxon rank sum test p value less than 0.05 were considered as DEGs.The intersection of DEGs between GBM and normal brain cortex, LGG and normal brain cortex, GBM and LGG was denoted as the critical DEGs in glioma. -