Machine Learning–Based Identification and Validation of PYCR1 and PYGM as Prognostic Biomarkers for Osteosarcoma

Abstract

Background:

Osteosarcoma (OS) is a malignant tumor originating in the bones, predominantly affecting children and adolescents, characterized by high aggressiveness and poor prognosis. Identifying new prognostic biomarkers is crucial for improving the diagnosis and treatment of OS.

Methods:

In this study, we collected gene expression data from 88 OS samples from the UCSC Xena platform and normal tissue expression data from 396 Genotype-Tissue Expression (GTEx) samples. Prognosis-related genes were first screened by univariate Cox regression and then further selected using the Least Absolute Shrinkage and Selection Operator (LASSO) regression. Based on these candidate genes, non-negative matrix factorization (NMF) was used for molecular subtype identification, and the Kaplan-Meier analysis was applied to compare survival among subtypes. Tumor microenvironment and immune cell infiltration analyses were performed to characterize differences between risk groups. In addition, the expression patterns of key genes were validated by quantitative real-time polymerase chain reaction (qRT-PCR), hematoxylin-eosin staining, immunohistochemistry, and immunofluorescence.

Results:

Pyrroline-5-carboxylate reductase 1 (PYCR1) was consistently upregulated in OS and was associated with poor prognosis. In contrast, glycogen phosphorylase, muscle-associated (PYGM) showed analysis-level-dependent expression patterns: it was downregulated at the bulk transcriptomic and tumor cell levels compared with normal controls, whereas within the OS cohort, relatively higher PYGM expression was observed in the high-risk group. Tumor microenvironment and immune cell infiltration analyses revealed significant immune differences between high- and low-risk groups. Histological and protein-level assays further confirmed the presence and cellular localization of PYCR1 and PYGM in OS tissues.

Conclusion:

This study systematically identified and validated PYCR1 and PYGM as potential prognostic biomarkers for OS using integrated statistical and machine learning approaches. The PYCR1 showed a consistently tumor-promoting expression pattern, whereas PYGM demonstrated context-dependent expression changes across bulk tissue, risk-stratified tumor samples, and tumor cell lines, highlighting the biological complexity of metabolic biomarkers in OS.

Introduction

Osteosarcoma (OS) is a malignant bone tumor and is one of the most common primary bone malignancies.^1-3 It predominantly affects children and adolescents, with the highest incidence observed in teenagers aged 10 to 20 years. ⁴ This aggressive tumor typically arises in the metaphyses of long bones, most commonly affecting the distal femur, proximal tibia, and proximal humerus. Despite significant advancements in treatment, including surgery and chemotherapy, the prognosis for metastatic or recurrent OS remains poor. ⁵ Recent developments in immunotherapy have offered new therapeutic possibilities, providing additional benefits to patients through a deeper understanding of immune responses and the identification of novel biomarkers. Advances in molecular profiling and model-based research have significantly enhanced the biological understanding of OS and its treatment. Targeted therapies, including antibody-drug conjugates, autologous cell therapies, and immune checkpoint inhibitors, are emerging as promising treatment options. ⁶ Therefore, identifying reliable prognostic biomarkers is crucial for early disease progression prediction and for optimizing personalized treatment strategies.

Cancer immunotherapy represents a strategy that harnesses and enhances the body’s immune system to combat tumors.^7-9 Although significant strides have been made in cancer immunotherapy, challenges persist. Studies have demonstrated that conventional type 1 dendritic cells (cDC1) form clusters with CD8+ T-cells within tumors, promoting the activation and expansion of TCF1+ stem-like CD8+ T-cells, which enhances anti-cancer immunity. ¹⁰ Furthermore, combining lipid nanoparticle-mRNA formulations with dendritic cell therapy (CATCH) has proven effective in overcoming immune suppression, activating dendritic cells, and enhancing specific T-cell immunity, which in turn inhibits tumor metastasis and recurrence. ¹¹ Additionally, cancer metabolism plays a critical role in the growth and survival of cancer cells. Metabolic reprogramming, characterized by enhanced macromolecule synthesis, altered energy metabolism, and the maintenance of redox balance, supports the rapid growth of cancer cells. Ubiquitination and deubiquitination also play essential roles in regulating this metabolic reprogramming. ¹² Although progress in targeting cancer metabolism has been slow, overcoming challenges related to metabolic vulnerabilities in both cancerous and non-cancerous cells within the tumor microenvironment remains critical for designing effective therapeutic strategies. ¹³ Understanding the roles of cancer immunity and metabolism in OS is particularly important, as it provides the potential to develop more effective treatment approaches. Regulating immune responses and metabolic reprogramming could offer new therapeutic avenues, improving the prognosis of OS. Consequently, combining immunotherapy with metabolic-targeted therapy holds promise for more comprehensive and effective treatment of OS.

Machine learning (ML) has proven to be a powerful tool in the study of prognostic biomarkers, offering the ability to extract valuable information from complex, multi-dimensional data and identify potential biological patterns and associations that can guide disease diagnosis and treatment.^14-16 The ML models can effectively screen for biomarkers related to prognosis, enhancing the accuracy and reliability of predictions. The use of multiple models further ensures the robustness of the findings, reducing biases and uncertainties associated with relying on a single method.

In this study, we employed 10 commonly used ML algorithms as complementary tools for feature evaluation and robustness assessment, including Random Forest, Support Vector Machine (SVM, radial basis kernel), XGBoost (xgbDART), Generalized Linear Model (GLM), Naive Bayes, Decision Tree (rpart/CART), K-Nearest Neighbors (KNNs), Logistic Regression, Least Absolute Shrinkage and Selection Operator (LASSO)/Elastic Net (glmnet), and Artificial Neural Network (ANN, multilayer perceptron).^17-21 These algorithms were not used as primary survival models, but rather as auxiliary methods to evaluate the stability and recurrence of candidate prognostic genes across multiple analytical frameworks. By utilizing these diverse models, we aimed to comprehensively assess the robustness of candidate gene selection in OS.

