3.1 Establishment of the inflammation-related risk scoring model

We conducted a univariate Cox regression analysis based on the inflammation gene set, and Fig. 1 illustrates our streamlined workflow for deriving the 11‑gene prognostic signature in a single integrated view: Fig. 1A presents the univariate Cox analysis of 13 inflammation‑related candidates (all p < 0.01), with hazard ratios indicating both protective (e.g., ITGAM, HR 0.862) and risk (e.g., TERT, HR 1.291) factors; Fig. 1B simplifies these candidates’ co‑expression network via a chord diagram, highlighting strong positive links (e.g., SAMHD1–ABCB4) and inverse relationships (e.g., PPARG–CD36); Fig. 1C tracks each gene’s LASSO coefficient across log(λ), emphasizing in bold the 11 genes that retain nonzero weights at the optimal λ_min (log(λ)=–3.20); Fig. 1D plots the 10‑fold cross‑validated partial‑likelihood deviance, with dashed lines marking λ_min and the more parsimonious λ_1se, thereby transparently defining the cutoff that yields our final, robust inflammation‑driven signature; and Fig. 1E presents GO and KEGG enrichment of those 11 genes across Biological Process (e.g., regulation of endothelial cell development, negative regulation of mRNA‑mediated gene silencing), Cellular Component (e.g., endocytic vesicle, platelet alpha granule), Molecular Function (e.g., scavenger receptor activity, transcription coregulator binding), and KEGG pathways (e.g., adipocytokine signaling, PPAR signaling, lipid and atherosclerosis), with bubble size denoting gene count and color indicating adjusted p‑value. Multivariate Cox regression analysis incorporating clinical traits like age, metastasis status, and gender indicated that metastasis, PPARG, TERT, and VEGFA were closely linked to osteosarcoma prognosis (Fig. 2A). A nomogram incorporating these four factors was constructed to predict the 1, 3, and 5-year prognosis of osteosarcoma patients (Fig. 2B). Calibration curves for 1, 3, and 5 years (Figs. 2C-E) demonstrated a close match between observed and predicted survival probabilities, indicating the nomogram’s strong prognostic value.

Fig. 1

Inflammation-Related Genes and Their Prognostic Significance in Osteosarcoma. A: Forest plot illustrating the hazard ratios of 13 inflammation-related genes identified through univariate Cox regression analysis, highlighting their association with osteosarcoma prognosis. Genes such as CD163, ITGAM, PPARG, and WAS are depicted as favorable prognostic factors, while CBS, TERT, ABCB4, CD36 are shown as adverse factors. B: Chord diagram representing significant positive correlations among selected genes (WAS, CD163, ITGAM, SAMHD1), indicating their potential synergistic roles in osteosarcoma prognosis. C-D: Lasso regression analysis outcomes, showcasing the selection process of 11 critical genes from the inflammation-related gene set. These genes’ expression patterns were further analyzed for their prognostic implications in osteosarcoma

Fig. 2

Construction and Validation of a Prognostic Nomogram for Osteosarcoma. A: Multivariate Cox regression analysis forest plot indicating the prognostic importance of clinical traits and selected genes (PPARG, TERT, VEGFA) alongside metastasis status in osteosarcoma patients. B: Prognostic nomogram incorporating metastasis status, PPARG, TERT, and VEGFA expression levels for predicting the 1-, 3-, and 5-year survival probabilities of osteosarcoma patients. C-E: Calibration curves for the nomogram predictions at 1, 3, and 5 years, demonstrating a close alignment between the predicted and observed survival probabilities, affirming the nomogram’s predictive accuracy

3.2 Single-cell sequencing analysis

Based on the GSE1624554 cohort, single-cell sequencing analysis was conducted, with results in Fig. 3A and B showing the smallest standard deviation at PC = 7 and a chosen clustering tree resolution of 1.2. UMAP plots in Figs. 3C-D illustrated the distribution of cell types such as myeloid cells, macrophages, stromal cells, NK cells, megakaryoblast, fibroblast cells, and endothelial cells. The heatmap in Fig. 3E displayed the enrichment levels of inflammation-related genes in these seven cell clusters, with genes like BCL2A1, S100A9, CCL3L1 significantly enriched in myeloid cells, and ACP5, SIGLEC15, ATP6V0D2 in macrophages, SFRP2, MEG3, CXCL14 in stromal cells. Figures 4A-B showed the enrichment levels of 11 inflammation-related genes associated with osteosarcoma across the seven cell types, with significant enrichment of CD163 in macrophages and myeloid cells, SAMHD1 in macrophages and myeloid cells, and TNFRSF1A in endothelial and fibroblast cells.

