Identification and verification of a polyamine metabolism-related gene signature for predicting prognosis and immune infiltration in osteosarcoma

Abstract

Background

Although an established correlation exists between tumor cell proliferation and elevated polyamine levels, research on polyamine metabolism in osteosarcoma (OS) remains limited. This study aimed to identify polyamine metabolism-related genes (PMRGs) associated with OS prognosis and develop a prognostic model, thereby offering novel insights into targeted therapies for patients with OS.

Methods

Datasets related to OS and PMRGs were sourced from publicly accessible databases. Candidate genes were initially identified through differential expression and weighted gene co-expression network analyses. Subsequently, prognostic genes were screened using univariate Cox and least absolute shrinkage and selection operator (LASSO) regression analyses, leading to the development of a risk model. Furthermore, a nomogram model was developed using variables selected through univariate Cox regression analysis. The relationship between the signature and immune landscape was also analyzed. Following the pre-processing of single-cell RNA sequencing data, a cell communication analysis was conducted based on the identified cell types. Finally, the expression levels of prognostic genes in clinical samples were verified using reverse transcription quantitative polymerase chain reaction, western blotting and immunohistochemistry.

Results

Ninety-six candidate genes were selected for univariate Cox and LASSO regression analyses, leading to the identification of eight prognostic genes: FAM162A, SIGMAR1, SQLE, PYCR1, DDI1, PAQR6, GRIA1, and TNFRSF12A. The risk model constructed from these genes demonstrated strong predictive accuracy and classified patients into two risk groups based on the median cut-off. A nomogram model was developed, incorporating the risk score as an independent prognostic factor. The high-risk cohort exhibited lower single-sample gene set enrichment analysis scores for 17 immune cell types and reduced expression levels of seven immune checkpoint-related genes. Furthermore, eight cell types were identified, among which endothelial cells, cancer-associated fibroblasts, osteoclasts, myeloid cells, and osteoblast OS cells showed significant interactions with NK/T, B, and plasma cells. Eight prognostic genes were confirmed to be overexpressed in OS tissues.

Conclusion

The identification of FAM162A, SIGMAR1, SQLE, PYCR1, DDI1, PAQR6, GRIA1, and TNFRSF12A as prognostic genes associated with PMRGs in OS provides valuable references for prognostic assessment and personalized treatment in patients with OS.

Background

Osteosarcoma (OS), the most prevalent primary malignant bone cancer, is characterized by high invasiveness, marked malignancy, a strong propensity for pulmonary metastasis, and frequent development of resistance to chemotherapy [1]. The advent of adjuvant chemotherapy has notably improved the overall 5-year survival rate of patients with localized OS, increasing it to approximately 60%. However, this survival rate decreases to only 20% in patients with recurrence or distant metastasis [2]. Clinical treatment for OS has remained stagnant for over three decades, primarily due to the significant intra- and inter-tumoral heterogeneity inherent to OS [3]. In recent years, the concept of “precision medicine” has been introduced, aiming to classify diseases by identifying biomarkers with prognostic or therapeutic response potential, thereby enabling personalized and precise disease treatment [4]. The identification of reliable prognostic markers and use of biomarkers to predict and monitor the prognosis and progression of OS are crucial steps in selecting appropriate patient populations and optimizing treatment strategies.

Polyamines, primarily putrescine, spermidine, and spermine, are small, positively charged, polycationic molecules that are widely present in various cells across organisms. These molecules are essential for processes such as protein and nucleic acid synthesis, cell differentiation and apoptosis, and free radical scavenging, all of which are essential for cell growth and proliferation [5, 6]. The dysregulation of polyamine metabolism is intricately linked to tumorigenesis. Elevated polyamine concentrations have been detected in various malignancies, including skin, lung, and prostate cancers [7]. The polyamine pathway functions as a downstream target for numerous oncogenes. In cancer cells, the overactivation of Myc, Raf/MEK/ERK, PI3K, and mTORC1 signaling pathways induces the overexpression of critical rate-limiting enzymes in polyamine biosynthesis, ensuring continuous polyamine production to support rapid proliferation [8, 9]. Owing to its close interaction with oncogenes, polyamine metabolism has emerged as an important target for anticancer therapies. Competitive inhibitors targeting key enzymes in polyamine synthesis have demonstrated promising anticancer effects, with some advancing into clinical trials [10]. Previous studies have demonstrated alterations in polyamine metabolism during OS. A metabolomic study of serum and urine samples from patients with OS revealed increased levels of putrescine in the urine, reflecting enhanced polyamine metabolism [11]. This change may be associated with the aggressiveness, growth rate, and treatment response of OS [12]. In vitro and in vivo experiments have demonstrated that inhibition of polyamine metabolism suppresses OS growth by inducing apoptosis [13]. The use of polyamine metabolic inhibitors in combination with chemotherapeutic agents, such as methotrexate (MTX) and doxorubicin, significantly enhances their antitumor effects [14]. Furthermore, a complex interplay exists between polyamine and glutamine metabolism, as well as epigenetic regulation, in OS. Glutamine metabolism and epigenetic regulation are closely linked to bone tumor progression [12]. Consequently, changes in polyamine metabolism in OS under high mechanical stress microenvironments warrant further exploration. As research progresses, dysregulation of polyamine metabolism and its immunosuppressive effects on the tumor microenvironment (TME) have also been identified [15]. Tumor cells and immunosuppressive myeloid cells consume significant amounts of polyamines and their metabolic precursor, l-arginine, exerting immunosuppressive effects that hinder T cell proliferation and activation [16, 17]. Polyamines also aid tumor cells in adapting to the hypoxic conditions of the TME [18], and their anti-inflammatory properties facilitate immune evasion by tumor cells [19]. The antitumor potential of polyamine-targeted therapies and the extensive role of polyamines in the TME suggest that identifying biomarkers associated with polyamine metabolism may advance personalized diagnostic, prognostic, and therapeutic strategies for OS. However, no studies have explored the expression patterns of polyamine metabolism-related genes (PMRGs) in OS, nor have polyamine metabolism-related prognostic models been developed. The effect of polyamine metabolism on the TME in OS remains unclear.