Our goal is to screen and evaluate prognostic biomarkers for OS by applying a broad range of ML models. Through an in-depth evaluation of these models’ predictive performance on OS data, we aim to uncover valuable biological patterns and associations, providing crucial support for disease diagnosis and personalized treatment. Ultimately, our research seeks to offer novel insights and methodologies for optimizing prognostic strategies and improving the clinical outcomes for OS patients.

The rest of this article is structured as follows. Section “Literature Review” reviews the literature. Section “Materials and Methods” describes the methodology. Section “Results” presents the results. Section “Discussion” discusses the findings. Section “Conclusion” summarizes the study.

Literature Review

Osteosarcoma is the most common primary malignant bone tumor in children and adolescents, characterized by high heterogeneity and high rates of metastasis and recurrence. Traditional diagnosis and treatment models have limitations such as strong subjectivity and limited accuracy. ²² In recent years, ML has achieved breakthrough progress in OS detection, classification, and prognostic prediction by virtue of its advantage in processing high-dimensional biomedical data. ²³ The following focuses on the latest research achievements and trends from the end of 2024 to 2025.

The construction of prognostic models is a core direction. Recent studies mostly integrate multidimensional data to enhance model performance: a 2025 study constructed a long-term survival prediction model via the LASSO algorithm and multivariate Cox regression based on macrophage polarization-related genes, and validated the function of the target gene BNIP3; ²⁴ a 2024 study published in Frontiers in Immunology focused on programmed cell death (PCD)-related genes, constructing an integrated prognostic model (Osteosarcoma Programmed Cell Death Score [OS-PCDS]) using the random survival forest (RSF) algorithm with a C-index of 0.943, forming a complete evidence chain; ²⁵ and a 2025 multicenter study combined conventional machine learning with deep learning to build a distant metastasis prediction model and develop an online nomogram. These models provide a scientific basis for personalized treatment. ²⁶ In recent years, the application of artificial intelligence in oncology has achieved breakthrough progress, yet substantial challenges remain regarding the acceptance of its practical clinical application. ²⁷

Accurate classification and early detection are crucial: A 2025 study showed that models such as convolutional neural networks (CNNs) and SVMs, which integrate X-ray imaging and clinical information, can improve the accuracy of distinguishing benign and malignant bone tumors, providing a new scheme for primary screening. ²⁸ In addition, ML is also applied to circulating tumor marker screening and metastatic lesion identification, facilitating early diagnosis and disease monitoring.^29,30

Existing studies have confirmed the core value of ML in OS research, but ML-based prognostic studies targeting metabolism-related genes (pyrroline-5-carboxylate reductase 1 [PYCR1], glycogen phosphorylase, muscle-associated [PYGM]) remain scarce. The roles of both in tumor progression have initially been reported. This study will identify and verify their feasibility as prognostic biomarkers through ML, filling the gaps in existing research and providing new evidence for the prognostic evaluation and therapeutic target development of OS.

Materials and Methods

Data Collection and Preprocessing

Gene expression data for OS were obtained from the UCSC Xena platform, which integrates multiple large-scale cancer genomics data sets, including the TARGET OS cohort. A total of 88 OS samples with available clinical information were included in the training cohort. Normal tissue expression data (n = 396) were obtained from the Genotype-Tissue Expression (GTEx) database. Because data sets derived from different repositories and sequencing platforms may introduce systematic bias, expression matrices were standardized prior to downstream analyses. Specifically, RNA-seq expression values were converted to TPM (Transcripts Per Million) and subsequently transformed using log2(TPM + 1) to reduce distributional bias. To minimize batch effects between data sets, the ComBat algorithm implemented in the “sva” R package was applied to correct batch effects between OS samples obtained from UCSC Xena and normal tissues from the GTEx database. For external validation, 3 independent OS data sets (GSE16091, GSE21257, and GSE146649) were downloaded from the Gene Expression Omnibus (GEO) database. For GEO microarray data sets, probe identifiers were mapped to gene symbols according to the corresponding platform annotation files. When multiple probes corresponded to the same gene, the average expression value was used as the gene expression level. The expression matrices were subsequently normalized using z-score scaling to ensure comparability across data sets. Batch correction was performed independently within each cohort. The GEO data sets were used solely for external validation and were not merged with the training data set to avoid potential data leakage. It should be noted that the normal reference data were obtained from the GTEx database and may not represent perfectly matched normal bone tissue from the same anatomical context as OS specimens. Therefore, comparisons between OS and normal tissues at the bulk transcriptomic level should be interpreted with caution, particularly for metabolism-related genes such as PYGM that may be influenced by tissue composition.

Data Processing and Gene Function Analysis

To ensure data quality, we first integrated the expression data of immune-related and metabolism-related genes. The immune-related gene set was obtained from the Gene Set Enrichment Analysis (GSEA)/Molecular Signatures Database (MSigDB) data set resource (Broad Institute), whereas the metabolism-related gene set was collected from predefined metabolism-associated gene panels curated from publicly available pathway databases and previous literature. Lowly expressed genes were removed to minimize background noise, and genes with mean expression values ⩽0.5 were excluded prior to differential analysis. Differential expression analysis between OS and normal samples was performed on a gene-by-gene basis using the Wilcoxon rank-sum test. P-values were further adjusted for multiple testing using the Benjamini-Hochberg method, and genes with |logFC| > 1 and false discovery rate (FDR) <0.05 were considered significantly differentially expressed. The distribution and statistical significance of all differentially expressed genes were visualized using heatmaps and volcano plots, providing an intuitive overview of global expression changes. The differentially expressed genes identified were then subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses to further elucidate the underlying biological processes and signaling pathways. Enrichment results with P-value < .05 and adjusted P-value (adjP) < .05 were considered statistically significant, ensuring the robustness of biological interpretation. Finally, GO and KEGG enrichment outcomes were displayed as circular and cluster diagrams, respectively, while radar plots were used to visually illustrate the distribution patterns of differentially expressed genes across various functional categories, thereby laying the foundation for subsequent mechanistic studies and functional validation.

