Abstract Background Osteosarcoma (OS) is an aggressive bone tumor with poor outcomes in advanced stages. Programmed cell death (PCD) plays a crucial role in tumor biology, but the prognostic value of integrating multiple PCD-related genes in OS remains underexplored. Methods Transcriptomic data from GEO and clinical data from TARGET were analyzed. Functional enrichment, N6-methyladenosine (m6A)-related genes analysis and immune checkpoint analyses were performed. Genes related to 18 types of PCD were intersected with OS-specific differentially expressed genes. Prognostic genes were identified by Kaplan-Meier and Cox regression. A LASSO model was developed to construct a survival prediction signature. Single-cell RNA-seq data were used to explore gene expression across cell types. Results There are 781 differentially expressed genes totally. M6A-related genes including ALKBH3, CBLL1 and immune checkpoints including HAVCR2, PDCD1 showed differential expression in OS. Twelve PCD-related genes were significantly associated with OS survival. A four-gene (FUCA1, GM2A, MAN2B1, COL13A1) prognostic model was established, showing strong predictive performance (3-year AUC = 0.801; 5-year AUC = 0.842). Lysosome-dependent cell death, apoptosis and anoikis emerged as key PCD pathways in OS. Single-cell analysis revealed COL13A1 expression in malignant cells while GM2A and FUCA1 were enriched in macrophages. Conclusion This study identifies a robust PCD-related prognostic model for OS and highlights key genes and pathways involved in tumor progression. The present findings determine potential biomarkers and therapeutic targets to improve OS prognosis and treatment. Supplementary Information The online version contains supplementary material available at 10.1007/s12672-025-02944-y. Keywords: Programmed cell death, OS, Overall survival, HAVCR2, COL13A1 Introduction Osteosarcoma (OS) is a malignant tumor that primarily affects adolescents and young adults. It is typically characterized by high aggressiveness and metastatic potential, posing a significant threat to patient survival [[24]1]. Although existing treatment methods (surgical resection, chemotherapy, radiotherapy, etc.) have improved the prognosis of patients to varying degrees, the treatment effect on critically ill patients is still unsatisfactory and often accompanied by serious complications [[25]2]. Therefore, there is an urgent need to further investigate the molecular mechanisms underlying tumorigenesis and development of OS disease and to identify more effective therapeutic targets and progression biomarkers. Programmed cell death (PCD), including apoptosis and autophagy, is a fundamental biological process essential for maintaining homeostasis in the body [[26]3]. In recent years, accumulating evidence has demonstrated that PCD is closely associated with the onset and progression of various diseases, particularly cancer [[27]4]. Numerous studies have shown that modulation of cell death pathways can significantly impact tumor cell survival and proliferation. The link between PCD and cancer has become increasingly evident, with research indicating that PCD not only influences the tumor microenvironment but also plays a role in mechanisms such as immune evasion [[28]5]. Experimental evidence suggests that PCD plays a critical role in the development of OS [[29]6]. For instance, BCL2, a molecule associated with apoptosis has been identified as a promising therapeutic target in OS and dysregulation of BCL-2 family members is associated with poor prognosis in OS patients [[30]7]. CHMP4C, a molecule related to pyroptosis, has also been implicated in OS progression [[31]8]. Li et al. reported that metformin induces cell cycle arrest, apoptosis and autophagy in human OS cells via the ROS/JNK signaling pathway, highlighting the role of autophagy in OS development [[32]9]. Another study found autophagy-related molecules such as TRIM68, PIKFYVE and DYNLL2 were associated with OS prognosis [[33]10]. Zhang et al. found that anoikis-related patterns are linked to the immune microenvironment of OS suggesting the involvement of anoikis in OS [[34]11]. In addition, another study identified several cuproptosis-related lncRNAs that are significantly correlated with OS prognosis [[35]12]. Collectively, these current findings highlight the significant role of PCD in the tumorigenesis and development of OS. However, most of them just focused on the role of a single form of PCD, overlooking the possibility of interactions among genes involved in different PCD pathways. Given this, the present research hypothesize that integrating multiple PCD-related genes into a comprehensive model may offer greater prognostic value. Therefore, based on previously identified PCD-related genes, this study aimed to construct a robust prognostic model capable of accurately predicting overall survival in OS patients. Methods and materials Data sources The transcriptomic data for Osteosarcoma (OS) was sourced from the Gene Expression Omnibus (GEO) database ([36]https://www.ncbi.nlm.nih.gov/geo/). The data of 3 normal tissue samples and 84 OS tissue samples is obtained from [37]GSE33382, and the data of 36 samples with survival information is obtained from [38]GSE21257. Additionally, the molecular data related to various forms of PCD were derived from the study by Zou et al. [[39]13]. In their research, twelve PCD patterns (apoptosis, necroptosis, pyroptosis, ferroptosis, cuproptosis, entotic cell death, netotic cell death, parthanatos, lysosome-dependent cell death, autophagy-dependent cell death, alkaliptosis and oxeiptosis) were analyzed for model construction, and a detailed list of molecules associated with each type of PCD was producted [[40]13]. The survival data of OS was obtained from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database ([41]https://www.cancer.gov/ccg/research/genome-sequencing/target). The TARGET program employs a comprehensive genomic approach to identify the molecular. Differential gene expression analysis (DEG) The dataset used in this study was [42]GSE33382 and the format of downloaded data was MINiML. DEG was conducted by Limma package (version 3.40.2) in R software. Adjusted P-values were calculated to correct for false-positive results. This study defined the adjusted P-value < 0.05 as a threshold to identify mRNAs with significantly differential expression in which log2(fold change) > 1 means upregulated genes while log2(fold change) < −1 means downregulated genes. Enrichment analysis To further explore the functions of potential targets, functional enrichment analysis was conducted. Gene Ontology (GO) analysis is a widely used method for annotating gene functions, specifically categorized into molecular function (MF), biological processes (BP) and cellular components (CC). Additionally, KEGG pathway enrichment analysis is a practical approach for examining gene functions and the associated high-level genomic information. The ClusterProfiler package in R software was used to analyze the GO terms and enriched KEGG pathways for potential mRNAs. DEG of m6A-related genes The source of m6A-related genes was obtained from the study conducted by Juan Xu et al., who performed molecular characterization and investigated the clinical significance of m6A regulators across 33 cancer types [[43]14]. Utilizing this dataset, this study explored differential expression of m6A-related genes between OS and normal tissues. The data of normal samples and OS samples was obtained from [44]GSE33382. The significance between these two groups was evaluated by Wilcoxon test. Immune checkpoint correlation analysis ITPRIPL1, SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3 and PDCD1LG2 are key genes associated with immune regulation. Analyzing the expression levels of these genes in both OS group and normal group can help this study to gain a more deep insight into the function of immune checkpoints [[45]15]. The data of normal samples and OS samples was obtained from [46]GSE33382. The significance of differential expression between these two groups was evaluated by Wilcoxon rank-sum test which enabling the exploration of potential associations between tumor development and these immune checkpoints. Intersection of upregulated genes and PCD-related genes in OS Upregulated genes were obtained from above DEG of [47]GSE33382. Genes involved in each type of PCD were obtained from previous study by Zou et al. [[48]13]. To get intersection of upregulated genes and PCD-related genes in OS, this study did Venn Diagram Analysis and Visualization by the online tool ([49]https://www.xiantaozi.com/). Overall survival analysis The prognostic data and gene expression of OS patients were downloaded from the TARGET database. According to the median value of each gene expression, all patients were divided into high- and low-expression groups. This study performed Kaplan-Meier (KM) survival analysis to assess survival difference between these two groups of intersecting genes, applying the log-rank test for statistical comparison. For the Kaplan-Meier curves, the p-value was calculated by the log-rank test while the hazard ratio (HR) and its corresponding 95% confidence interval (CI) were determined through univariate Cox regression analysis. When p-value < 0.05, genes were significantly associated with overall survival prognosis. Prognostic model Building These intersecting genes significantly associated with overall survival prognosis were used to construct the prognostic model. Lasso regression is a regularization method that achieves a more refined model by constructing a penalty function, which compresses some coefficients and sets others to zero. As a result, Lasso retains the advantages of subset shrinkage and serves as a biased estimation method for handling multicollinearity in data. It enables variable selection while performing parameter estimation, effectively addressing the multicollinearity problem in regression analysis. The model is optimized when lambda is minimized. Feature selection is performed by the Lasso regression algorithm, with a 10-fold cross-validation to construct the model [[50]16]. A larger AUC value and a smaller Log-rank p-value indicate better predictive performance. During model selection, this study extracts the AUC results from the survival ROC package, prioritizing genes with AUC value greater than 0.7. If gene model with an AUC value exceeding 0.7 can be found, then the model with the highest AUC value would be chosen [[51]17]. All analyses were conducted by the glmnet package in R software. The prognostic model was first constructed using TARGET data and then validated in [52]GSE21257 datasets. Distribution of prognostic model genes in cells TISCH ([53]http://tisch.comp-genomics.org) is a comprehensive database dedicated to single-cell RNA sequencing data [[54]18]. Through this platform, the present research explored the distribution of prognostic model-related genes across different cell types within OS tissues. All statistical analysis was conducted by R software, version v4.0.3. Results were considered statistically significant when the p-value was less than 0.05. Results Differential gene expression Through differential analysis of Osteosarcoma (OS) expression profile data from GEO, this study identified 399 upregulated genes and 382 downregulated genes. Figure [55]1A presents a volcano plot illustrating the distribution of differentially expressed genes. Figure [56]1B shows the top 50 of significantly upregulated genes and top 50 of downregulated genes. The full list of gene names is provided in Supplementary Table 1. Enrichment analysis of upregulated genes revealed that the most significant biological process (BP) was antigen processing, the most significant molecular function (MF) was extracellular matrix structural constituent and the most significant cellular component (CC) was collagen-containing extracellular matrix (Fig. [57]1C). KEGG pathway enrichment identified several pathways, including Human T-cell leukemia virus 1, phagosome and cell adhesion molecules (Fig. [58]1D). Enrichment analysis of downregulated genes found that the most enriched biological processes were related to reproductive structure development, the most significant molecular function was receptor-ligand activity and the most enriched cellular component was collagen-containing extracellular matrix (Fig. [59]1E). KEGG enrichment revealed pathways such as transcriptional misregulation in cancer, cellular senescence and cytokine-cytokine receptor interaction (Fig. [60]1F). Fig. 1. [61]Fig. 1 [62]Open in a new tab Differential expression analysis in OS. (A) shows a volcano plot illustrating genes with differential expression between OS group and normal group. Upregulated genes locate upper right area with log2(fold change) value > = 1 and downregulated genes locate upper left area with log2(fold change) value <=−1. (B) presents the top 50 upregulated genes and top 50 downregulated genes which highlighting the most significant expression changes. (C) and (D) mean GO and KEGG enrichment analysis of upregulated genes. (E) and (F) mean GO and KEGG enrichment analysis of downregulated genes DEG analysis of m6A-related genes This study found that the expression of some m6A-related genes in OS group significantly higher than that in normal group (also defined as control group). Specifically, ALKBH3 (P = 0.006), CBLL1 (P = 0.014), IGF2BP2 (P = 0.014), METTL3 (P = 0.012), RBMX (P = 0.035), YTHDC2 (P = 0.004) and YTHDF2 (P = 0.049) showed significant difference. However, other m6A-related genes did not remarkably exhibit difference (Fig. [63]2A-B). Fig. 2. [64]Fig. 2 [65]Open in a new tab Expression profiles of m6A-related genes between tumor tissues and normal tissues. The x-axis represents different m6A-related genes while the y-axis shows their expression levels (A-B). Colors indicate different sample groups (red means tumor while blue means normal). The number of asterisks denotes statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ns = not significant. Statistical difference between two groups was evaluated by Wilcoxon test Immune checkpoint correlation Several genes belong to immune checkpoints exhibited differential expression between the OS group and normal group (also defined as control group), including HAVCR2 (P = 0.004), ITPRIPL1 (P = 0.004), LAG3 (P = 0.006), PDCD1 (P = 0.006), PDCD1LG2 (P = 0.035) and TIGIT (P = 0.037). No significant difference was observed for the rest immune checkpoints (Fig. [66]3). Fig. 3. [67]Fig. 3 [68]Open in a new tab Expression profiles of immune checkpoint genes between tumor tissues and normal tissues. The x-axis represents different sample groups (red means tumor while blue means normal), and the y-axis indicates the expression levels of immune checkpoint genes. Asterisks indicate statistical significance: *p < 0.05, **p < 0.01, ***p < 0.001, ns = not significant. Statistical difference between two groups was assessed by Wilcoxon rank-sum test Intersection of upregulated genes and PCD-related genes in OS By intersecting the upregulated genes in OS and genes involved in Programmed Cell Death (PCD), this study found that most intersecting genes were related to Lysosome-dependent cell death (with a total number of 21). This was followed by Apoptosis (15 intersecting genes) and Anoikis (13 intersecting genes). Fewer intersections were observed in other PCD types. Details of these genes are listed in Table [69]1 while the visualization is shown as Fig. [70]4A. Figure [71]4B presents the results of a subsequent enrichment analysis of these intersected genes. This analysis revealed that these genes not only being associated with their respective PCD types but also linked to pathways such as Lipid and Atherosclerosis, positive regulation of response to external stimulus and myeloid leukocyte activation. Table 1. Intersection of upregulated genes and PCD-related genes in OS Programmed Cell Death High-Expression Genes in OS Apoptosis BTK, CD14, CD74, CRIP1, CTSC, E2F2, FGFR3, HMGB2, MLLT11, MMP9, MSX1, PARP1, TNFSF10, TXNDC12, TYROBP Pyroptosis GSDMD IL18 Ferroptosis CYBB MT1G Autophagy LZTS1 RRAGD Necroptosis TNFSF10 PARP1 Cuproptosis - Parthanatos - Entotic cell death CYBB Netotic cell death - Lysosome-dependent cell death ACP5, AP1S2, BTK, CD300A, CD68, CD84, CTSC, CTSH, CTSK, CTSZ, FUCA1, GM2A, GNPTAB, GUSB, LAPTM4B, LAPTM5, LGMN, LYN, MAN2B1, PPT1, VAMP8 Alkaliptosis - Oxeiptosis - NETosis GSDMD, PLA2G7, MMP9 Immunogenic cell death CALR Anoikis CSPG4, CALR, CDKN3, MMP13, MMP9, CXCR4, TNFSF10, S100A4, PARP1, HAVCR2, CD36, COL13A1, MMP11 Paraptosis CDKN3 Methuosis - Entosis TNFSF10, KIF2C [72]Open in a new tab A dash (-) is used to indicate the absence of intersecting genes Fig. 4. [73]Fig. 4 [74]Open in a new tab Intersection of upregulated genes and PCD-related genes in OS. (A) shows a Venn diagram illustrating the intersection of upregulated genes and PCD-related genes in OS. (B) presents the results of the enrichment analysis of these intersecting genes Overall survival analysis in OS After these intersecting genes were determined, this study did overall survival analysis for them using OS dataset downloaded from the TARGET database. All information of the results is provided in Supplementary Table 2. This study found that 12 genes were significantly associated with overall survival prognosis in OS with p-value < 0.05 (Table [75]2). These genes are CD14 (HR = 0.485, 95% CI: 0.251 to 0.939, P = 0.032), FGFR3 (HR = 2.049, 95% CI: 1.063 to 3.947, P = 0.032), TNFSF10 (HR = 0.453, 95% CI: 0.230 to 0.892, P = 0.022), TYROBP (HR = 0.383, 95% CI: 0.195 to 0.753, P = 0.005), CTSK (HR = 0.486, 95% CI: 0.251 to 0.941, P = 0.032), FUCA1 (HR = 0.367, 95% CI: 0.184 to 0.730, P = 0.004), GM2A (HR = 0.486, 95% CI: 0.250 to 0.943, P = 0.033), LAPTM5 (HR = 0.458, 95% CI: 0.236 to 0.886, P = 0.020), MAN2B1 (HR = 0.398, 95% CI: 0.205 to 0.772, P = 0.006), VAMP8 (HR = 0.352, 95% CI: 0.177 to 0.700, P = 0.003), HAVCR2 (HR = 0.359, 95% CI: 0.183 to 0.706, P = 0.003), COL13A1 (HR = 3.003, 95% CI: 1.502 to 6.002, P = 0.002). The survival analysis results of these 12 genes are visualized as Fig. [76]5A-L while HR values and confidence intervals are displayed in a forest plot (Fig. [77]6). Table 2. The list of 12 prognosis-related genes Genes P value HR Low 95%CI High 95%CI COL13A1 0.002 3.003 1.502 6.002 VAMP8 0.003 0.352 0.177 0.700 HAVCR2 0.003 0.359 0.183 0.706 FUCA1 0.004 0.367 0.184 0.730 TYROBP 0.005 0.383 0.195 0.753 MAN2B1 0.006 0.398 0.205 0.772 LAPTM5 0.020 0.458 0.236 0.886 TNFSF10 0.022 0.453 0.230 0.892 CD14 0.032 0.485 0.251 0.939 FGFR3 0.032 2.049 1.063 3.947 CTSK 0.032 0.486 0.251 0.941 GM2A 0.033 0.486 0.250 0.943 [78]Open in a new tab Fig. 5. [79]Fig. 5 [80]Fig. 5 [81]Open in a new tab Kaplan–Meier survival curves of intersecting genes. This figure only displays the Kaplan–Meier (KM) survival curves of 12 genes significantly associated with overall survival prognosis with p-value < 0.05. A risk table is included below each plot, indicating the number of patients remaining at risk at specific time points. These genes presented are (A) CD14, (B) TYROBP, (C) MAN2B1, (D) VAMP8, (E) LAPTM5, (F) FUCA1, (G) GM2A, (H) CTSK, (I) TNFSF10, (J) HAVCR2, (K) COL13A1 and (L) FGFR3 Fig. 6. [82]Fig. 6 [83]Open in a new tab Forest plot of 12 intersecting genes associated with overall survival in OS. This figure shows a forest plot illustrating the hazard ratios (HRs), 95% confidence intervals (CIs) and p-values of the genes positively associated with overall survival in OS Overall survival prognostic model 12 prognosis-related genes were used to construct Lasso regression model in this study (Fig. [84]7A-B). The prediction model included a total of 95 OS patients, with an average 2.7 survival years. Finally, four most significant features (FUCA1, GM2A, MAN2B1, COL13A1) were identified as predictors for the prognostic risk model (Fig. [85]7C-E). Fig. 7. [86]Fig. 7 [87]Open in a new tab Construction and validation of the prognostic risk model based on selected gene features. (A) means a Lasso analysis used to construct a four-gene model. (B) is the gene screening of Lasso regression. (C) shows the distribution of risk score and overall survival status in the training set. (D) displays the Kaplan-Meier plot of risk score and overall survival in the training set. (E) represents the ROC curve of risk score in the training cohort. (F) shows the distribution of risk score and overall survival status in the validation set. (G) means the Kaplan-Meier plot of risk score and overall survival in the validation set. (H) is a ROC curve of risk score in the validation cohort The risk model was evaluated by using Kaplan-Meier survival curves. The Log-rank test was applied for group comparison. When conditions, including HR = 5.66, 95% CI: 2.58 and 12.417, P = 1.53E-05, meet together, the model is a risk model. The final prognostic model was presented under the equation: Riskscore = (−0.1455)*FUCA1+(−0.5933)*GM2A+(−0.1465)*MAN2B1+(0.4545)*COL13A1. Furthermore, the ROC analysis at different time points showed the following results: one-year overall survival prediction AUC = 0.668, 95% CI: 0.494 to 0.842; 3-year overall survival prediction AUC = 0.801, 95% CI: 0.719 to 0.883; 5-year overall survival prediction AUC = 0.842, 95% CI: 0.762 to 0.922. Both the 3-year and 5-year overall survival predictions had an AUC value greater than 0.8, indicating that the model performs exceptionally well in these two predictions. Figure [88]7F-H illustrates the prognostic performance of above four molecules by overall survival analysis. And the data source was downloaded from [89]GSE21257 dataset. Notably, the model constructed with these four molecules demonstrated a significant association with poor prognosis (HR = 3.57, 95% CI: 1.14–11.14, P = 0.028). In addition, the model yielded an AUC value of 0.72 for five-year survival prediction. Analysis of distribution of prognostic model genes in OS To further assess the roles of above four molecules in OS, this study conducted the single-cell analysis to identify which cell types predominantly express these molecules. Figure [90]8A shows the cell clusters where eight cell types (Fibroblasts, Mono/Macro, Osteoblasts, CD4Tconv, CD8Tex, Plasmocytes, Endothelial and Malignant) are annotated. The result shows that COL13A1 is primarily expressed in Malignant cells (Fig. [91]8B), GM2A and FUCA1 are predominantly expressed in Mono/Macro cells (Fig. [92]8C and D) and MAN2B1 is mainly expressed in Osteoblasts (Fig. [93]8E). Expression details can be found in Supplementary Table 3. Fig. 8. [94]Fig. 8 [95]Open in a new tab Single-cell analysis of gene expression in OS. (A) shows the overall cellular distribution based on single-cell RNA sequencing data. (B) represents COL13A1 expression distribution. (C) represents GM2A expression distribution. (D) represents FUCA1 expression distribution. (E) represents MAN2B1 expression distribution Discussion In this study, intersection of upregulated genes and PCD-related genes in OS was identified. This approach led to the discovery of 12 prognosis-related molecules, including CD14, FGFR3, TNFSF10, TYROBP, CTSK, FUCA1, GM2A, LAPTM5, MAN2B1, VAMP8, HAVCR2 and COL13A1. A prognostic model for overall survival in OS was constructed by the LASSO regression model based on the four most significant molecular features (FUCA1, GM2A, MAN2B1 and COL13A1), which resulting in the following Risk Score formula: Riskscore = (− 0.1455) × FUCA1 + (− 0.5933) × GM2A + (− 0.1465) × MAN2B1 + (0.4545) × COL13A1. Additionally, this study identified several immune checkpoint molecules (HAVCR2, ITPRIPL1, LAG3, PDCD1, PDCD1LG2 and TIGIT) as well as m6A–related genes (ALKBH3, CBLL1, IGF2BP2, METTL3, RBMX, YTHDC2 and YTHDF2) that may affect the tumorigenesis and progression of OS. These molecules represent potential therapeutic targets and prognostic biomarkers for OS. Among these 12 molecules associated with OS prognosis identified in the present study, CTSK, FUCA1, GM2A, LAPTM5, MAN2B1 and VAMP8 are related to lysosome-dependent cell death (LCD). LCD is a form of cell death triggered by increased lysosomal membrane permeability, which leads to the release of hydrolytic enzymes such as cathepsins into the cytoplasm, initiating a cascade of events that ultimately result in cell death [[96]19]. Li et al. demonstrated that Verteporfin slows the progression of spontaneous OS in Trp53/Rb1-deficient, CTSK-expressing cells by inhibiting the Hippo pathway, highlighting the potential importance of CTSK-expressing cells in OS [[97]20]. Wang et al. found that FUCA1 may play an important role in tumor-associated macrophages in OS tissues and also involved in immune response, signal transduction, TNF regulation and phagocytosis [[98]21]. That is consistent with the present findings that FUCA1 is positively associated with favorable prognosis. Additionally, Li et al. identified LAPTM5 as a hub gene in OS [[99]22]. The function of VAMP8 has also been validated by Yang et al. Their study reported that VAMP8 suppresses metastasis through the DDX5/β-catenin signaling pathway, which indicating its potential as a therapeutic target [[100]23]. Although the roles of GM2A and MAN2B1 in OS have not been extensively studied, this study revealed that they are highly expressed in tumor tissues and closely associated with prognosis with high expression of GM2A in monocyte-macrophage populations. Collectively, these findings suggest that lysosome-dependent cell death may play a crucial role in the development and progression of OS and its related molecules would be potential therapeutic targets and prognostic biomarkers. This study shows another form of programmed cell death (PCD), apoptosis, plays a significant role in OS prognosis with four associated molecules (CD14, FGFR3, TNFSF10 and TYROBP). Apoptosis is the principal pathway of regulated cell death [[101]24] and one of the most knowledgeable PCD mechanisms in OS [[102]25]. Among these molecules, CD14 has been widely reported in previous studies of OS. For instance, the importance of CD14⁺&CD16⁺ macrophages in OS progression was determined by Cillo et al. [[103]26]. Another study shows that CD14 is inherently linked to immune responses [[104]27]. Joung CI, et al. shows that CD14⁺ monocytes and soluble CD14 in synovial fluid are associated with osteoarthritis progression [[105]28]. However, there is no study explores the direct evidence connecting CD14 to OS right now. Studies about the function of apoptosis-related genes (FGFR3, TNFSF10 and TYROBP) in OS remain limited. Notably, TNFSF10, a member of the TNF superfamily, has demonstrated roles in various tumors. Han et al. have reported its involvement in antiviral immune regulation in triple-negative breast cancer [[106]29] while Charbonneau et al. have discussed its association with the NF-κB pathway in ovarian cancer [[107]30]. These findings indicate potential mechanisms by which apoptosis-related genes may influence OS biology. Further investigation is still needed to elucidate their exact functions in this context. Anoikis has recently emerged as a highly promising research phenotype in OS [[108]31]. The results of this study shows two Anoikis-related genes (HAVCR2 and COL13A1) are significantly associated with OS prognosis. HAVCR2 locates on chromosome 5q33.2 and encodes the T-cell immunoglobulin mucin-3 (TIM-3) protein, which has been proposed as a potential immune checkpoint target in various tumors [[109]32]. Yang et al. reported a higher expression of HAVCR2 in OS risk groups than in low-risk patients [[110]33]. Consistent with this, this study reveals that high HAVCR2 expression is associated with better overall survival, supporting its role as a favorable prognostic marker in OS. Notably, the findings of this study suggests that COL13A1 expression is strongly correlated with OS malignancy, as it is primarily expressed in highly malignant tumor cells. Several previous studies also supported its potential role as a therapeutic target and prognostic biomarker in OS [[111]34–[112]36]. However, the underlying molecular mechanism of COL13A1 in OS progression remains unclear. Although our study primarily focused on the types of PCD significantly associated with OS prognosis, this does not imply that other forms of PCD lack research value in this context. For instance, ferroptosis [[113]37] and necroptosis [[114]38] have been extensively investigated and the corresponding results showed a broader influence of PCD-related molucules on the tumorigenesis and development of OS. The present study determined the potential biomarkers and therapeutic targets in OS by analyzing the integration of PCD-related genes and highly expressed genes in OS along with constructing a prognostic mode. Nonetheless, there are certain limitations. Given that OS is a rare malignancy with relatively low incidence, it is difficult to collect an independent clinical cohort to further validate the comprehensive prognostic model. Furthermore, this study lacks in-depth exploration of the underlying molecular mechanisms. To sum up, the present study firstly determined the relationship beteween PCD-related genes and OS, and the findings would contribute the future research, including the exact mechanism of immune regulation, RNA methylation and other cellular processes in OS. Electronic supplementary material [115]Supplementary Material 1^ (1.6MB, xlsx) Acknowledgements