This study investigated the predictive value of PMRGs in relation to the prognosis and immune infiltration status of OS, aimed at understanding the effects of key genes on OS, and constructed a prognostic PMRG risk signature. In addition, the relationships between key genes and various cells in the OS TME were explored by revealing the immune landscape. This research offers novel insights into prognostic evaluation and targeted therapy development for patients with OS, offering a robust theoretical foundation for future investigations in this field.

Materials and methods

Data retrieval

Raw data were initially sourced from the TARGET database (https://ocg.cancer.gov/programs/target), particularly including 87 OS tumor tissue samples. Of these, 85 samples contained complete survival and clinical information along with RNA expression profiles. Additional OS-related datasets, GSE99671 (GPL20148) and GSE21257 (GPL10295), were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), along with the single-cell RNA sequencing (scRNA-seq) dataset GSE162454 (GPL24676). The GSE99671 dataset included 18 tumoral bone samples and 18 paired non-tumoral bone samples, with 3 samples collected before and 15 after chemotherapy. The GSE21257 dataset comprised 53 biopsy samples of OS taken before chemotherapy. The GSE162454 dataset contained tissue samples from six patients with primary OS who had not received neoadjuvant chemotherapy. The Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/) provided the REACTOME_METABOLISM_OF_POLYAMINES gene set, consisting of 59 PMRGs.

Weighted gene co-expression network analysis

In the GSE99671 dataset, PMRG scores for OS and control samples were initially calculated using the single-sample gene set enrichment analysis (ssGSEA) algorithm from the ‘GSVA’ package (ver. 1.46.0) [20], based on their expression levels. The Wilcoxon test was used to compare differences in PMRG scores between OS and control samples (p value < 0.05). After excluding genes with zero expression values, 18,603 genes in the GSE99671 expression matrix were used to construct the weighted gene co-expression network analysis (WGCNA) network.

WGCNA was performed using the R package ‘WGCNA’ (ver. 1.46.0) to identify key module genes associated with PMRGs [21]. Clustering analysis was conducted using the ‘Hclust’ function to detect any outlier samples. The ‘pickSoftThreshold’ function was utilized to analyze network topology and determine the optimal soft threshold by calculating the scale-free fitting index and mean connectivity (ranging from 1 to 20). A minimum of 50 genes per module was set, and a hierarchical clustering tree was generated to identify the co-expression modules. Pearson’s correlation analysis was then applied to calculate correlation coefficients and corresponding p values between these modules and PMRG scores (|cor|> 0.3, p < 0.05). Key modules were identified based on their consistent positive or negative correlations with PMRG scores, and the genes within these key modules were designated as key module genes.

Differential expression analysis

Differential expression analysis between OS and control samples in the GSE99671 dataset was performed using the R package ‘DESeq2’ (ver. 1.42.0), with the criteria of |Log2 fold change (FC)|> 1 and an adjusted p value < 0.05 to identify differentially expressed genes (DEGs) [22]. Candidate genes associated with both PMRGs and OS were then identified by intersecting the key module genes with these DEGs.

Functional enrichment analysis of candidate genes

To explore the functions and pathways associated with the candidate genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the ‘clusterProfiler’ package in R (ver. 4.7.1.003), with significance defined as p value < 0.05 [23].

Establishment of a risk model

Several analyses were performed to develop a risk model for predicting overall survival in patients from the TARGET-OS dataset. Univariate Cox regression analysis, performed using the ‘coxph’ function from the R package ‘survival’ (ver. 3.5-3), identified genes significantly associated with OS (hazard ratio [HR] ≠ 1, p value < 0.05). The proportional hazards (PH) assumption test confirmed that the model assumptions were satisfied (p value > 0.05). Next, least absolute shrinkage and selection operator (LASSO) regression analysis was applied using the ‘glmnet’ package (ver. 4.1-4) to refine the selection of prognostic genes [23]. A risk model was subsequently constructed based on these identified prognostic genes, with each patient’s risk score calculated as the sum of the product of the regression coefficients (coef) and the relative expression levels (expr) of each gene: $\text{risk score}=\sum_{\text{i}=1}^{\text{n}}\text{coef }({\text{gene}}_{\text{i}})\times \text{expr}{ (\text{gene}}_{\text{i}})$. Using the median risk score as a threshold, patients in the TARGET-OS dataset were classified into the high-risk and low-risk groups. The prognostic performance of the risk model was evaluated using Kaplan–Meier (K–M) survival analysis and the log-rank test (p value < 0.05). Additionally, receiver operating characteristic (ROC) curves were generated to assess the model’s predictive accuracy, with an area under the curve (AUC) ≥ 0.6 considered indicative of satisfactory predictive ability. Finally, to validate the robustness and generalizability of the risk model, its performance was evaluated using an independent dataset, GSE21257.

Correlation analysis of clinicopathological features

To elucidate the relationships between various clinicopathological factors (age, sex, metastasis status, and stage) and risk scores in the TARGET-OS dataset, Wilcoxon rank-sum tests were conducted to assess differences in risk scores across diverse clinical subgroups (p value < 0.05). Additionally, variations in the expression of prognostic genes within these clinical subgroups were analyzed (p value < 0.05).

Construction of the nomogram model

Clinical data from patients with OS were analyzed to determine whether risk scores served as the sole independent prognostic factor. Owing to the limited sample size for the stage variable, it was excluded from subsequent analyses. Initially, univariate Cox regression analysis (p value < 0.05) and PH assumption tests (p value > 0.05) were performed using the ‘Survival’ package (ver. 3.5-3) to identify independent prognostic factors among risk scores and other clinicopathological features (age, sex, and metastasis status). Subsequently, a nomogram model incorporating these independent prognostic factors was developed using the R package ‘rms’ (ver. 6.5-0). To evaluate the predictive ability of the model, calibration and ROC curves were generated, followed by decision curve analysis (DCA).

Effects of prognostic genes on cell proliferation

To investigate the effects of prognostic genes on OS cell proliferation and explore their potential as therapeutic targets, gene dependency in OS cell lines was analyzed using the R package ‘DepMap’ (ver. 1.15.0), a resource for genome-wide CRISPR-Cas9 knockout and RNAi knockdown screening across hundreds of cancer cell lines. A negative value indicates that gene knockdown inhibits cell survival and proliferation.

Gene set enrichment analysis

To investigate enriched pathways in the high- and low-risk cohorts, GSEA was performed using the ‘clusterProfiler’ package in R. Differential expression analysis was first conducted using the ‘DESeq2’ package, where the log2FC of DEGs was calculated and ranked from largest to smallest to create a gene list corresponding to the two risk cohorts. The KEGG gene set (c2.cp.kegg.v2023.1. Hs.symbols.gmt) was used as the background for GSEA, with the enrichment criteria set at |NES|> 1 and an adjusted p value < 0.05.

Immune infiltrating analysis and drug sensitivity analysis

The ssGSEA scores for 28 immune cells in the two risk cohorts were calculated using the ssGSEA algorithm within the ‘GSVA’ package. Subsequently, the Wilcoxon rank-sum test was used to compare differences in immune cell populations and the expression of immune checkpoint-related genes between the two cohorts in the TARGET-OS dataset (p value < 0.05). The immune checkpoint-related genes analyzed included BTLA, CD200R1, HAVCR2, ICOS, LAIR1, LGALS9, PDCD1LG2, TNFSF14, TNFSF15, and TNFSF4 [24]. Relationships between differentially abundant immune cells were further assessed using Spearman’s correlation analysis (|cor|> 0.3, p value < 0.05).

Additionally, drug sensitivity in high- and low-risk patients was evaluated by calculating half-maximal inhibitory concentrations (IC50s) using the ‘oncoPredict’ package (ver. 0.2) in R [25]. Differences in IC50 values between the two cohorts were determined using the Wilcoxon signed-rank test (p value < 0.05).

Pre-processing of scRNA-seq data

Single-cell transcriptomes of OS samples were downloaded from GSE162454, and the ‘Seurat’ package (Ver. 4.3.0) was used to create Seurat objects from the single-cell data [26]. Cells with < 200 genes, genes detected in < 3 cells, cells expressing < 101 genes, > 6000 genes, or with > 20% mitochondrial genes were filtered out. Data were normalized using the LogNormalize method. Subsequently, 2000 highly variable genes were identified using the ‘FindVariableFeatures’ function with the vst method. The ScaleData and RunPCA functions were then applied to standardize the data following principal component analysis (PCA), with the optimal number of linear dimensions set to 20.

Cell annotation and cell communication

Clustering analysis was performed using the ‘FindNeighbors’ and ‘FindClusters’ functions, with the resolution parameter set at 0.1 to identify distinct cell clusters. A uniform manifold approximation and projection (UMAP) technique was then employed to visualize these cell clusters. For cell annotation, marker genes cited in the literature were used to assign cell types [27]. In the GSE162454 dataset, the interactional relationships among cell types were elucidated using the ‘createCellChat’ function from the R package ‘CellChat’ (ver. 1.6.1) to create a ‘CellChat’ object. The ‘computeCommunProb’, ‘computeCommunProbPathway’, and ‘aggregateNet’ functions were subsequently employed to construct the cell communication network [28]. The expression patterns of prognostic genes across different cell types were visualized using UMAP and illustrated using violin plots.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

RNA was extracted from ten samples, consisting of five pairs of OS and normal bone tissues, collected from the Shanghai Sixth People’s Hospital Affiliated to Shanghai Jiao Tong University School of Medicine. The inclusion criteria were as follows: (1) confirmed diagnosis of OS based on postoperative histopathology and (2) lesions located in the extremities. The exclusion criteria were as follows: (1) comorbidity with other diseases such as fever, infection, and lymph node disease and (2) recurrent patients. The detailed clinical information is provided in Supplementary Table 1. This study was performed with informed consent from the donors and was approved by the Ethics Committee (approval number: 2021-133). Total RNA was isolated using TRIzol reagent (Ambion, USA) according to the manufacturer’s protocol, and RNA concentration and quality were assessed via a NanoPhotometer N50. Subsequently, cDNA was synthesized using a SureScript First-Strand cDNA Synthesis Kit (Servicebio, China). Relative mRNA expression levels were quantified using RT-qPCR on a CFX Connect Thermal Cycler (Bio-Rad, USA), with GAPDH as the internal reference gene. The data were analyzed using the 2^−ΔΔCT method, with detailed primer information provided in Supplementary Table 2.

Western blotting

Total protein was extracted from eight pairs of fresh tissues preserved in liquid nitrogen. Tissue samples were lysed in radioimmunoprecipitation assay lysis buffer (Servicebio, China) with a cocktail of proteinase and phosphatase inhibitors. The protein concentrations were determined using bicinchoninic acid assay kit (Solarbio, China). Equal amounts of protein samples were separated by 8% or 10% sodium dodecyl sulfate‒polyacrylamide gel electrophoresis and transferred to polyvinylidene difluoride membranes (Millipore, Germany). The membranes were blocked with 5% skim milk powder for 1h. The membranes were incubated with primary antibodies (anti-DDI1: 1:1000 dilution, Bioss, China, catalog number: bs-14212R; anti-FAM162A: 1:1000 dilution, Affinity, USA, catalog number: DF15168; anti-GRIA1: 1:2000 dilution, Cusabio, USA, catalog number: CSB-PA002723; anti-PAQR6: 1:3000 dilution, GeneTex, USA, catalog number: GTX122692; anti-PYCR1: 1:2000 dilution, Affinity, USA, catalog number: DF12051; anti-SIGMAR1: 1:2000 dilution, Cusabio, USA, catalog number: CSB-PA209031; anti-SQLE: 1:2000 dilution, Bioss, China, catalog number: bs-4218R; anti-TNFRSF12A: 1:3000 dilution, Affinity, USA, catalog number: DF8023; anti-β-actin: 1:25,000 dilution, Proteintech, China, catalog number: 66009-1-Ig) at 4 °C overnight. The membrane was washed with TBST three times for 15 min each, followed by different secondary antibodies (HRP-conjugated goat anti-rabbit IgG: 1:3000 dilution, Servicebio, China, catalog number: GB23303; HRP-conjugated goat anti-mouse IgG, 1:5000 dilution, Servicebio, China, catalog number: GB23301) for 1 h at room temperature. After washing the membrane, the bands were detected with an enhanced chemiluminescence detection system using a multifunctional digital imager to obtain images.

Immunohistochemistry

Five pairs of tissues were embedded. Slides were dewaxed in xylene and rehydrated through 100%, 95%, 80% and 70% gradient alcohol. After rinsing with PBS, the sections were incubated in citrate buffer in a microwave for 15 min to facilitate antigen retrieval. After cooling to room temperature, the residual buffer was discarded. Endogenous peroxidase activity was quenched by treatment with a 3% hydrogen peroxide solution. The sections were then blocked in 5% BSA for 30 min. After incubation with primary antibodies (anti-DDI1: 1:100 dilution, Bioss, China, catalog number: bs-14212R; anti-FAM162A: 1:50 dilution, Affinity, USA, catalog number: DF15168; anti-GRIA1: 1:100 dilution, Cusabio, USA, catalog number: CSB-PA002723; anti-PAQR6: 1:100 dilution, GeneTex, USA, catalog number: GTX122692; anti-PYCR1: 1:50 dilution, Affinity, USA, catalog number: DF12051; anti-SIGMAR1: 1:100 dilution, Cusabio, USA, catalog number: CSB-PA209031; anti-SQLE: 1:100 dilution, Bioss, China, catalog number: bs-4218R; anti-TNFRSF12A: 1:50 dilution, Affinity, USA, catalog number: DF8023) at 4 °C overnight, the slides were incubated for 30 min with secondary antibodies (Servicebio, China). The color was developed using a Diaminobenzidine chromogenic solution. Then the slices were stained with hematoxylin, dehydrated, transparentized, and sealed using neutral gum. FAM162A, SIGMAR1, SQLE, PYCR1, DDI1, PAQR6, GRIA1, and TNFRSF12A are positively stained with the presence of yellowish-brown fine particles in the cytoplasm of the cells. Five 400 × high-power fields were randomly selected for each section, the percentage of positive cells were scored in each field, and the average value was calculated.

Quantification and statistical analysis

R software (version 4.3.1) and GraphPad Prism (Version 10.1.0) contributed to the statistical work of this study. Wilcoxon’s test was employed to determine inter-group differences. Student’s t test was used to compare the differences between the two groups. A p value < 0.05 was considered statistically significant.

Results

Ninety-six candidate genes associated with OS and PMRGs were identified

Initially, discrepancies in PMRG scores between OS and control samples in the GSE99671 dataset were compared, revealing that the PMRG score in the OS group was significantly higher than that in the control group (Fig. S1). Cluster analysis confirmed the absence of outliers in the GSE99671 dataset (Fig. 1A). The optimal soft threshold was determined to be 6, where the R² value approached approximately 0.85, and the mean connectivity was close to zero (Fig. 1B). This analysis led to the identification of 25 co-expression modules (Fig. 1C), with MEblue (cor = 0.93, p value < 0.001) and MEpink (cor = − 0.76, p value < 0.001) showing the strongest positive and negative correlations with PMRGs scores, respectively. These modules were designated as key modules. The intersection of 639 genes in MEpink and 3,268 genes in MEblue resulted in 3907 key module genes (Fig. 1D).

Fig. 1

Identification and functional analysis of candidate genes associated with OS and PMRGs. A Sample clustering analysis of the GSE99671 dataset. B Screening of the optimal soft threshold. C Cluster dendrogram and module overview. D Correlations between merged modules and OS. E–F Volcano plot (E) and heatmap (F) showing DEGs between OS and normal tissues. G–H GO (G) and KEGG (H) enrichment analysis of the candidate genes. OS, osteosarcoma; PMRGs, polyamine metabolism-related genes; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes

Subsequently, differential expression analysis identified 711 DEGs associated with OS, comprising 198 upregulated and 513 downregulated genes (Fig. 1E–F). By intersecting the key module genes with the DEGs, 96 genes associated with both OS and PMRGs were identified and considered candidate genes (Fig. S2).

Functional exploration of candidate genes

GO and KEGG analyses were performed to elucidate the functions associated with candidate genes. These genes were enriched in 606 GO terms, including 449 biological processes, 29 cellular components, and 128 molecular functions. GO analysis revealed that the candidate genes might play roles in the body’s response to hypoxic conditions, the synthesis and metabolism of organic and amino acids, various REDOX reactions, and membrane pressure differentials. Additionally, their functions were intricately linked to adrenaline secretion (Fig. 1G). Furthermore, candidate genes were enriched in 15 KEGG pathways. Consistent with the GO analysis, these pathways were involved in amino acid synthesis and metabolism, as well as the synthesis of folate and unsaturated fatty acids. Notably, these genes were also implicated in the HIF-1 signaling pathway, which is closely related to the response to hypoxic environments (Fig. 1H). In summary, the primary functions of the candidate genes are centered on the response to hypoxia, as well as the synthesis and metabolism of organic and amino acids.

Performance of the risk model remained robust

Prognostic genes associated with OS prognosis were initially identified using univariate Cox and LASSO regression models. Univariate Cox regression analysis identified 11 genes, of which squalene epoxidase (SQLE), family with sequence similarity 162 member A (FAM162A), sigma-1 receptor (SIGMAR1), fatty acid desaturase 2, progestin and adipoQ receptor 6 (PAQR6), pyrroline-5-carboxylate reductase-1 (PYCR1), coiled-coil domain containing 158, glutamate receptor 1 (GRIA1), cystathionine beta-synthase, and DNA damage-inducible 1 (DDI1) were associated with an increased risk of OS (HR > 1, p value < 0.05), whereas tumor necrosis factor receptor superfamily member-12A (TNFRSF12A) was associated with a reduced risk of OS (HR < 1, p value < 0.05). The PH assumption test confirmed that the proportional hazards assumption in the univariate regression model was valid (Fig. 2A). When the minimum lambda value was set to 0.0495 through tenfold cross-validation, eight prognostic genes were identified in the LASSO model: FAM162A, SIGMAR1, SQLE, PYCR1, DDI1, PAQR6, GRIA1, and TNFRSF12A (Fig. 2B–C). The risk coefficients for these prognostic genes are provided in Table 1. The risk score equation was calculated as follows:

Fig. 2

Identification of genes associated with OS prognosis and establishment of prognostic risk models. A Univariate Cox regression analyses of gene expression and overall survival of OS. B LASSO regression coefficient path plot. C Cross-validation curve for LASSO regression. D–E Risk score (D) and overall survival (E) status in the training cohort. F Heatmap of the distribution of the eight genes showing the expression differences of PMRGs in the high- and low-risk groups. G Difference in the survival outcomes between the two groups. H ROC curve showing the prediction accuracy of the risk signature. OS, osteosarcoma; PMRGs, polyamine metabolism-related genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver operating characteristic

Table 1 Lasso regression coefficients

$$\begin{aligned} {\text{risk score }} & = SQLE \times \, 0.{381}0{444 } + FAM162A \times \, 0.{5238747 } + SIGMAR1 \times \, 0.{1131558 } \\ & \quad + PAQR6 \times \, 0.{17598}0{1 } + PYCR1 \times \, 0.{1496255 } + GRIA1 \times \, 0.{3193388 } \\ & \quad + TNFRSF12A \times \, \left( { - 0.{1226966}} \right) \, + DDI1 \times \, 0.{77}0{2}0{21}. \\ \end{aligned}$$

In the TARGET-OS dataset, 43 high-risk and 42 low-risk patients with OS were classified based on a median risk score of 4.213956 (Fig. 2D). Patients with higher risk scores had a greater likelihood of mortality, indicating a potential correlation between elevated risk scores and adverse outcomes (Fig. 2E). Additionally, TNFRSF12A was predominantly expressed in the low-risk cohort, whereas the seven other prognostic genes were more prominently expressed in the high-risk cohort (Fig. 2F).

Kaplan–Meier survival curves further highlighted that patients in the low-risk cohort had longer survival times (Fig. 2G). The AUC for 1-, 3-, and 5-year survival were 0.872, 0.817, and 0.799, respectively, demonstrating the robust performance of the risk model (Fig. 2H). Patients with OS in the GSE21257 dataset were similarly classified into 27 high-risk and 26 low-risk individuals, based on a median risk score of 17.6239508. Consistent with the findings from the TARGET-OS dataset, results from the GSE21257 dataset further supported the validity of the risk model (Fig. S3A–E).

Discrepancies of risk score and prognostic genes’ expression in clinicopathological subgroups

Among the clinicopathological subgroups, the risk score was notably higher in the metastatic subgroup (p < 0.05), whereas no significant differences were observed in other clinicopathological subgroups (Fig. S4A–D). Additionally, DDI1 and PAQR6 were significantly elevated in patients aged < 18 years and > 18 years, respectively (Fig. S4E). No significant differences were found in the expression of prognostic genes between sex subgroups (Fig. S4F). Notably, SIGMAR1 exhibited elevated expression in both the metastatic and stage 1/2 subgroups, whereas a similar expression pattern was observed for FAM162A in the metastatic subgroup (Fig. S4G–H).

Predictive efficacy of the model incorporating all prognostic genes was calculable

To determine whether the risk score served as an independent prognostic factor for OS, univariate Cox regression analysis and the PH assumption test were conducted, integrating sex, age, metastasis status, and risk score. A significant association was observed between the risk score and overall survival in OS (Fig. 3A). Consequently, a nomogram model incorporating all prognostic genes was developed. Each patient with OS was assigned a total point score, with a higher score indicating a lower probability of survival at the 1-, 3-, and 5-year intervals (Fig. 3B). Notably, the calibration curves, with 45° lines corresponding to the 1-, 3-, and 5-year intervals, suggested that the nomogram model achieved an ideal prediction (Fig. 3C). The AUCs for the 1-, 3-, and 5-year intervals were 0.8044, 0.8025, and 0.8005, respectively, further demonstrating the predictive accuracy of the model (Fig. 3D). Additionally, DCA indicated that the net benefits of the nomogram model exceeded those of other individual factors, underscoring the model’s strong predictive performance (Fig. 3E).

Fig. 3

Construction of the nomogram predicting patients’ survival in OS. A Univariate Cox regression analysis of clinical characteristics and risk scores. B Construction of the nomogram for survival prediction of patients with OS. C Calibration curve of the nomogram. D ROC curves showing the prediction accuracy of the nomogram. E DCA of the nomogram. OS, osteosarcoma; ROC, receiver operating characteristic; DCA, decision curve analysis

Relationship between prognostic genes and OS cell lines

Cancer-dependency maps are frequently used to identify the reliance of cancer cells on specific genes or genomic alterations. In this study, the dependencies of prognostic genes in OS cell lines were assessed using RNA interference screening from cancer dependency maps. Notably, the knockout of DDI1 and GRIA1 enhanced cell line proliferation, whereas the knockout or downregulation of the other six prognostic genes inhibited cell proliferation (Fig. 4A).

Fig. 4

Differences in molecular and immune infiltration characteristics between the high- and low-risk groups. A Validating the impact of prognostic genes on OS cell proliferation based on the cancer dependency map. B GSEA of KEGG in the high- and low-risk groups. C Heatmap depicting immune cell infiltration in the high- and low-risk groups. D Differences in the levels of 28 types of immune cell infiltration between the high- and low-risk groups. E Correlation of 28 types of immune cells. F Differences in gene expression of major immune checkpoints between the high- and low-risk groups. *: p < 0.05; **: p < 0.01; ***: p < 0.001; ****: p < 0.0001, ns: not significant. OS, osteosarcoma; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes

Functional annotation of risk score related to two risk cohorts

GSEA was employed to explore the KEGG-related pathways and the underlying mechanisms associated with the genes in the high- and low-risk cohorts. In total, 74 KEGG pathways were identified, including those related to the one-carbon pool by folate, RNA degradation, autoimmune thyroid disease, and other pathways (Fig. 4B).

Immune system of high-risk patients might be relatively underactive

Immune cell infiltration in the two risk cohorts was analyzed, with Fig. 4C presenting the ssGSEA scores for 28 immune cell types across both groups. Seventeen immune cell types exhibited significant differences between cohorts, all showing higher ssGSEA scores in the low-risk cohort. These included activated B cells, CD56bright NK cells, central memory CD4( +) T cells, central memory CD8( +) T cells, effector memory CD8( +) T cells, immature B cells, macrophages, monocytes, myeloid-derived suppressor cells, NK cells, NK T cells, neutrophils, plasmacytoid dendritic cells (DCs), regulatory T cells, T follicular helper cells, T helper type 1 cells, and T helper type 17 cells (Fig. 4D). Furthermore, significant correlations were identified among most of the discrepant immune cells, except for T helper type 17 cells, NK T cells, plasmacytoid DCs, and CD56bright NK cells (|cor|> 0.3, p value < 0.05) (Fig. 4E). Among the 10 immune checkpoint-related genes examined, the expression levels of PDCD1LG2, TNFSF4, BTLA, CD200R1, HAVCR2, LAIR1, and TNFSF15 were significantly downregulated in the high-risk cohort (p value < 0.05) (Fig. 4F). Furthermore, drug sensitivity analysis between the two risk cohorts revealed that among the top 20 drugs (NSC319726_461, Torin 2_371, CD532_449, NPK76_II_72_1_257, SB52334_304, SB505124_476, TL_2_105_211, kb NB 142_70_407, KIN001_270_345, BMS_345541_203, CX_5461_300, PHA_793887_301, Tubastatin A_265, Dasatinib_51, WH_4_023_56, TGX221_94, AZD8186_1444, WZ_1_84_59, HG6_64_1_159 and NSC_207895_269) with significantly different IC50 values, Dasatinib_51, WH_4_023_56, TGX221_94, AZD8186_1444, WZ_1_84_59, and HG6_64_1_159 exhibited higher IC50 values in the high-risk cohort (Fig. S5).

Eight cell types were merged by 12 cell clusters

Following quality control, 24,611 genes were identified as expressed in 42,786 cells from the GSE162454 dataset. The distribution of nFeature RNA, nCount RNA, and the percentage of mitochondrial genes were visualized using violin plots (Fig. S6A). Fig. S6B highlights 2000 highly variable genes, marked in red. UMAP dimensionality reduction was performed using the top 20 principal components (Figs. S6C, 5A–B). Subsequent UMAP analysis and cell annotation allowed for the integration of 12 cell clusters into eight distinct cell types based on the expression patterns of marker genes within these clusters. These cell types included osteoblastic OS cells, osteoclasts (OCs), myeloid cells, cancer-associated fibroblasts (CAFs), NK/T cells, plasma cells, endothelial cells, and B cells (Fig. 5C–E).

Fig. 5

Identification of eight cell types in OS based on single-cell RNA-seq data. A–B PCA identified the top 20 PCs at P < 0.05. C UMAP plot of single cells illustrating 12 classified cell clusters. D Marker genes for each cell type. E UMAP plot of eight cell types in OS. OS, osteosarcoma; PCA, principal component analysis; PCs, principal components; UMAP, uniform manifold approximation and projection

Analysis of relationships among eight cell types

Analysis of the interactional relationships among these cell types revealed that endothelial cells, CAFs, OCs, myeloid cells, and osteoblastic OS cells exhibited a high number of significant interactions with NK/T cells, plasma cells, and B cells (Fig. 6A–B). Further examination of prognostic gene expression across the eight cell types indicated that in the GSE162454 dataset, DDI1 was not detected. However, SQLE, FAM162A, SIGMAR1, TNFRSF12A, and PYCR1 showed high expression levels in all cell types. Additionally, PAQR6 was primarily expressed in CAFs and osteoblastic OS cells, whereas GRIA1 was almost exclusively expressed in osteoblastic OS cells (Figs. 6C–I, S7).

Fig. 6

Intercellular interactions within osteosarcoma microenvironment. A–B Interconnections and nodes of cell–cell communication. C–I UMAP plot showing the expression patterns of prognostic genes across various cell types. OS, osteosarcoma; UMAP, uniform manifold approximation and projection

Prognostic genes were overexpressed in OS

RT-qPCR results confirmed that, except for GRIA1 and SIGMAR1, which did not show significant differences between OS and normal samples (p value > 0.05), the remaining six prognostic genes (FAM162A, SQLE, PYCR1, DDI1, PAQR6, and TNFRSF12A) were significantly overexpressed in OS (p value < 0.05) (Fig. 7A). Western blotting analysis further confirmed that the protein expressions levels of all eight prognostic genes were markedly elevated in OS tissues compared with those in normal tissues (Fig. 7B–C). Consistent with these findings, immunohistochemical staining demonstrated strong positive signals for FAM162A, SIGMAR1, SQLE, PYCR1, DDI1, PAQR6, GRIA1, and TNFRSF12A in OS tissues (Fig. 8A–B).

Fig. 7

The mRNA and protein expression of prognostic genes in OS samples. A The mRNA expression of prognostic genes relative to that of β-actin in OS and normal tissues. B–C The protein expression of prognostic genes relative to that of β-actin in OS and normal tissues. *: p < 0.05; **: p < 0.01; ****: p < 0.0001, ns: not significant. OS, osteosarcoma

Fig. 8

IHC verification of the expression level of prognostic genes in OS and normal tissue. *: p < 0.05; **: p < 0.01. OS, osteosarcoma

Discussion

The mechanisms underlying the occurrence and progression of OS are highly complex, and no global consensus on its etiology has yet been reached. OS is characterized by a high incidence of chemotherapy resistance and early hematogenous metastasis, both of which significantly affect patient survival [29]. Developing an effective prognostic model for accurately stratifying patients with OS is essential for optimizing treatment [30]. As the volume of OS-related data in public databases continues to grow, researchers have more opportunities to utilize open gene expression datasets to develop predictive tools based on specific genes. Polyamines, naturally occurring organic polycations, play vital biological roles in both eukaryotic and prokaryotic organisms. In normal cells, polyamine homeostasis is tightly regulated via feedback mechanisms within the polyamine metabolic network. However, this homeostasis is frequently disrupted in tumor cells [31]. Tumor cells synthesize polyamines that can enter the bloodstream and are excreted in urine, making polyamine levels in the urine and serum important biomarkers for monitoring tumor progression in patients with cancer [32]. Targeting polyamine metabolism, either by inhibiting polyamine synthesis or promoting its degradation, can deplete intracellular polyamines, thereby inhibiting tumor cell proliferation and inducing apoptosis [33]. Consequently, the polyamine metabolic pathway is gaining attention as a potential target for tumor prognosis assessment and treatment. While prognostic models based on PMRGs have been developed for lung, renal, and colorectal cancers [34,35,36], studies exploring the relationship between polyamine metabolism and OS remain limited.

In this study, multiple OS datasets from the TARGET and GEO databases were analyzed using bioinformatics techniques, leading to the identification of 96 candidate genes associated with both OS and PMRGs. Using LASSO-Cox regression analysis, eight genes were identified as strong predictors of OS prognosis, and a prognostic risk model was developed based on these genes. The model demonstrated high predictive accuracy for 1-, 3-, and 5-year survival probabilities, and its robustness and reliability were further confirmed through validation using independent testing datasets.

Most of the genes included in this model have been previously identified as prognostic biomarkers for various malignancies. Huang et al. developed a glycolysis-related gene risk model incorporating FAM162A, which effectively assessed the prognosis of OS [37]. Liu et al. identified 18 target genes associated with the prognosis of patients with bladder cancer, including PYCR1 [38]. Additionally, copy number variations in PAQR6 have been recognized as independent prognostic factors in patients with muscle-invasive bladder cancer [39]. From a molecular perspective, FAM162A acts as a HIF-1α-responsive protein, with its expression significantly elevated in cervical cancer cells. This protein potentially promotes tumorigenesis by influencing hypoxia-induced signaling pathways and mitochondrial apoptosis [40]. SIGMAR1 regulates membrane electrical activity in response to extracellular matrix stimuli, thereby promoting angiogenesis and enhancing tumor cell invasiveness [41]. SQLE plays a pivotal role in cholesterol metabolism, facilitating pancreatic cancer growth by alleviating endoplasmic reticulum stress, inhibiting apoptosis, and activating the Src/PI3K/Akt signaling pathway [42]. Consistent with our findings, Wang et al. observed significantly elevated SQLE expression in OS tissues, correlating with chemotherapy resistance and poor prognosis. Their experiments confirmed that knocking down SQLE in OS cells reduced cholesterol levels, suppressed the PI3K/Akt pathway, inhibited proliferation, and enhanced chemosensitivity, reinforcing the validity of our model [43]. PYCR1, a key enzyme in proline synthesis, is crucial for collagen production in the extracellular matrix of CAFs. Knocking down PYCR1 in CAFs reduces tumor collagen production, inhibits tumor growth, and suppresses metastasis [44]. PAQR6, a member of the PAQR family, functions as a membrane progesterone receptor, promoting prostate cancer cell proliferation and migration via the MAPK signaling pathway [45]. GRIA1, a critical subunit of alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionate-type glutamate receptor, is overexpressed in glioblastoma cells and is associated with increased cell-surface expression of beta1 integrin, focal adhesion kinase phosphorylation, and focal adhesion complex formation, thereby promoting tumor invasion [46, 47]. Notably, TNFRSF12A exhibits dual roles in different cancers: it promotes angiogenesis, cell proliferation, and migration and induces cachexia in lung cancer [48], whereas in endometrial cancer cells, it induces apoptosis and helps overcome platinum resistance in high-grade serous ovarian cancer [49, 50]. Xu et al. reported TNFRSF12A overexpression in OS cells, highlighting its role in attenuating the anticancer effects of miR-149-5p, which contrasts with the protective role suggested by our risk model [51], Therefore, further investigation is necessary to elucidate the specific mechanisms by which TNFRSF12A functions in OS. The role of DDI1 in cancer development and progression remains largely unknown. Our findings provide valuable insights for future OS studies.

To explore the mechanisms underlying the varied prognoses in patients with OS, GSEA was performed to identify pathways enriched in the high- and low-risk groups. The analysis revealed significant enrichment of highly expressed genes in the one-carbon pool via the folate pathway. Folate is crucial for the generation of S-adenosylmethionine (SAM) and deoxythymidine monophosphate through one-carbon metabolism and the methionine cycle. SAM functions as a universal methyl donor in cellular processes and is critical for polyamine biosynthesis [52]. In humans, the reduced folate carrier (RFC) is essential for cellular folate uptake and also mediates MTX uptake [53]. Unfortunately, over 60% of OS samples exhibit decreased RFC expression at the initial biopsy [54], and mutations in the RFC gene are frequently observed in patients with OS [55]. These alterations lead to diminished RFC transport capacity, which has been linked to increased resistance to MTX treatment in OS [56]. Recent research indicates that in tumor cells with low RFC expression, cytoplasmic one-carbon metabolism mediated by SHMT1 serves as the primary source of one-carbon units [57]. A deeper understanding of folate and one-carbon unit uptake and utilization in the high-risk groups may provide insights into predicting patient sensitivity to MTX and could inform the development of novel therapeutic targets for OS.

Given the suppressive effect of polyamines on tumor immunity, immune infiltration levels were assessed in both the high- and low-risk groups. These findings revealed that the low-risk group exhibited higher levels of immune infiltration, indicating enhanced immune activity. Notably, a significant disparity in NK cell infiltration was observed between the two cohorts. As innate lymphoid cells, NK cells exhibit inherent cytotoxicity against tumor cells, thereby inhibiting their proliferation and metastasis [58]. Polyamines exert bidirectional regulatory effects on NK cells. As natural immune suppressors, polyamines can reduce NK cell cytolytic activity, facilitating immune evasion by tumor cells [59]. Additionally, polyamines inhibit the expression of NK1.1 receptors, perforin, and interferon-γ, impairing NK cell-mediated recognition of tumor cells [60]. Conversely, polyamine biosynthesis enhances interleukin-2 production, which induces NK cell proliferation and augments their cytolytic activity, thereby boosting NK cell cytotoxicity [61]. OS cells fully express adhesion molecules CD54 and CD58, rendering them susceptible to NK cell recognition and binding [62]. Furthermore, OS cells exhibit reduced expression of ligands for NK cell inhibitory receptors and increased expression of ligands for activating NK cell receptors, contributing to NK cell activation [63, 64]. Consequently, NK cells in patients with OS typically retain intact cellular activity and functionality against tumor cells [65], which may explain the improved outcomes observed in the low-risk group. Recent advances in immunotherapy have significantly enhanced the therapeutic efficacy for malignant tumors, with various immunotherapeutic approaches showing substantial clinical progress. Basic research and clinical trials focusing on immunotherapy for OS are increasing [66]. Immunological checkpoint analysis revealed higher expression levels of most checkpoints in the low-risk group, suggesting that these patients may have heightened sensitivity to immunotherapy. The ongoing development and refinement of scRNA-seq technology has provided a robust foundation for studying tumor heterogeneity [67]. By leveraging existing single-cell datasets, this study described the TME landscape of OS under the influence of polyamine metabolism and identified eight distinct cell clusters: osteoblastic OS cells, OCs, myeloid cells, CAFs, NK/T cells, plasma cells, endothelial cells, and B cells. Prognostic gene expression profiles were analyzed across these cell types, and SQLE, FAM162A, SIGMAR1, TNFRSF12A, and PYCR1 showed high expression levels across various cells. PAQR6 was predominantly expressed in CAFs and osteoblastic OS cells, whereas GRIA1 was expressed almost exclusively in osteoblastic OS cells.

This study had some limitations. First, all data were sourced from public datasets, which may limit the diversity of the data. The predictive role of the risk model in OS prognosis requires validation in prospective cohorts. Second, the correlation between risk score and polyamine concentration may not be direct, necessitating further investigation to confirm these findings. Finally, although this study validated gene expression in clinical samples, additional research is required to explore the specific mechanisms involving polyamines and PMRGs in OS pathogenesis.

Conclusion

This study identified eight prognostic genes associated with both PMRGs and OS and developed a risk model to provide insights into the prognostic assessment of patients with OS. In addition, this gene signature effectively characterizes the immune infiltration status in patients with OS, potentially providing deeper insights into the complexity of the OS TME. These novel findings have valuable clinical implications for tumor-targeting therapies for OS.