Gene Screening and Subtype Classification

We performed univariate Cox regression analysis on the gene expression data to identify significant prognostic genes. These genes were further classified into molecular subtypes using non-negative matrix factorization (NMF).^31,32 The candidate number of clusters was evaluated across k = 2 to k = 4, and the final subtype number was selected by jointly considering clustering stability, consensus matrix structure, and biological interpretability. Based on these criteria, k = 2 was chosen for the downstream subtype analysis. The NMF analysis generated a consensus matrix heatmap, and subtype classifications were determined. A subtype heatmap displayed gene expression patterns across different subtypes. The Kaplan-Meier survival analysis was then conducted to evaluate survival differences among subtypes, visualized through survival curves. These steps helped identify genes significantly associated with prognosis and their corresponding tumor subtypes. Differential expression effect sizes were reported as log2 fold change (log2FC), and group comparisons were performed using the Wilcoxon rank-sum test for transcriptomic data. For survival modeling, hazard ratios (HRs), regression coefficients, and 95% confidence intervals (CIs) were calculated using the Cox proportional hazards regression.

Tumor Microenvironment and Immune Cell Infiltration Analysis

To analyze the tumor microenvironment and immune cell infiltration, we integrated and processed gene expression data, retaining only tumor samples after filtering out low-expression genes. The ESTIMATE package was used to calculate tumor microenvironment scores, visualized through heatmaps. Differences in tumor microenvironment among subtypes were explored using violin plots. The MCPcounter package was utilized to assess immune cell infiltration in tumor samples, calculating infiltration scores. Boxplots were created to compare immune cell infiltration differences among subtypes, combining tumor microenvironment scores, immune cell infiltration scores, and subtype information. A comprehensive heatmap was generated to display the tumor microenvironment and immune cell infiltration characteristics across subtypes.

Univariate Cox Screening, Least Absolute Shrinkage and Selection Operator Selection, Multivariate Cox Modeling, and Machine Learning–Assisted Gene Evaluation

Gene expression and survival data were obtained from the UCSC Xena and GEO databases. After preprocessing, univariate Cox proportional hazards regression was first performed in the training cohort to identify genes significantly associated with overall survival. For each candidate gene, the HR, 95% CI, and P-value were calculated and are provided in Supplemental Table S1. Genes meeting the predefined screening criterion were subsequently entered into the LASSO regression for dimensionality reduction. Based on the retained prognostic genes, a multivariate Cox proportional hazards model was then constructed to calculate the risk score for each patient. The regression coefficient, HR, 95% CI, and P-value for each variable in the final multivariate model are summarized in Supplemental Table S2. Survival prediction in this study was therefore primarily based on Cox regression modeling rather than on conventional classification algorithms.

To further improve the robustness of key gene identification, 10 commonly used ML algorithms (Random Forest, SVM, XGBoost, GLM, Naive Bayes, Decision Tree, KNNs, Logistic Regression, LASSO, and ANN) were applied as complementary feature evaluation tools. These algorithms were used to rank the relative importance of candidate genes across different methods, rather than to directly model time-to-event outcomes. Genes that were repeatedly identified as important across multiple algorithms were retained as core prognostic features.

All prognostic gene screening procedures, including univariate Cox regression, LASSO selection, and ML-assisted feature evaluation, were conducted in the training cohort only. The GEO cohorts were used exclusively for external validation of the final Cox-based prognostic model and were not involved in feature selection or model development, in order to minimize the risk of data leakage.

Survival Analysis, Risk Assessment, and Visualization

We employed comprehensive statistical methods utilizing R packages to conduct survival analysis, risk assessment, and data visualization. Initially, we extracted risk data from UCSC XENA and GEO databases. The Kaplan-Meier survival curves were used to compare survival differences between high- and low-risk groups, with significance evaluated by calculating P-values. Survival curves were then plotted to illustrate the 95% CIs and risk tables. To visualize patient data stratified by risk scores, risk curves and survival status plots were generated, and heatmaps were employed to display gene expression patterns across different risk groups. The predictive performance of the model was assessed using receiver operating characteristic (ROC) curves, plotting 1-year, 3-year, and 5-year ROC curves for both UCSC XENA and GEO data sets, with area under the curve (AUC) values computed accordingly. Additionally, we visualized the coefficients and directions of significant genes from the multivariate Cox regression model using bar plots. These steps provided a robust framework for evaluating and visualizing survival prediction performance and gene expression characteristics across diverse data sets and models. Time-dependent ROC curves were calculated using censored survival data at 1, 3, and 5 years, and the corresponding AUC values were used to assess discrimination of the Cox-based prognostic model over time. In this study, the cutoff for risk stratification was defined a priori as the median risk score in the training cohort. Patients with risk scores greater than or equal to the median were assigned to the high-risk group, whereas those with scores below the median were assigned to the low-risk group. The same fixed cutoff was subsequently applied to the testing cohort and external validation cohorts to ensure consistency of risk classification and to avoid data-driven optimization of group thresholds in the validation setting.

Nomogram Construction, Calibration Curve Plotting, and Immune Cell Infiltration Analysis

We began by retrieving and merging risk scores and clinical data from the UCSC XENA database to construct a nomogram using the Cox regression analysis, aiming to predict 1-year, 3-year, and 5-year survival rates. Calibration curves were generated to evaluate the nomogram’s predictive accuracy, with C-index values calculated to measure model performance. Subsequently, immune cell infiltration results were integrated with risk scores, allowing for a comparative analysis of immune cell infiltration differences between high- and low-risk groups. Violin plots were utilized to depict the distribution of immune cells within different risk groups, and correlation matrices and plots were generated to elucidate the relationships and significance among various immune cells. Calibration performance of the nomogram was assessed by comparing predicted and observed survival probabilities at predefined time points. Because of the limited sample size, calibration results were interpreted as exploratory rather than definitive. Clinical covariates included in the multivariable Cox analysis comprised age, sex, metastasis status, and the Cox-derived risk score. Only cases with complete information for these variables were included in the multivariable analysis. Samples with missing clinical data were excluded from this step using a complete-case analysis strategy, and no imputation was performed.

