Exploring RCD pathway heterogeneity in glioma
To characterize the heterogeneity of cellular death pathways in gliomas, we utilized GSVA to estimate the enrichment scores for RCD pathways in each sample within the TCGA cohort. Through clustering analysis, we determined that the optimal number of clusters was 4 (Fig. 2A and S1A-G). Subsequently, based on the enrichment scores of the RCD pathways, we classified 667 glioma tumors into 4 distinct heterogeneous subtypes (Fig. 2B). Cluster 1 exhibited relatively specific high expression of alkaliptosis and netotic cell death, while Cluster 2 showed relatively specific high expression of ferroptosis and entotic cell death. Clusters 3 and 4 displayed distinct expression characteristics in Mpt-driven necrosis, disulfidepotosis, parthanatos, autophagy, and cuproptosis. Moreover, compared to the other three clusters, Cluster 1 demonstrated the poorest clinical prognosis, with the lowest overall survival and disease-specific survival rates (Fig. 2C,D). Subsequently, we conducted functional enrichment analysis for the four clusters. In the HALLMARKS database, patients in Cluster 1 were found to be more enriched in Mtorc1 signaling, reactive oxygen species pathway, hypoxia, EMT, and DNA repair. In the Reactome database, cases in Cluster 1 showed greater enrichment in events such as Cyclin A B1 B2 Associated Events During G2 M Transition, Unwinding of DNA, and release of apoptotic factors from the mitochondria, among others (Fig. 2E, S1H).
Clinical and immune profiling of four clusters
In the subset of grade 2 gliomas with relatively favorable prognosis, Cluster 3 and Cluster 4 show significant representation, with Cluster 3 being predominant. Notably, in Grade 4 gliomas, Cluster 1, associated with the poorest prognosis, holds the largest share, reaching 80%. Among other clinical features, Cluster 1 consistently dominates in clinical phenotypes associated with poorer prognosis, such as glioblastoma, IDH wild-type, 1p/19q non-codeletion, MGMT unmethylated, and others (Fig. 3A,B and S2A-F).
Subsequently, we utilized the CIBERSORT algorithm to calculate the proportions of immune cells in each cluster. The results revealed variations in several cell populations, including macrophage M2, monocytes, T cells CD4 memory resting, and plasma cells (Fig. 3C and S2G). Furthermore, utilizing four distinct algorithms, we evaluated the differences in cellular composition within the microenvironment of each cluster. In the EPIC algorithm, Cluster 1 exhibited the highest relative expression of Cancer-Associated Fibroblasts (CAFs) and endothelial cells. According to the xCELL algorithm, Cluster 1 showed elevated expression of astrocytes and CD4 + memory T cells. The MCPcounter algorithm indicated a relatively higher presence of endothelial cells and fibroblasts in Cluster 1. In the QUANTISEQ algorithm, Cluster 1 demonstrated a higher relative presence of macrophages and T cells (Fig. 3D).
Additionally, we utilized the ESTIMATE algorithm to compute the tumor immune-related scores. Compared to the clusters with better prognosis, namely Cluster 3 and Cluster 4, Cluster 1 exhibits higher ESTIMATE Score, Immune Score, and Stromal Score. Moreover, Cluster 1 has the lowest tumor purity (Fig. 3E,F).
Genetic landscape analysis of four clusters
We proceeded to elucidate the genetic variation characteristics of the tumors. Initially, we computed the MATH scores for the four clusters, revealing that Cluster 1 exhibits the lowest MATH score among them (Fig. 4A). Conversely, Cluster 1 demonstrates the highest TMB score (Fig. 4B). Subsequent scrutiny of the prognostic implications of these clusters indicated that the ClusterA-TMB high group is associated with the most unfavorable prognosis, while the ClusterNonA-TMBlow group exhibits the most favorable survival outcomes (Fig. 4C). Within the glioma patient cohort from TCGA, approximately 91.19% (611/670) displayed genetic mutations. Notably, IDH1 mutations were most prevalent (60%), with Cluster 1 predominantly devoid of IDH1 mutations. Additional mutations encompassed TP53 (44%), ATRX (29%), CIC (16%), and others (Fig. 4D). Lastly, we provided an account of gene amplifications or deletions across chromosomal positions. Cluster 1 exhibited heightened mutation rates across all chromosomes compared to the other three clusters, with amplifications predominantly localized to chr4, chr7, and chr12, while deletions were principally concentrated on chr1, chr4, and chr9 (Fig. 4E-H).
Exploring prognostic signatures across four clusters
To further explore the association of genes with prognosis across different glioma clusters, we initially conducted a comprehensive analysis of gene expression differences between Cluster 1, indicative of the poorest prognosis, and the other three clusters (Fig. 5A-C). The Principal Component Analysis (PCA) plot demonstrated distinct and independent patient distribution landscapes across the four clusters (Fig. 5D). Subsequently, we identified a set of differentially expressed genes, comprising 234 commonly upregulated genes and 75 commonly downregulated genes (Fig. 5E,F). Leveraging these 309 differential genes, we employed a range of sophisticated machine learning algorithms to develop a robust prognostic model for gliomas, utilizing the TCGA dataset for training and CGGA325 and CGGA693 datasets for validation. Among the diverse algorithms evaluated, the CoxBoost + SuperPC combination demonstrated superior performance, leading us to select 12 genes for the construction of the RCD-related prognostic signature (Fig. 5G,H).
Comprehensive prognostic evaluation of RCD-related signature
We conducted a comprehensive evaluation of the prognostic significance of the RCD-related model built upon the 12 identified genes. In the TCGA training set, patients classified as RRScore-high group exhibited markedly lower average survival rates compared to their low-risk counterparts, with all areas under the curve (AUC) values surpassing 0.85 (Fig. 6A). The CGGA325 and CGGA693 validation sets further confirmed the shorter overall survival in the RRScore-high group, corroborated by AUC values indicating robust prognostic assessment (Fig. 6B,C). Expanding the validation scope, we incorporated the GSE55918 and GLASS datasets as additional verification sets. In these datasets, patients in the RRScore-high group consistently displayed significantly reduced survival times relative to those in the RRScore-low group, as elucidated by Kaplan-Meier curves highlighting distinct prognostic disparities (Fig. 6F,G). Furthermore, we performed univariate and multivariate risk assessments for the prognostic signature. The outcomes identified age, grade, IDH status, 1p/19q codeletion, and riskscore as independent prognostic factors (Fig. 6D,E). Integrating the riskscore with all variables, we developed a nomogram to quantitatively depict the prognostic index associated with glioma patients (Fig. 6H). The study findings underscored a c-index of 0.871, underscoring its robust capacity for prognostic prediction (Fig. 6I).
Exploring the clinical implications of RCD-related signature in glioma
The prognostically relevant genes exhibit positive correlations with alkaliptosis, apoptosis, ferroptosis, and various other forms of cell death, while showing negative correlations with MPT-driven necrosis and parthanatos. Notably, PTPRT displays a contrasting trend to the associations (Fig. 7A). We further explored the potential clinical implications of the signature in terms of therapeutic response. Through drug sensitivity correlation analysis, we identified 13 compounds as potential therapeutic agents for gliomas. Compounds such as Dasatinib and Gemcitabine demonstrated high IC50 values in the RRScore-low group, whereas Daporinad and Lapatinib exhibited elevated IC50 values in the RRScore-high group (Fig. 7B). The relationships among glioma prognostic genes, RCD, and predicted drug pathways are depicted in Fig. 7C.
Additionally, the RRScore-high group demonstrated relatively lower Microsatellite Instability (MSI) scores but higher Exclusion scores, IFNG scores, and TIDE scores (Fig. 7D). Lastly, we conducted an analysis of the correlation with immune modulators, encompassing co-stimulators, co-inhibitors, ligands, cell adhesion molecules, receptors, antigen presentation components, and others. Most of the immunomodulators in the high RRScore group exhibited relatively high expression levels, with only a few exceptions such as VTCN1 and IL12A, which showed no significant correlation (Fig. 7E).
Single cell analysis of RCD-related signature
We conducted an in-depth analysis at the single-cell level to explore programmed cell death and its potential impact on glioma tumor progression. Cells with high riskscore predominantly clustered within malignant cells, with a subset also distributed among immune cells in CPTAC dataset (Fig. 8A). Employing published datasets consisting of both LGG and GBM samples [26], we observed a lower average communication number in the high-risk group, indicating potential differences in cell communication (Fig. 8B) upon stratifying patients based on their scores. Notably, interaction strength was relatively higher in malignant cells of the high-risk group, while myeloid cells exhibited stronger interactions in the low-risk group (Fig. 8C). Subsequent analysis revealed signaling alterations between different cell types, including pathway changes in PTN, MK, EGF, MIF, and PSAP in various malignant cell subtypes, along with alterations in the CCL pathway in myeloid cells (Fig. 8D). Comparative analysis of relative information flow between high and low-risk groups highlighted predominant alterations in pathways like ncWNT, CHEMERIN, CD137, and VEGF in the high-risk group, whereas changes in CX3C, CSF3, EDN, and PROS pathways were notable in the low-risk group (Fig. 8E). Furthermore, we compared signaling pathway alterations within different cell types in each group. In the high-risk group, pathways such as CCL3, CCL5, IGF, and TGFB were significantly upregulated (Fig. 8F). Lastly, we identified potential target genes for myeloid cells, including SPP1, HMOX1, SOCS3, CXCL2, and NUPR1, which could serve as targets for therapeutic interventions (Fig. 8G).
Spatial characterization of RCD-related signature
We explored the spatial characteristics of the RCD-related signature within the glioma microenvironment. We computed RRScore for various cell types and observed that diff.-like and prolif.stem-like cancer cells exhibited the highest RRScore across all cell types, while stem-like malignant cells and oligodendrocytes displayed the lowest RRScore (Fig. 9A). Spatial cell type deconvolution analysis revealed distinct spatial distributions of transcriptome-defined cell types in glioma. Cells with high RRS in the diff.-like cancer subtype displayed unique spatial niches compared to those with low RRS, forming a hypoxic, EMT-high tumor core. This core, reminiscent of a histologically defined pseudopalisading zone [27], exhibited regulatory cell death patterns. Notably, myeloid cells were excluded from the tumor core in diff.-like RRS-high cells, despite having a similar spatial distribution overall (Fig. 9B-D).
malignant phenotype experiments of tumor cells
We initially performed TIMP1 staining on pathological specimens from glioma patients, revealing a relative increase in TIMP1 protein expression in high-grade gliomas (Fig. 10A). Subsequently, we conducted cellular experiments to validate these findings. Upon downregulating TIMP1 expression (Fig. 10B), a notable decrease in cell proliferation capacity was observed (Fig. 10C,D). Furthermore, results from colony formation assays indicated a reduction in both the size and number of cell colonies following TIMP1 knockdown (Fig. 10E). Similarly, findings from Transwell and scratch assays demonstrated that TIMP1 promotes the migration and invasion capabilities of glioma cells (Fig. 10F).