Fig. 3

Single-cell RNA Sequencing Analysis of Osteosarcoma Tissue. A-B: Analysis of single-cell RNA sequencing data from the GSE1624554 cohort, highlighting the optimal principal component (PC = 7) for minimal standard deviation and the chosen resolution (1.2) for cluster analysis. C-D: UMAP plots depicting the distribution of various cell types within the osteosarcoma microenvironment, including myeloid cells, macrophages, stromal cells, NK cells, megakaryoblasts, fibroblast cells, and endothelial cells, illustrating the cellular heterogeneity. E: Heatmap showing the enrichment levels of inflammation-related genes across seven identified cell clusters, revealing significant gene expression patterns in myeloid cells, macrophages, and stromal cells, which may contribute to the tumor microenvironment’s complexity

Fig. 4

Distribution of Inflammation-Related Genes Across Identified Cell Clusters. A-B: Analysis of the expression levels of 11 key inflammation-related genes within the seven cell types identified in the osteosarcoma microenvironment. The figures highlight the significant enrichment of CD163 and SAMHD1 in macrophages and myeloid cells, and TNFRSF1A in endothelial and fibroblast cells, underscoring their potential roles in modulating immune responses within the tumor microenvironment

3.3 Prognostic and functional enrichment analysis between high and low-risk groups

Osteosarcoma patients were divided into high and low-risk groups based on their risk scores. The Kaplan-Meier (KM) curve in Fig. 5A indicated significant prognostic differences between the groups, with the low-risk group showing significantly better prognosis than the high-risk group. The ROC curve in Fig. 5B demonstrated that the risk scoring model had AUC values of 0.808, 0.883, and 0.879 for predicting 1, 3, and 5-year survival, respectively. The survival status scatter plot in Fig. 5C revealed that patients in the low-risk group tended to have a survival status, while those in the high-risk group were more likely to be deceased. Decision Curve Analysis (DCA) in Figs. 5D-F showed that the risk scoring model had a higher net benefit for predicting 1, 3, and 5-year overall survival (OS) when the threshold probability exceeded 0.1, outperforming other factors like metastasis, PPARG, TERT, and VEGFA. Furthermore, differentially expressed genes (DEGs) between high and low-risk groups were identified and subjected to GO and KEGG enrichment analyses. GO enrichment analysis in Fig. 6A revealed that DEGs were mainly enriched in pathways like positive regulation of cytokine production, external side of plasma membrane, and leukocyte-mediated immunity. KEGG enrichment analysis showed DEGs were predominantly enriched in pathways such as B cell receptor signaling pathway, osteoclast differentiation, and primary immunodeficiency (Fig. 6B).

Fig. 5

Prognostic Analysis and Risk Assessment in Osteosarcoma Patients. A: Kaplan-Meier survival curves comparing overall survival between high and low-risk groups, demonstrating significant prognostic differences with better outcomes in the low-risk group. B: Receiver operating characteristic (ROC) curves for the risk scoring model, showing area under the curve (AUC) values for 1, 3, and 5-year survival predictions, indicating high predictive accuracy. C: Scatter plot of survival status against risk score, highlighting the distribution of survival outcomes within high and low-risk groups, with deceased patients predominating in the high-risk group. D-F: Decision curve analysis (DCA) for 1, 3, and 5-year overall survival predictions, illustrating the net benefit of the risk scoring model across different threshold probabilities and its superior predictive power compared to other clinical factors

Fig. 6

Functional Enrichment Analysis of Differentially Expressed Genes Between high and low-Risk Osteosarcoma Groups. A: Gene Ontology (GO) enrichment analysis of DEGs showcasing significant enrichment in biological processes related to the positive regulation of cytokine production, the external side of the plasma membrane, and leukocyte-mediated immunity, highlighting the immunological underpinnings that may contribute to the differences in prognosis between the high and low-risk groups. B: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis illustrating the predominant enrichment of DEGs in critical pathways such as B cell receptor signaling, osteoclast differentiation, and primary immunodeficiency, suggesting potential mechanisms through which these genes influence osteosarcoma progression and patient outcomes