Gene Expression Data Processing and Visualization

Differentially expressed genes were uploaded to the STRING database to obtain interaction data, which was then visualized using radar plots to depict gene expression distribution. By integrating clinical data and risk profiles, we extracted expression levels of immune checkpoint-related genes, presenting differential expression between high- and low-risk groups through box plots. Furthermore, correlation matrices and plots were employed to conduct a thorough analysis of gene expression and relationships.

Quantitative Real-Time Polymerase Chain Reaction

Total RNA was isolated from OS cell lines (HOS and 143B) and the human normal osteoblast cell line (hFOB1.19) using the TRIzol reagent (Invitrogen, Waltham, Massachusetts), following standard protocols. The concentration and purity of RNA were assessed with a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts). Reverse transcription was performed using the PrimeScript RT reagent kit (Takara Bio, Japan) with 1 µg of RNA, according to the manufacturer’s instructions. Quantitative real-time polymerase chain reaction (qRT-PCR) was carried out using SYBR Green Master Mix (Applied Biosystems, Waltham, Massachusetts) on a StepOnePlus Real-Time PCR System (Applied Biosystems), with specific primers designed for PYCR1, PYGM, and the internal control glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The relative expression of target genes was normalized to GAPDH and calculated using the 2^^−ΔΔCt method. All reactions were performed in triplicate, and data were presented as mean ± standard deviation (SD). Statistical analysis was conducted using 1-way analysis of variance (ANOVA), with a threshold of P < .05 indicating statistical significance. The cell lines used in this study are the human normal osteoblast cell line hFOB1.19 (Xiamen Immocell Biotechnology Co, Ltd, China, RRID: CVCL_3708), the OS cell lines HOS (Xiamen Immocell Biotechnology Co, Ltd, RRID: CVCL_0312), and 143B (Xiamen Immocell Biotechnology Co, Ltd, RRID: CVCL_2270).

Hematoxylin and Eosin Staining

The OS tissue samples were fixed in 10% neutral formalin, dehydrated through a graded ethanol series, cleared with xylene, and embedded in paraffin to prepare 5 μm thick sections. These sections were deparaffinized with xylene and rehydrated through a graded ethanol series. Hematoxylin staining was performed for 5 minutes, followed by tap water rinsing, differentiation with 1% hydrochloric acid alcohol, and bluing with tap water. The sections were then stained with 0.5% eosin for 1 minute and rinsed with tap water. After staining, sections were dehydrated through graded ethanol, cleared with xylene, and mounted with neutral gum. The morphological features of the tumor tissues were examined and photographed under a light microscope. The tissue specimens used for histological and protein expression validation (hematoxylin and eosin [HE] staining, immunohistochemistry, and immunofluorescence) were obtained from archived paraffin-embedded samples stored in the institutional pathology repository of the First Affiliated Hospital of Guangxi Medical University. These samples were anonymized pathological specimens used solely for experimental validation of protein expression patterns. Because the specimens were de-identified before analysis, detailed patient clinical information such as age, sex, tumor stage, treatment history, and follow-up data were not available for this study. Normal control tissues were obtained from non-tumor orthopedic surgical specimens in the same pathology archive. All samples were used only for qualitative and semi-quantitative assessment of protein expression.

Immunohistochemical Staining

Sections were deparaffinized in xylene, rehydrated through graded ethanol, and subjected to antigen retrieval using citrate buffer (pH 6.0) under high pressure. Following cooling to room temperature, sections were incubated with 3% hydrogen peroxide to block endogenous peroxidase activity, followed by phosphate-buffered saline (PBS) washing. To block nonspecific binding, 10% normal goat serum was applied for 30 minutes at room temperature. Specific primary antibodies (anti-PYCR1 and anti-PYGM) were added and incubated overnight at 4°C. The next day, sections were washed with PBS, incubated with biotinylated secondary antibodies for 30 minutes at room temperature, and washed again with PBS. Horseradish peroxidase (HRP)-labeled streptavidin was applied and incubated for 30 minutes at room temperature, followed by PBS washing. 3,3’-Diaminobenzidine (DAB) was used for chromogenic detection, and the reaction was terminated with tap water. Hematoxylin counterstaining was performed, followed by dehydration through graded ethanol, clearing with xylene, and mounting with neutral gum. Immunohistochemical results were examined and photographed under a light microscope to assess target protein expression. The antibodies used in this study were rabbit polyclonal anti-PYCR1 antibody (Catalog No: 13108-1-AP, Proteintech Group, Wuhan, China; IHC working dilution: 1:2000) and rabbit polyclonal anti-PYGM antibody (Catalog No: biosb-5012R, Beijing, China; IHC working dilution: 1:2000).

Immunohistochemical staining results were evaluated using a semi-quantitative H-score approach based on staining intensity and the proportion of positive cells. Two independent observers who were blinded to the experimental grouping evaluated the staining results. In cases of disagreement, the final score was determined by joint review. For each tissue specimen, at least 2 sections were examined to ensure staining consistency. Negative control slides were processed in parallel without primary antibody incubation to verify staining specificity.

Immunofluorescence Staining

Immunofluorescence staining was performed to evaluate the spatial expression of PYCR1 and PYGM in OS and normal control tissues. After fixation in 4% paraformaldehyde, tissue sections were washed with PBS, permeabilized with 0.1% Triton X-100 for 10 minutes, and blocked with 5% bovine serum albumin (BSA) for 30 minutes at room temperature. Sections were then incubated overnight at 4°C with primary antibodies against PYCR1 and PYGM (working dilution: 1:500). After washing with PBS, the sections were incubated with species-appropriate fluorescent secondary antibodies (working dilution: 1:500) for 1 hour at room temperature in the dark. Nuclei were counterstained with 4”,6-diamidino-2-phenylindole (DAPI), and images were captured using a fluorescence microscope under identical exposure settings for all groups to enable qualitative comparison. All immunofluorescence staining experiments were repeated independently to confirm the reproducibility of the observed expression patterns.

