To conduct our study on sphingolipid metabolism, a total of 160 genes were retrieved from the InnateDB (https://www.innatedb.com/), which are known to be implicated in this metabolic pathway. Additionally, transcriptomic datasets along with comprehensive clinical profiles for osteosarcoma patients were procured from the TARGET initiative (https://ocg.cancer.gov/programs/target). In total, transcriptomic profiles were obtained for 88 individuals diagnosed with osteosarcoma. However, complete clinical records were available for only 85 of these patients. Consequently, transcriptomic data from the remaining three patients were excluded from the study, thereby enabling subsequent analyses to proceed with the cohort of 85 osteosarcoma patients. The U-2 OS cell line, characterized by its epithelial morphology, originates from a sarcoma of moderate differentiation. This line exhibits adherent growth characteristics and displays significant chromosomal alterations. The osteosarcoma cell line U2OS (KCB200962YJ) was procured from Kmcellbank, serving as a crucial biological model for our analyses.

To evaluate the prognostic influence of distinct cuproptosis-SPGs, univariate Cox regression analyses were undertaken, linking their differential expression with survival durations (31). Genes demonstrating a significant association with survival metrics were identified as prospective elements for subsequent LASSO regression (32, 33). Selection of the optimal penalty parameter (λ) for the model was achieved through 10-fold cross-validation, utilizing a criterion of minimality. Computation of the final risk score was based on normalized expression values and coefficients of the selected variables. The formula is presented as follows: B4GALNT1*0.0314 + SGMS2*0.1847 + ABCA2*0.3275 - B3GALT4*0.0206 -ST3GAL2*0.1409 - NPC2*0.1068 - APOE*0.0882.

Following the development of the model, individuals in the TARGET-OS cohort were stratified into high-risk and low-risk categories, centered on the median score from the training set. Differences in survival rates between these categories were assessed using Kaplan-Meier survival plots and log-rank tests, facilitated by the “survival” and “survivalminer” software packages. Additionally, the model’s predictive accuracy was evaluated over 1, 3, and 5 years, using receiver operating characteristic (ROC) analysis (34, 35), implemented via an R package.

To assess the clinical relevance of the prognostic framework, decision curve analysis (DCA) curves were generated using the “ggDCA” package. This methodology facilitates the identification of patient risk and quantifies the net benefit of the predictive model across varying threshold probabilities.

For functional enrichment assessment, the “clusterProfiler” R package was employed to analyze functional enrichment based on differentially expressed genes (DEGs) between high-risk and low-risk groups. Additionally, single-sample Gene Set Enrichment Analysis (ssGSEA) was performed for each individual in both cohorts (36). Pathways that were significantly enriched were depicted through heatmaps illustrating the average ssGSEA scores.

The algorithms CIBERSORT and ESTIMATE were employed to determine the relative abundance and infiltration levels of immune cell subtypes, respectively (37). Spearman’s method was utilized to assess the correlation between the Riskscore and the degree of immune infiltration as per ESTIMATE. Additionally, the immune scores for TIL subgroups were derived using ssGSEA. Furthermore, ssGSEA in the GSVA packages was applied to quantify immune features within each sample.

For the processing of single-cell RNA sequencing data, we employed the R package Seurat (38, 39), version 4.1.0, to facilitate the log-normalization and standardization of single-cell transcriptomic profiles (40). Rigorous quality control measures were implemented, evaluating cells based on gene expression within the count matrix and the proportion of mitochondrial gene counts. Cells expressing more than 1500 genes were excluded, along with those exhibiting over 5% mitochondrial gene count. Subsequently, data normalization across the library size was conducted using the “normalize_total” function from the Scanpy toolkit. This normalization facilitated the generation of a logarithmically transformed data matrix, which was utilized in subsequent analytical processes. Ultimately, eight distinct cell types were annotated.

To investigate variations in drug responsiveness between low-risk and high-risk cohorts, we utilized the R package “pRRophetic” (41). Subsequently, differentially expressed genes (DEGs) distinguishing these groups were integrated into the Connectivity Map to pinpoint therapeutic agents that could potentially benefit patients at higher risk. Following this, the U2OS cell line was treated with roscovitine to assess the impact of cuproptosis-SPGs on drug sensitivity via cytotoxic assays.

The impact of cuproptosis-SPGs on the migratory abilities of U2OS cells was evaluated utilizing both wound-healing and transwell migration assays (42). In the wound-healing assay, siRNA was employed to disrupt U2OS cells over a period of 72 hours. Subsequently, a deliberate “wound” was created within the cell monolayer, and cell migration rates were assessed after 48 hours. Additionally, siRNA-treated U2OS cells were placed in the upper chamber of a transwell setup. The lower chamber was supplied with 15% fetal bovine serum to serve as a chemoattractant for the U2OS cells. The number of cells present in the lower chamber after 36 hours was quantified to gauge cellular migration capabilities (43).

In this study, the analytical procedures were performed using RStudio, release 4.1.3. To assess overall survival differences, we utilized the log-rank test along with Kaplan-Meier survival analysis. For processing Western blot data, software such as ImageJ and GraphPad Prism were applied. Pearson’s method of correlation analysis was employed to execute all correlation assessments. Evaluation of disparities between the groups was conducted through the unpaired, two-tailed Student’s t-test, with a significance threshold set at P<0.05.

We conducted univariate Cox regression analysis on 158 cuproptosis-SPGs, identifying 32 of them as prognostically significant genes in osteosarcoma. Of these, 12 were associated with lower survival rates ( Figure 1A ). Subsequent to this identification, we employed LASSO regression modeling based on the same set of 32 cuproptosis-SPGs. This process culminated in the selection of seven core cuproptosis-SPGs for the construction of a prognostic model for osteosarcoma patients ( Figures 1B, C ). Multivariable Cox regression analysis elucidated the risk coefficients for these seven SPGs, categorizing patients into high and low-risk groups ( Figure 1D ). Kaplan-Meier survival plots revealed that patients in the high-risk category exhibited significantly poorer prognoses compared to those scored as low-risk ( Figure 1E ). The predictive capacity of the cuproptosis-SPGs model was evaluated using ROC curves. The model demonstrated impressive predictive accuracy with an AUC of 0.874 at two years and 0.837 at five years ( Figure 1F ).

Following the construction of the prognostic model, we further assessed its overall accuracy and credibility using the C-index ( Figure 2A ). Calibration curves additionally demonstrated the predictive capacity of the model, showing good consistency with actual outcomes over 1, 3, and 5 years ( Figure 2B ). Moreover, osteosarcoma patients classified as high risk exhibited an increased cumulative hazard ( Figure 2C ), suggesting that Cuproptosis-SPGs may serve as a potential marker for the progression of osteosarcoma patients. Decision curve analysis indicated that the drug decisions based on nomograms provided greater benefits for osteosarcoma patients compared to other individual indicators ( Figure 2D ).

The tumor microenvironment plays a pivotal role in the development and progression of cancer. In this study, patients with osteosarcoma were stratified into high-risk and low-risk categories utilizing the Cuproptosis-SPGs model. Differential gene expression analysis was conducted for both groups, followed by functional enrichment studies. This approach identified 79 differentially expressed genes (DEGs), the majority of which were upregulated ( Figure 3A ). Functional enrichment of these DEGs highlighted extensive involvement in immune response processes and sphingolipid metabolism between the high and low-risk groups ( Figures 3B, C ). Among these, four Gene Ontology (GO) processes showed significant alterations ( Figure 3D ). The copper ion binding process and sphingomyelin metabolic process are linked through the extracellular structure organization, suggesting that the extracellular structure may act as a bridge in the sphingolipid-cuproptosis network ( Figure 3E ).

To elucidate the interplay between cuproptosis metabolism and sphingolipid metabolism, correlations between genes involved in cuproptosis and sphingolipid pathways were investigated. As depicted in Figure 4A , a notable, albeit negative, correlation persists between these gene sets. Noteworthy, the correlation between FDX1 and B4GALNT1 was found to be the most significant (P < 0.01) ( Figure 4B ). Additionally, the heatmap in Figure 4C illustrates the correlation levels among cuproptosis genes, with DAL and DLD exhibiting the highest correlation coefficient (R = 0.52).

Examination at the level of individual cells has unveiled a significant presence of T cell and monocyte/macrophage subsets within the immune context of osteosarcoma ( Figures 5A, B ). By employing UMAP analysis, we identified 28 unique clusters across eight cellular subtypes, revealing diverse subpopulations of cells ( Figures 5C, D ). Employing heatmap visualization methods, distinctive patterns of interaction between osteoblasts and fibroblasts were emphasized ( Figure 5E ). To delve into the expression patterns of SRGs within the tumor microenvironment in osteosarcoma, we explored the expression levels of the top six ranked SRGs in each cell subtype ( Figures 5F, G ). A pivotal finding was the notable decrease in peptide antigen binding pathways ( Figures 6A, C ), coupled with a notable increase in ATP synthase activity in T cells ( Figures 6B, D ). Analysis of SRGs showed a consistent rise in APOE expression among different cell types ( Figure 7A ). Given the variety of T cell subsets discovered in osteosarcoma, additional studies were initiated to explore the interactions between these lymphocytes and the tumor cells ( Figures 7B–D ). Investigative efforts into the transcription factors in tumor cells highlighted potential pathways influencing differential expression ( Figure 7E ).

In patients with osteosarcoma, there exists a stark contrast in the levels of Cuproptosis-SPGs between those at high risk and those at low risk ( Figure 8A ). An estimation of cellular proportions indicates a diminished presence of CD8+ T cells in the high-risk group ( Figure 8B ). Analyses of the tumor immune microenvironment (TME) reveal significant variances not only in CD8 T cells but also in NK cells and Tregs across different risk categories of osteosarcoma patients ( Figures 8C, D ). These differences were corroborated by an analysis linking risk scores with immune cell correlations ( Figure 9A ). Given that risk scores are derived from Cuproptosis-SPGs, further investigation into the association between Cuproptosis-SPGs and the infiltration levels of CD8 T cells, NK cells, and Tregs was conducted. The results demonstrate a significant negative correlation between ABCA2 and the infiltration of CD8 T cells ( Figure 9B ), suggesting that ABCA2 may be a critical Cuproptosis-SPG gene influencing anti-tumor immunity.

The emergence of chemoresistance presents a significant challenge in the treatment of tumors. Our predictive analyses suggest that patients with high-risk osteosarcoma exhibit enhanced sensitivity to roscovitine, metformin, and bortezomib compared to those with low-risk osteosarcoma ( Figure 10 ). Given the substantial weight of ABCA2 in the risk assessment, further investigations were conducted to explore its role in mediating drug sensitivity in osteosarcoma. After transfection with siRNA, osteosarcoma U2OS cells were exposed to varying concentrations of roscovitine for 48 hours. Cytotoxicity assays revealed that, relative to controls, cells with ABCA2 knockdown demonstrated reduced cytotoxic response to roscovitine ( Figure 11A ), with more pronounced differences at concentrations between 100–200 µM ( Figure 11B ). These findings align with the predictions from the drug sensitivity model, indicating that osteosarcoma patients with elevated ABCA2 expression are more susceptible to roscovitine. Future research should elucidate the mechanisms by which ABCA2 influences drug sensitivity through its interactions with pathways such as cuproptosis and sphingolipid metabolism.

Following the validation of the drug sensitivity model using siRNA, the precision of the Cuproptosis-SPGs prognostic model was assessed. It was observed that 48 hours post siRNA treatment, the number of U2OS cells traversing the compartments was diminished ( Figure 11C ) and exhibited statistically significant differences ( Figures 11D, E ). Consistent with the transwell assay results, wound healing experiments indicated that inhibition of ABCA2 curtailed the migration distance of U2OS cells ( Figure 11F ), thus underscoring the role of ABCA2 in facilitating U2OS cell migration ( Figures 11G, H ).

Tumor cells evade apoptosis by remodeling their metabolism. We assessed the impact of the Cuproptosis-SPG gene ABCA2 on apoptosis in osteosarcoma cells using flow cytometry. It was observed that U2OS cells with disrupted ABCA2 were more susceptible to apoptosis, exhibiting an increased apoptotic rate of 16.39% compared to the control group ( Figures 11I, J ). Interestingly, after ABCA2 knockdown, autophagy levels in U2OS cells were also suppressed, and this was associated with hindered expression of ERK proteins ( Figures 11K, L ). These findings suggest that Cuproptosis-SPG may regulate the apoptotic process in osteosarcoma cells through both the autophagy and the MAPK pathways.

We further tried to explore the potential utility of Cuproptosis-SPG across additional tumor contexts. Our investigations incorporated an expression analysis of ABCA2 in normal tissues ( Supplementary Figure S1A ). The findings indicate a predominant expression of ABCA2 in immune cells, such as Gdh T cells ( Supplementary Figure S1B ). Notably, within a pan-cancer analysis, ABCA2 exhibited a significantly elevated expression in CHOL compared to normal tissues ( Supplementary Figure S1C ). Despite this, the correlation between ABCA2 expression levels and prognosis in CHOL was minimal ( Supplementary Figure S1D ).