Osteosarcoma (OS) is the most common primary malignant bone tumor in orthopedic oncology, and it strikes disproportionately during childhood and adolescence. The prognosis for localized disease has improved with modern chemotherapy regimens, but patients who develop metastases face a 5-year survival rate below 20%. Despite decades of research, reliable early diagnostic biomarkers for OS remain limited, and the molecular mechanisms that drive its aggressive behavior, particularly in the context of immune dysregulation and the tumor microenvironment, are incompletely understood.
A new form of cell death: This study is framed around cuproptosis, a recently characterized copper-dependent cell death mechanism that differs fundamentally from apoptosis, necroptosis, and ferroptosis. Copper ions at concentrations exceeding the cell's homeostatic threshold bind directly to lipid-acylated components of the tricarboxylic acid (TCA) cycle, causing pathological aggregation of acylated proteins and loss of iron-sulfur cluster proteins. A defined set of 10 genes has been identified as cuproptosis regulators, including CDKN2A, FDX1, DLD, DLAT, LIAS, GLS, LIPT1, MTF1, PDHA1, and PDHB. Whether and how these genes are dysregulated in osteosarcoma was, at the time of this study, essentially unknown.
The immune context: Immune dysregulation is increasingly recognized as a central feature of OS. Osteosarcoma cells, immune cells, osteoblasts, and osteoclasts interact within the tumor microenvironment, and immune checkpoint inhibition has emerged as a promising therapeutic avenue. This study positions cuproptosis-related genes not only as potential diagnostic markers but also as potential links to the immune infiltration patterns that shape OS behavior. The authors explicitly aim to identify biomarkers that can guide both early diagnosis and immunotherapy selection.
The study's primary objective is to use three parallel artificial intelligence methods, along with standard bioinformatic screening tools, to converge on the most specific and accurate cuproptosis-related diagnostic biomarkers for osteosarcoma from genome-wide expression data, then validate those findings at the protein level (immunohistochemistry), the transcript level (RT-qPCR), and in a real-world clinical blood test dataset of 25,665 cases.
The study draws from three independent public genomic repositories to build its training and validation sets. Genome-wide mRNA expression data for the training set were downloaded from the UCSC Xena database, comprising 88 osteosarcoma samples, and paired with 396 normal skeletal muscle samples from the GTEx (Genotype-Tissue Expression) database. These two datasets were normalized and processed to remove inter-batch differences before being merged. The validation set was sourced separately from the GEO (Gene Expression Omnibus) database, specifically the GSE42352 dataset, which contains 84 osteosarcoma samples and 3 osteoblast controls (87 samples total). Probe numbers from both training and validation platforms were converted into standard gene symbols using Strawberry Perl scripts.
Differential expression analysis: The limma package in R (version 4.1.3) was used for genome-wide differential expression analysis with cutoff values of |logFC| greater than or equal to 1.5 and adjusted p-value below 0.05. From 54,751 total genes analyzed, this yielded 4,811 differentially expressed genes. The top 100 were visualized as a heatmap and the full set as a volcano plot. The expression of all 10 cuproptosis-related genes was extracted from the expression matrix and cross-correlated using corrplot and circlize packages.
Random forest: The randomForest algorithm is a supervised ensemble learning method that generates multiple decision trees simultaneously and aggregates their outputs to improve accuracy and reduce overfitting. Among its advantages for this application are its ability to handle large numbers of input features, its built-in variable importance ranking, and its tolerance for mixed data types. The cuproptosis genes were ranked by importance score, and the top-ranked genes were carried forward.
SVM-RFE: Support Vector Machine with Recursive Feature Elimination (SVM-RFE) is a wrapper-based feature selection algorithm. SVM classifiers seek the maximum-margin hyperplane that correctly separates two classes; RFE iteratively removes the lowest-ranked features and retrains until the optimal feature subset is found. SVM-RFE was implemented using the e1071, kernlab, and caret packages in R. The analysis found that optimal diagnostic efficacy was achieved when 4 cuproptosis gene variables were retained: PDHA1, PDHB, CDKN2A, and DLD.
LASSO regression: Least Absolute Shrinkage and Selection Operator (LASSO) regression achieves variable selection by penalizing model complexity, compressing coefficients of non-informative features to exactly zero. This produces sparse, interpretable models. The LASSO analysis identified that optimal diagnostic performance was achieved with a 6-variable model. Taking the intersection of five separate gene lists (differential expression, cuproptosis subset, random forest, SVM-RFE, and LASSO), two genes consistently appeared across all methods: PDHA1 and CDKN2A, confirmed as the final diagnostic targets using a Venn diagram.
The random forest analysis produced a ranked importance plot of the 10 cuproptosis-related genes. PDHA1 and CDKN2A occupied the top two positions in this ranking, meaning that among all cuproptosis-associated candidate genes, these two provided the greatest discriminative power for distinguishing osteosarcoma samples from normal controls in the training set. Random forest importance scores reflect how much each gene's inclusion or exclusion degrades classification accuracy across the ensemble of decision trees, making this a robust measure of predictive value.
SVM-RFE performance curve: The SVM-RFE analysis swept across different numbers of retained features and plotted classification accuracy as a function of feature count. The curve identified that the optimal model contained 4 variables: PDHA1, PDHB, CDKN2A, and DLD. Below 4 features, accuracy dropped noticeably, and adding more features beyond 4 produced diminishing or no improvements. This result confirmed that a small, curated feature set is sufficient for high-performance classification and that PDHA1 and CDKN2A are among the minimal necessary variables.
LASSO lambda selection: In LASSO regression, the regularization parameter lambda controls the degree of coefficient shrinkage. As lambda increases, more coefficients are driven to zero and fewer features remain in the model. The cross-validation curve identified that 6 variables produced the best balance of model fit and parsimony. The genes retained after LASSO shrinkage overlapped substantially with those from the other two methods, reinforcing the convergence on PDHA1 and CDKN2A.
The Venn diagram intersection of five independent gene lists (genome-wide differentially expressed genes, cuproptosis-specific differentially expressed genes, random forest top genes, SVM-RFE selected genes, and LASSO retained genes) yielded exactly two genes present in all five categories: PDHA1 and CDKN2A. This multi-method convergence strengthens confidence that these genes are genuinely informative for OS diagnosis rather than artifacts of any single algorithmic approach. The authors describe this convergent AI strategy as a deliberate hedge against method-specific bias.
Receiver Operating Characteristic (ROC) curves were constructed for both genes in both the training set and the independent validation set, providing four total performance assessments. The area under the ROC curve (AUC) is the primary metric, with 1.0 representing perfect discrimination and 0.5 representing chance performance. The training set results were notably strong: PDHA1 achieved an AUC of 100% in the training set, indicating perfect separation of OS from normal controls within the training data. CDKN2A reached a training set AUC of 93.9%.
Validation set performance: External validation on the independent GSE42352 dataset, which the model was not trained on, produced more modest but still clinically meaningful results. PDHA1 achieved a validation set AUC of 88.5%, while CDKN2A achieved 68.7%. The drop from training to validation performance is consistent with the general principle that training metrics overestimate real-world generalizability, and the authors acknowledge this gap. However, they note that AUC values substantially above the 0.6 threshold (their benchmark) in both datasets support the clinical relevance of these genes as diagnostic markers.
Interpreting the asymmetry between the two genes: PDHA1 maintained much stronger diagnostic performance in the validation set (88.5%) compared to CDKN2A (68.7%). This asymmetry may reflect differences in how consistently these genes are dysregulated across different patient populations, tissue sources, and microarray platforms. CDKN2A's validation AUC, while lower, still substantially exceeds chance and may retain additive diagnostic value when combined with PDHA1 in a multi-gene panel. The authors suggest that these two genes should be considered as complementary components of a diagnostic signature rather than standalone markers.
Gene Set Enrichment Analysis (GSEA) was performed on both diagnostic genes to identify which biological processes and molecular pathways are most strongly associated with their expression patterns in osteosarcoma. For PDHA1, Gene Ontology (GO) enrichment was concentrated in chromatin assembly or disassembly, mitochondrial gene expression, and mitochondrial translation. The KEGG pathway enrichment for PDHA1 mapped to the ERBB signaling pathway, pyrimidine metabolism, the ribosome, and the spliceosome. These findings connect PDHA1 to fundamental cellular energy metabolism, since PDHA1 encodes the E1-alpha subunit of the pyruvate dehydrogenase complex, which converts pyruvate to acetyl-CoA at the entry point to the TCA cycle.
PDHA1 in cancer metabolism: Prior research cited in the discussion establishes that PDHA1 activity is frequently downregulated in cancer through phosphorylation, which limits glucose-derived carbon entry into the TCA cycle and redirects metabolism toward aerobic glycolysis (the Warburg effect). PDHA1 has been specifically linked to prostate cancer lipid synthesis through its nuclear localization in tumor cells, and its enzymatic activity is inhibited by phosphorylation in cells with fumarate hydratase (FH) defects. The connection to the ERBB signaling pathway is particularly relevant given that ERBB2/HER2 and related kinases phosphorylate metabolic enzymes to rewire energy metabolism in tumor cells.
CDKN2A pathways: CDKN2A GO enrichment was concentrated in phagocytosis, regulation of translation initiation, plasma membrane signaling receptor complex, and T-cell receptor complexes. Its KEGG pathways included natural killer cell-mediated cytotoxicity, the NOD-like receptor (NLR) signaling pathway, and the T cell receptor signaling pathway. These enrichments place CDKN2A at the interface of cell cycle regulation and immune signaling, consistent with its established role as a tumor suppressor whose loss disrupts the G1/S checkpoint. Its presence in T-cell and NK-cell pathway enrichments connects CDKN2A expression to the immune surveillance mechanisms that osteosarcoma exploits for immune evasion.
Protein-protein interaction (PPI) networks: Both PDHA1 and CDKN2A were analyzed using the STRING database at an elevated confidence cutoff (score 3.2) to build high-confidence interaction networks. Both genes showed connections to multiple hub proteins within their respective networks. Radar charts were constructed to visualize the correlation between each diagnostic gene and its network neighbors, confirming dense connectivity and suggesting that dysregulation of PDHA1 or CDKN2A would ripple through broader regulatory networks rather than operating in isolation.
To characterize immune cell infiltration in osteosarcoma, the study applied CIBERSORT, a computational deconvolution algorithm that estimates the proportional composition of 22 distinct immune cell types from bulk tissue gene expression profiles. CIBERSORT was run on the 88 OS samples and 396 normal controls from the training set. The analysis revealed that the overall immune cell composition of the OS group differed substantially from normal controls, with multiple immune cell populations showing statistically significant differential infiltration.
Key immune cell differences: Among the 22 cell types analyzed, the following showed significant differential expression between OS and normal controls: neutrophils, eosinophils, naive B cells, resting CD4 memory T cells, CD8 T cells, and naive CD4 T cells. Most notably, neutrophils were elevated in the OS group. This finding is biologically significant because tumor-associated neutrophils (TANs) can exhibit both pro-tumorigenic and anti-tumorigenic phenotypes depending on the local cytokine environment, and elevated neutrophil-to-lymphocyte ratio (NLR) is an established adverse prognostic indicator in multiple cancers including bone sarcomas.
Gene-immune cell correlations: Scatter plots were generated to examine correlations between PDHA1 and CDKN2A expression and each of the six differentially infiltrating immune populations. PDHA1 expression showed significant positive correlations with neutrophils, eosinophils, naive B cells, CD8 T cells, and naive CD4 T cells. CDKN2A expression showed a significant positive correlation with resting CD4 memory T cells, while showing significant negative correlations with neutrophils, eosinophils, naive B cells, CD8 T cells, and naive CD4 T cells. The inverse relationship between CDKN2A expression and neutrophil infiltration is particularly intriguing, suggesting that CDKN2A loss in OS may create a tumor microenvironment that is more permissive to neutrophil recruitment, or alternatively that neutrophil-rich microenvironments suppress CDKN2A expression through paracrine mechanisms.
The immune cell findings from CIBERSORT were subjected to large-scale clinical validation using 25,665 routine blood records collected from the First Affiliated Hospital of Guangxi Medical University between January 2012 and January 2022. Of these, 1,727 were confirmed osteosarcoma patients and 23,938 were healthy controls. Four specific parameters from the complete blood count were analyzed: absolute neutrophil count, neutrophil percentage, absolute lymphocyte count, and lymphocyte percentage. All four parameters showed statistically significant differences between OS and control groups by two independent samples t-test.
Blood test findings: Neutrophil percentage was significantly higher in OS patients compared to healthy controls. Conversely, both absolute lymphocyte count and lymphocyte percentage were significantly lower in the OS group than in healthy controls. These patterns are consistent with an elevated neutrophil-to-lymphocyte ratio (NLR) in osteosarcoma, a phenomenon previously reported in other solid tumors as a marker of systemic inflammation and adverse prognosis. The clinical blood test results corroborate the CIBERSORT-predicted immune cell differences, providing real-world biological plausibility to the computational findings.
Immunohistochemistry: Tissue sections from OS specimens and matched paracancerous (adjacent normal) tissues excised during surgery were stained with validated antibodies against PDHA1 (Proteintech, catalog 18068-1-AP) and CDKN2A (Bioss, catalog bs-0740R). The staining results confirmed that PDHA1 protein expression was significantly higher in paracancerous tissue than in osteosarcoma tissue, meaning PDHA1 is downregulated at the protein level in tumor tissue. CDKN2A protein expression was significantly higher in osteosarcoma compared to paracancerous tissue, meaning CDKN2A is upregulated in OS. Both findings were directionally consistent with the bioinformatic differential expression analysis.
RT-qPCR in OS cell lines: ATCC-sourced OS cell lines (HOS and MG63) and the normal osteoblast cell line hFOB1.19 were used for transcript-level validation. Real-time quantitative PCR using the 2 to the negative delta-delta-Ct method confirmed that PDHA1 mRNA expression in control cells (hFOB1.19) was significantly lower than in OS cell lines HOS and MG63. This appears to be the opposite direction from the IHC result, which showed PDHA1 protein higher in paracancerous tissue than in OS tissue. The authors describe this as consistent with their analysis, though the apparent discrepancy between IHC and PCR directions for PDHA1 warrants careful interpretation regarding post-transcriptional regulation.
Acknowledged limitations: The authors identify two primary limitations. The first is sample size. Although 88 OS cases from UCSC Xena and 396 GTEx normal controls constitute a meaningful training set, the authors note this is small relative to large-sample clinical studies. The validation set from GEO (GSE42352) adds an independent cohort, but with only 84 OS cases and 3 osteoblast controls, external generalizability remains constrained. The second limitation is insufficient laboratory validation depth. The study validated PDHA1 and CDKN2A via immunohistochemistry and RT-qPCR, but functional experiments to determine how these genes mechanistically contribute to OS tumorigenesis and immune evasion were not performed. The absence of functional knockdown or overexpression studies means that while PDHA1 and CDKN2A are associated with OS, their causal roles remain to be established.
The PDHA1 biology in broader context: Pyruvate Dehydrogenase E1 Subunit Alpha 1 encodes the catalytic subunit of the pyruvate dehydrogenase complex (PDC), which sits at the metabolic nexus between glycolysis and the TCA cycle. PDHA1 dysregulation has been documented in prostate cancer (where it drives lipid synthesis via compartmentalized mitochondrial-nuclear shuttling), in FH-deficient cancers (where phosphorylation inhibits PDHA1 to limit TCA carbon entry), and in cancers driven by EGF receptor signaling (where lysine acetylation at K202 of PDHA1 suppresses a regulatory phosphatase). In all cases, impaired PDHA1 activity promotes aerobic glycolysis and the metabolic reprogramming that supports tumor growth. Its role in cuproptosis is logical given that cuproptosis specifically targets lipid-acylated TCA cycle proteins, and PDHA1 is a core TCA component.
CDKN2A as a multifunctional tumor suppressor: CDKN2A encodes p16INK4a (and p14ARF from an alternative reading frame), a canonical CDK4/6 inhibitor that enforces the G1/S cell cycle checkpoint. Its deletion or silencing is among the most common events in human cancer across multiple tumor types, including melanoma, pancreatic cancer, glioma, and now confirmed here in osteosarcoma contexts. The enrichment of CDKN2A in T-cell receptor and NK-cell signaling pathways identified by GSEA in this study aligns with emerging evidence that CDKN2A loss not only drives proliferation but also reshapes the immunosuppressive tumor microenvironment by altering senescence signaling and cytokine secretion patterns.
Future work should include larger multicenter OS datasets to confirm diagnostic AUC performance, functional cell line and animal model studies to determine whether PDHA1 and CDKN2A are causal drivers or downstream consequences of OS transformation, and prospective clinical validation of blood neutrophil-to-lymphocyte ratio as a simple, accessible prognostic tool. Integration of these cuproptosis biomarkers with copper-chelation-based therapeutic approaches, currently being explored in other malignancies, represents a particularly promising translational direction for osteosarcoma.