Results

Identification of Differentially Expressed Genes and Molecular Subtypes in Osteosarcoma

Differential expression analysis based on the Wilcoxon rank-sum test and Benjamini-Hochberg multiple-testing correction revealed marked disparities in the expression profiles of immune- and metabolism-related genes between OS and normal tissues. As shown in the heatmap (Figure 1A), unsupervised clustering based on these differentially expressed genes clearly distinguished tumor samples from normal controls. The volcano plot (Figure 1B) further illustrated the overall distribution and directionality of gene expression changes, with a clear segregation of upregulated and downregulated genes. Functional annotation demonstrated that these genes were primarily enriched in biological processes closely associated with the tumor immune microenvironment and metastatic potential, including leukocyte migration, cell chemotaxis, and granulocyte chemotaxis, as visualized by the GO chord diagram (Figure 1C). The KEGG pathway analysis, presented as a circular cluster plot (Figure 1D), showed significant enrichment of these genes in several cancer- and immunity-related pathways, such as the chemokine signaling, JAK-STAT signaling pathway, and Ras signaling pathways. To further explore molecular heterogeneity within OS, unsupervised clustering using NMF was conducted. The NMF rank survey (Figure 1E) indicated optimal clustering stability when the number of clusters ranged from 2 to 4. The consensus clustering heatmap (Figure 1F) visualized the consistency of sample clustering across different subtype assignments. Collectively, these results highlight pronounced molecular heterogeneity and distinct biological features within OS, providing a robust foundation for subsequent investigations into subtype-specific mechanisms.

Figure 1.

Analysis of differentially expressed genes and molecular subtypes in OS. Panel A shows a heatmap clustering of OS and normal tissues based on immune- and metabolism-related differentially expressed genes (OS, n = 88; normal, n = 396), with clear distinction between the 2 groups. Panel B presents a volcano plot illustrating genes with elevated or reduced expression in tumors vs normal tissues. Differential expression was assessed using the Wilcoxon rank-sum test, and multiple-testing correction was performed using the Benjamini-Hochberg method. Genes with |logFC| > 1 and FDR < 0.05 were considered significant. Panel C is a GO enrichment chord diagram, indicating that the DEGs are primarily involved in tumor immunity-related processes such as leukocyte migration and cell chemotaxis. Panel D is a KEGG circular cluster diagram, showing the enrichment of DEGs in classic tumor-associated pathways including chemokine, JAK-STAT, and Ras signaling. Panel E displays the NMF rank survey, determining the optimal number of clusters and validating subtype stability. Panel F shows the NMF consensus clustering heatmap, revealing OS molecular subtypes and their heterogeneity.

Survival Differences and Tumor Microenvironment Characteristics Among Osteosarcoma Molecular Subtypes

Based on the molecular subtypes identified above, we systematically evaluated differences in survival outcomes and tumor microenvironment characteristics between subgroups. The consensus clustering heatmap (Figure 2A) clearly classified OS samples into 2 major subtypes (C1 and C2), with high clustering consistency observed among samples. The Kaplan-Meier survival analysis (Figure 2B) demonstrated that patients with the C2 subtype exhibited a significantly improved overall survival compared to those with the C1 subtype. Comparative analysis of immune cell infiltration between subtypes revealed distinct patterns, as shown in the heatmap (Figure 2C), with pronounced differences in various immune cell populations and tumor microenvironment scores. Quantitative assessment using boxplots (Figure 2D) indicated that the C2 subtype displayed higher infiltration levels of multiple immune cell types, suggesting a more active immune microenvironment. The tumor microenvironment score heatmap (Figure 2E) further detailed differences between subtypes in terms of stromal score, immune score, and the ESTIMATE score. Violin plots (Figure 2F) provided a visual summary, illustrating that C2 had significantly higher stromal, immune, and ESTIMATE scores than C1. Collectively, these findings reveal substantial biological differences between OS molecular subtypes in terms of survival outcomes, immune cell infiltration, and tumor microenvironmental features, thereby providing a molecular basis for risk stratification and the development of targeted therapeutic strategies.

Figure 2.

Survival and tumor immune microenvironment analysis of OS molecular subtypes. Panel A displays a consensus clustering heatmap, identifying 2 molecular subtypes of OS samples. Panel B shows the Kaplan-Meier survival curves, with the C2 subtype exhibiting a significantly better survival rate than C1. Panel C presents a heatmap of immune cell infiltration and tumor microenvironment scores, highlighting substantial immune differences between subtypes. Panel D is a boxplot of immune cell infiltration, showing generally higher infiltration in the C2 subtype. Panels E and F (heatmap and violin plot, respectively) demonstrate that the C2 subtype has significantly higher stromal, immune, and ESTIMATE scores compared to C1.

Multi-dimensional Selection of Prognostic Risk Genes and Identification of Key Features

Based on the initial screening, the LASSO regression was applied to further refine the set of genes associated with prognosis in OS. The LASSO coefficient trajectories and cross-validation curves (Figure 3A and B) supported parsimonious feature selection and minimized overfitting, ultimately retaining a subset of candidate prognostic genes. To reduce potential bias arising from reliance on any single method, 10 ML algorithms were subsequently used as complementary approaches to evaluate feature importance among these candidate genes. Rather than directly modeling survival endpoints, these algorithms were used to assess the relative contribution and recurrence of candidate genes across multiple analytical frameworks. Feature importance ranking (Figure 4B) showed that, although each algorithm assigned different weights to individual genes, PYCR1 and PYGM were consistently retained across multiple models. The UpSet plot (Figure 4E) further demonstrated that PYCR1 and PYGM were recurrently identified by all 10 algorithms, supporting their robustness as core prognostic features in OS.

Figure 3.

LASSO regression and prognostic risk model construction and evaluation. Panel A displays the LASSO coefficient trajectory, and panel B shows the cross-validation for optimal λ selection. Panels C and D show the Kaplan-Meier survival curves in the training and validation cohorts, confirming the model’s ability to distinguish prognosis between high- and low-risk patients. Panel E presents the regression coefficients of the multi-gene model. Panels F and G display the ROC curves at various time points for the training and validation sets, indicating strong prognostic predictive ability.