3.4 Immune infiltration analysis between high and low-risk groups

We analyzed the level of immune cell infiltration between high and low-risk groups. The results in Figs. 7A-B showed that compared to the high-risk group, the low-risk group exhibited higher levels of immune cell infiltration and enrichment in immune-related pathways, such as APC co-inhibition, APC co-stimulation, B cells, CCR, CD8 + T cells, checkpoint, cytolytic activity, and DCs. Heatmaps in Figs. 7C-E displayed the differences in immune cell and immune-related pathway enrichment levels between high and low-risk groups, based on the expression levels of TERT, VEGFA, and PPARG.

Fig. 7

Immune Infiltration Analysis Between high and low-risk groups. A-B: Comparative analysis of immune cell infiltration and immune pathway enrichment between high and low-risk groups, showing higher levels of immune cell presence and pathway activity (e.g., APC co-inhibition, CD8 + T cells, cytolytic activity) in the low-risk group, suggesting a more robust immune response against the tumor. C-E: Heatmaps displaying the differential enrichment levels of immune cells and immune-related pathways based on the expression of TERT, VEGFA, and PPARG, illustrating the variation in immune landscape between high and low-risk groups

3.5 Pan-cancer analysis

We further investigated the expression differences of TERT, VEGFA, and PPARG genes in pan-cancer and normal tissues, as well as the relationship of these genes with immune scores and prognosis. Results in Figs. 8A-C showed that PPARG was expressed at low levels in cancers such as BRCA, CESC, COAD, HNSC, and at high levels in KIRC, KIRP, LIHC, STAD, among others. The TERT gene was highly expressed in most cancer types, while VEGFA was expressed at low levels in BRCA, CHOL, COAD, HNSC, and at high levels in PRAD, THCA, KIRP, and others. The correlation results between TERT, VEGFA, PPARG expression levels and immune scores (Figs. 8D-F) indicated that PPARG was positively correlated with immune scores in cancers like ACC, BRCA, DLBC, GBM, and negatively correlated in BLCA, COAD, STAD, UVM. VEGFR showed negative correlations with immune scores in cancer types such as ACC, BLCA, BRCA, CESC, and positive correlations in PCPG, LGG, DLBC. In contrast, TERT was negatively correlated with immune scores in most cancer types. The correlation results between TERT, VEGFA, PPARG expression levels and prognosis (Figs. 8G-I) showed that PPARG acted as an oncogene in GBM, LGG, PAAD, LIHC, and as a tumor suppressor in KIRC, KIPAN, UVM. Similarly, TERT acted as an oncogene in GBM, LGG, PAAD, LIHC, and as a tumor suppressor in THYM. VEGFA acted as an oncogene in most cancer types, such as GBM, LGG, KIPAN, KIRP, CESC.

Fig. 8

Pan-Cancer Analysis of TERT, VEGFA, and PPARG. A-C: Expression patterns of TERT, VEGFA, and PPARG across various cancer types and their correlation with immune scores, highlighting the diverse roles of these genes in different cancer contexts, with implications for their involvement in immune regulation and cancer progression. D-F: Correlation analysis between the expression of TERT, VEGFA, PPARG, and immune scores in different cancers, showing varying patterns of positive and negative correlations, underscoring the complex interplay between these genes and the immune microenvironment. G-I: Analysis of the prognostic impact of TERT, VEGFA, and PPARG expression across multiple cancer types, revealing their oncogenic or tumor-suppressive roles in specific contexts, which could inform their potential as therapeutic targets

3.6 Molecular docking analysis

We annotated drug functions based on differentially expressed genes between high and low-risk groups (Fig. 9) and performed molecular docking analysis between TERT and the top-ranked drug Temozolomide, concluding with a binding energy of -6.8 kcal/mol.

Fig. 9

Molecular Docking Analysis and Drug Repurposing. Drug function annotation based on DEGs between high and low-risk groups, followed by molecular docking analysis between TERT and Temozolomide, identifying a significant binding energy of -6.8 kcal/mol, suggesting the potential for repurposing Temozolomide in the treatment of osteosarcoma based on molecular compatibility