Figure 4.

Machine learning–assisted feature evaluation of prognostic candidate genes. Panel A shows the residual distributions across different algorithms. Panel B presents the feature importance ranking generated by each method. Panel C shows the model classification summaries used for auxiliary comparison among algorithms. Panel D presents the cumulative residual distributions. Panel E is an UpSet plot highlighting PYCR1 and PYGM as recurrent candidate genes identified across all 10 algorithms. Panels F to N show ROC-based auxiliary performance summaries of the algorithms in the context of feature evaluation. These analyses were used only for complementary assessment of feature stability rather than direct survival prediction.

Construction and Evaluation of the Prognostic Risk Model

Based on the core genes identified from the intersection of multiple ML models, a multi-gene prognostic risk score model was constructed and systematically evaluated in both the training and independent validation cohorts. The Kaplan-Meier survival analysis (Figure 3C and D) demonstrated that this model effectively stratified patients into high- and low-risk groups, with the high-risk group exhibiting significantly worse overall survival; the difference between groups was statistically significant. Visualization of regression coefficients (Figure 3E) illustrated the specific contribution and direction of each gene included in the risk score. Further ROC curve analysis showed that the AUC values of the PYCR1/PYGM-related risk score model at 1 year, 3 years, and 5 years were 0.897, 0.917, and 0.872 in the UCSC XENA training cohort, and 0.968, 0.868, and 0.849 in the GEO external validation cohort, respectively (Figure 3F and G), reflecting strong predictive accuracy and robustness in both the training and validation sets. These findings collectively support the potential clinical utility and prognostic value of this multi-gene model for risk assessment in OS patients. The univariable Cox regression results for candidate genes and the multivariable Cox regression results integrating the final prognostic genes are summarized in Supplemental Tables S1 and S2, respectively, including HRs, 95% CIs, and P-values.

Prognostic Stratification and Clinical Tool Construction of the Risk Model in Different Cohorts

In the UCSC XENA cohort (Figure 5A1 to A3), the risk model demonstrated robust stratification and prognostic prediction capability. Specifically, panel A1 shows that patients in the high-risk group exhibited shorter survival times, with the majority of death events clustered in this group, indicating the model’s effectiveness in identifying individuals with poor prognosis. The risk score distribution in A2 clearly delineated high- and low-risk groups, with sharp increases in risk score corresponding to a rapid decline in survival outcomes. The heatmap in A3 revealed distinct expression patterns of model-included genes between high- and low-risk groups, supporting the molecular basis of risk stratification. Overall, the model’s performance in the UCSC XENA cohort confirmed its ability to discriminate patient survival risk.

Figure 5.

Prognostic stratification and nomogram construction based on the risk model. Panels A1 to A3 display the risk distribution, survival status, and a heatmap of key gene expression in the UCSC XENA cohort, all indicating lower survival and distinctive gene expression patterns in high-risk patients. Panels B1 to B3 show the consistent results in the external GEO validation cohort. Panel C presents a calibration curve, and panel D shows a nomogram tool for individualized survival probability prediction.

In the GEO validation cohort (Figure 5B1 to B3), the model also demonstrated high robustness and generalizability. B1 indicated that high-risk patients in this independent data set likewise experienced shorter survival, with death events closely associated with higher risk scores. The risk score curve in B2 showed similarly clear stratification in the validation cohort, reflecting the model’s reliable discriminative power. The heatmap in B3 also displayed marked differences in core gene expression between risk groups, further supporting the model’s reproducibility and stability across cohorts. Together, these findings indicate that the risk model enables effective prognostic stratification not only in the training set but also in external, independent data sets.

Furthermore, the calibration curve (Figure 5C) showed an overall agreement between predicted and observed survival probabilities. However, given the limited sample size and the potential optimism of model fitting, the nomogram performance should be interpreted cautiously and requires further validation in larger independent cohorts. Finally, Figure 5D presents a nomogram that integrates risk score with clinical features such as metastasis status, age, and gender, enabling individualized risk assessment and conversion of total risk points into survival probability, thus supporting precise prognosis prediction and clinical decision-making. The corresponding multivariable Cox regression statistics for the clinical covariates included in the nomogram are also provided in Supplemental Table S2. Collectively, the multi-gene risk model and its clinical tool exhibited outstanding prognostic stratification ability and significant potential for clinical translation.

Expression Patterns and Functional Validation of Key Risk Genes Pyrroline-5-Carboxylate Reductase 1 and Glycogen Phosphorylase, Muscle-Associated

This study further focused on PYCR1 and PYGM, 2 key risk genes consistently identified across multiple statistical and ML models, and examined their expression patterns at different analytical levels. Radar plots (Figure 6A1 and A2) illustrated the overall expression profiles of PYCR1 and PYGM within the candidate gene set, showing that PYCR1 was generally highly expressed, whereas PYGM displayed relatively lower bulk expression. Within the OS cohort, subgroup analysis revealed that both PYCR1 (Figure 6B1) and PYGM (Figure 6B2) were significantly more highly expressed in the high-risk group than in the low-risk group, indicating that relatively elevated expression of these genes within tumor samples was associated with unfavorable prognosis. At the cellular level, qRT-PCR results (Figure 6C1 and C2) further showed that PYCR1 expression was markedly higher in OS cell lines (HOS and 143B) than in the normal osteoblast cell line (hFOB1.19), whereas PYGM expression was significantly lower in tumor cells than in normal osteoblasts. These findings suggest that PYGM exhibits context-dependent expression patterns; although its expression is relatively higher in high-risk OS samples than in low-risk samples, it is reduced in tumor cells compared with normal osteoblasts. One possible explanation is that bulk tumor tissues contain not only malignant cells but also stromal, vascular, and immune components, which may contribute to the overall PYGM signal. Collectively, PYCR1 showed a relatively consistent tumor-promoting expression pattern across data sets and experiments, whereas PYGM demonstrated analysis-level-dependent variation, underscoring the biological complexity of metabolism-related markers in OS.

Figure 6.

Expression grouping and cellular validation of key genes PYCR1 and PYGM. Panels A1 and A2 show radar plots; panels B1 and B2 demonstrate significantly increased expression of both genes in the high-risk group. Panels C1 and C2 display qRT-PCR results, further confirming the upregulation of PYCR1 and downregulation of PYGM in tumor cell lines. Statistical comparisons were performed using 1-way ANOVA for qRT-PCR experiments, and P < .05 was considered statistically significant.

Validation of Pyrroline-5-Carboxylate Reductase 1 and Glycogen Phosphorylase, Muscle-Associated Expression by Histopathology and Immunohistochemistry

The HE staining (Figure 7A1, A2, B1, and B2) revealed that OS tissues (A1, A2), under both low (100×) and high (400×) magnification, displayed marked cellular atypia, including enlarged and irregularly distributed nuclei, pale cytoplasmic staining, and abundant vascularization and inflammatory cell infiltration in the stroma, all consistent with the classic histopathological features of OS. In contrast, the normal control group (B1, B2) exhibited intact bone structure, regular lamellar arrangement, and no evidence of abnormal cells, indicating normal histological architecture. Immunohistochemical analysis further compared the protein expression patterns of PYCR1 and PYGM between OS and normal tissues. The PYCR1 expression was markedly higher in OS tissues (C1) than in normal controls (C2), with positive signals primarily localized to the cytoplasm and, to a lesser extent, the nuclei of tumor cells. For PYGM, immunoreactivity was detectable in OS tissues (D1) and, in some specimens, appeared stronger than that in normal control tissues (D2). However, given the heterogeneity of tissue composition and the discrepancy between bulk tissue-, risk group-, and cell line–based analyses, the PYGM protein-level findings should be interpreted with caution. Rather than indicating a uniformly increased expression in all OS contexts, these data suggest that PYGM may show context-dependent protein distribution in clinical tissue specimens.

Figure 7.

Pathological and immunohistochemical validation of PYCR1 and PYGM expression. Panels A1 and A2 are HE-stained images of OS; panels B1 and B2 are HE-stained normal tissues, illustrating structural differences. Panels C1 and D1 show the representative immunohistochemical staining of PYCR1 and PYGM in OS tissues, whereas C2 and D2 show the staining in normal control tissues. The PYCR1 displayed stronger staining in OS, while PYGM showed detectable but heterogeneous staining across specimens.

Immunofluorescence-Based Expression of Pyrroline-5-Carboxylate Reductase 1 and Glycogen Phosphorylase, Muscle-Associated in Osteosarcoma Tissues

To further evaluate the spatial distribution of PYCR1 and PYGM proteins, immunofluorescence staining was performed in OS and normal control tissues. The DAPI was used for nuclear staining (blue), while PYCR1 and PYGM were labeled with red fluorescence (Figure 8). In normal tissues, PYCR1 fluorescence was weak, whereas OS tissues showed markedly increased PYCR1 signal, consistent with its upregulated expression pattern. For PYGM, fluorescence signals were detectable in both OS and normal tissues, with variable intensity across specimens. In some OS tissues, PYGM staining appeared relatively strong, suggesting that PYGM protein distribution in clinical samples may not fully mirror the expression pattern observed in tumor cell lines. This discrepancy may reflect the influence of non-tumor cellular components and tissue heterogeneity in bulk specimens. Overall, the immunofluorescence results support the presence and cytoplasmic localization of PYCR1 and PYGM in OS tissues, while also indicating that PYGM should be interpreted as a context-dependent marker rather than a uniformly upregulated gene in OS.

Figure 8.

Immunofluorescence analysis of PYCR1 and PYGM expression. Panels A1 to A3 and C1 to C3 are normal tissues; B1 to B3 and D1 to D3 are OS tissues. The DAPI marks nuclei, and red fluorescence marks PYCR1/PYGM. Results show that PYCR1 is markedly upregulated in OS tissues and predominantly localized in the cytoplasm, whereas PYGM is detectable in OS tissues with heterogeneous staining intensity, suggesting context-dependent protein expression.

Discussion

In this study, we performed a comprehensive analysis of OS gene expression data using integrated statistical survival modeling and ML-assisted feature evaluation, ultimately identifying PYCR1 and PYGM as key prognostic biomarkers. We collected gene expression data from both OS and normal samples obtained from the UCSC XENA and GEO databases. Initially, the univariate Cox regression analysis was used to screen for significant prognostic genes, which were then refined through LASSO regression. To further elucidate the tumor’s characteristics, we employed the NMF method to classify the tumors into 2 primary subtypes, C1 and C2. The Kaplan-Meier survival analysis revealed significant survival differences between these subtypes. The ROC curve analysis showed that the 1-year, 3-year, and 5-year AUC values of the PYCR1/PYGM-related risk scoring model were 0.897, 0.917, and 0.872 in the UCSC XENA training cohort, and 0.968, 0.868, and 0.849 in the GEO external validation cohort, confirming the stability of the core gene screening results. Moreover, Tumor Microenvironment and immune cell infiltration analyses indicated notable variations between high- and low-risk groups. To validate these findings, HE staining, immunohistochemistry, and immunofluorescence were employed, confirming the robust upregulation of PYCR1 in OS tissues and the detectable but context-dependent expression pattern of PYGM in OS tissues and highlighting their specific localization within tumor cells. These results provide a systematic evaluation of the biological significance and clinical relevance of PYCR1 and PYGM, underscoring their potential as prognostic biomarkers in OS. In this study, we further enhanced the reporting transparency of the prognostic model by presenting full results of univariate and multivariate Cox regression analyses, clearly defining the cutoff value for risk stratification, and specifying the processing methods for missing clinical covariates.

Our findings indicate that PYCR1 and PYGM hold significant prognostic value in OS. Previous studies have demonstrated the pivotal role of PYCR1 in various cancers. In breast cancer, PYCR1 has been shown to promote collagen synthesis in cancer-associated fibroblasts (CAFs), which in turn drives tumor invasiveness. Reducing PYCR1 levels in CAFs results in decreased tumor collagen production, tumor growth, and metastatic spread. ³³ Under hypoxic conditions, PYCR1 is phosphorylated by nuclear IGF1R, which promotes its binding to ELK4, thereby inhibiting gene transcription and sustaining cell growth. ³⁴ In bladder cancer, PYCR1 enhances cell proliferation and invasion. Delivery of si-PYCR1 into bone marrow mesenchymal stem cell (BMSC)-derived exosomes significantly reduces malignancy and aerobic glycolysis in bladder cancer cells by inhibiting the epidermal growth factor receptor (EGFR)/phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) pathway. ³⁵ Furthermore, the dysregulation of N6-methyladenosine (m6A) modification, a process involving the m6A demethylase fat mass and obesity-associated protein (FTO), has been implicated in bladder cancer development by stabilizing PYCR1 transcripts. ³⁶ These studies not only emphasize the importance of PYCR1 at the gene expression level but also highlight its involvement in post-transcriptional modifications, providing insight into the multi-layered regulatory mechanisms of PYCR1 in OS.

In our study, PYGM showed a complex and context-dependent expression pattern rather than a unidirectional change. At the bulk transcriptomic level and in qRT-PCR assays comparing OS cell lines with normal osteoblasts, PYGM tended to be downregulated, suggesting reduced tumor cell–intrinsic expression. However, within the OS cohort, relatively higher PYGM expression was observed in the high-risk group, and protein-based assays detected PYGM signals in clinical tissue specimens. These findings may be explained by differences in analytical resolution, as bulk tissues include stromal, vascular, inflammatory, and other non-malignant cellular components that are absent in cell line models. In addition, the normal reference tissues derived from GTEx may not be perfectly matched to OS-adjacent normal bone, which may further contribute to the observed discrepancy. Therefore, PYGM is better interpreted as a context-dependent metabolism-related biomarker in OS rather than a uniformly downregulated or upregulated gene. This aligns with other research that underscores the significant roles of PYCR1 and PYGM in various cancers. For instance, glycerol-3-phosphate dehydrogenase (GPI), a key enzyme in glycolysis, has been targeted by oleuropein to inhibit glycolysis, demonstrating significant anti-tumor activity. ³⁷ This observation parallels our finding that PYGM is involved in glycogen metabolism, suggesting that regulation of glucose metabolic pathways may be an effective strategy for inhibiting tumor growth. Additionally, PYGM is implicated in insulin and glucagon signaling pathways, insulin resistance, necroptosis, immune response, and phototransduction. ³⁸ This multifaceted role of PYGM in OS further underscores its potential as a therapeutic target. Furthermore, recessive mutations in the PYGM gene have been linked to McArdle disease, which leads to muscle glycogen phosphorylase deficiency, causing exercise intolerance, muscle weakness, and acute rhabdomyolysis. ³⁹ These findings highlight the critical role of PYGM in muscle and other tissues, reinforcing its relevance to our OS research. Additionally, studies on head and neck squamous cell carcinoma (HNSCC) have revealed that PYGM and TNNC2 are significantly downregulated in HNSCC, with this aberrant expression correlating with prognosis. ⁴⁰ This suggests that PYGM could hold prognostic value not only in OS but also in other cancer types, positioning it as a potential biomarker or therapeutic target.

We employed multiple ML models to provide a more comprehensive and systematic approach for biomarker screening. The PYCR1 and PYGM were identified as prognostic biomarkers for OS for the first time, thanks to the robust and integrative methods used. By comparing our multi-model approach to studies relying on single methods, we have reduced potential biases and increased the reliability of our screening results. As prognostic biomarkers, PYCR1 and PYGM hold significant clinical potential. They can be utilized for early diagnosis, prognostic prediction, and guiding personalized treatment strategies for OS patients. The expression levels of these genes can help identify high-risk patients, enabling the development of more tailored and aggressive treatment plans. Furthermore, these biomarkers may serve as potential targets for targeted therapies, offering a foundation for the development of novel therapeutic approaches.

However, our study does have some limitations. First, the relatively small sample size may affect the generalizability of our findings. Second, despite employing various ML models, the selection of models and parameter settings may still introduce bias. Additionally, our study primarily relies on gene expression data and lacks validation at the protein level or through functional experiments. Future research should expand the sample size, incorporate more diverse data sources, and explore functional validation to further refine and validate our findings.

Conclusion

This study identified and experimentally validated PYCR1 and PYGM as potential prognostic biomarkers for OS using integrated Cox-based survival analysis and ML-assisted feature screening.

Supplemental Material

sj-docx-1-onc-10.1177_11795549261452567 – Supplemental material for Machine Learning–Based Identification and Validation of PYCR1 and PYGM as Prognostic Biomarkers for Osteosarcoma

Supplemental material, sj-docx-1-onc-10.1177_11795549261452567 for Machine Learning–Based Identification and Validation of PYCR1 and PYGM as Prognostic Biomarkers for Osteosarcoma by Guoyong Xu, Chong Liu, Jiang Xue, Jiarui Chen, Zhuan Zou, Sen Mo, Zhongxian Zhou and Xinli Zhan in Clinical Medicine Insights: Oncology

Supplemental Material

sj-docx-2-onc-10.1177_11795549261452567 – Supplemental material for Machine Learning–Based Identification and Validation of PYCR1 and PYGM as Prognostic Biomarkers for Osteosarcoma

Supplemental material, sj-docx-2-onc-10.1177_11795549261452567 for Machine Learning–Based Identification and Validation of PYCR1 and PYGM as Prognostic Biomarkers for Osteosarcoma by Guoyong Xu, Chong Liu, Jiang Xue, Jiarui Chen, Zhuan Zou, Sen Mo, Zhongxian Zhou and Xinli Zhan in Clinical Medicine Insights: Oncology