Osteosarcoma (OS) is the most common primary malignant bone tumor, occurring most frequently around the knee and the proximal humerus. Its incidence follows a bimodal age distribution with peaks in the 10-19 and 60-79 age groups. Although rare at 2-4.2 cases per million population per year, OS carries a sobering prognosis: the 5-year survival rate is approximately 60% for patients without metastases and falls to roughly 20% in patients with metastasis or disease recurrence. Despite advances in surgical technique and multiagent chemotherapy regimens including methotrexate, doxorubicin, and cisplatin, the overall prognosis for OS has not substantially improved over the past several decades.
The role of hypoxia: Tumor hypoxia, a state of reduced oxygen availability within the tumor microenvironment, is a well-established hallmark of most solid tumors. In hypoxic conditions, hypoxia-inducible factors (HIF), particularly HIF-1alpha, become activated and bind to hypoxia-response elements (HRE) in target genes. This transcriptional activation drives a cascade of downstream effects spanning angiogenesis, metabolic reprogramming, cell survival, immune evasion, and metastasis. Critically, hypoxia also promotes chemotherapy resistance by reducing reactive oxygen species production and upregulating drug efflux pumps, which is particularly relevant in OS where chemoresistance is a major driver of treatment failure.
The knowledge gap: Despite extensive research into hypoxia across cancer types, the specific interaction between hypoxia signaling and OS biology had not been comprehensively characterized at the molecular level prior to this study. The authors recognized that identifying hypoxia-regulated gene signatures in OS could simultaneously address two clinical challenges: early accurate diagnosis (distinguishing OS from normal bone tissue and from benign bone lesions) and prognostic stratification (identifying patients at high risk of recurrence or death who might benefit from intensified therapy or investigational agents targeting hypoxia pathways).
This 2024 study published in Scientific Reports presents an integrative bioinformatics approach combining gene expression analysis from multiple GEO datasets, weighted gene co-expression network analysis (WGCNA), Gene Set Enrichment Analysis (GSEA), and LASSO machine learning to derive and experimentally validate two hypoxia-related gene biomarkers, STC2 and TMEM45A, for OS diagnosis and prognosis.
The study assembled data from four Gene Expression Omnibus (GEO) datasets plus the GDC TARGET-OS cohort, covering different aspects of the analysis. The primary discovery dataset was GSE99671, containing gene expression profiles from 18 normal bone tissue samples and 18 OS samples. GSE16088 provided an external validation set with 4 normal individuals and 16 OS patients, while GSE225588 (6 normal, 6 OS) offered additional validation. GSE21257, comprising 53 OS patients without normal controls, was used specifically for prognostic modeling. The TARGET-OS cohort from the UCSC Xena database provided the primary cohort for survival analysis and prognostic model construction, with clinical variables including age, sex, and metastasis status.
Data preprocessing and quality control: Expression data in GSE99671 were subjected to principal component analysis (PCA) and hierarchical clustering to assess intergroup differences and identify outlier samples. This quality control step resulted in the removal of 3 samples from the 36-sample GSE99671 dataset, leaving 33 samples for downstream WGCNA analysis. The DESeq2 package was used for differential expression analysis with cutoff criteria of P less than 0.05 and fold change greater than 1, yielding 218 upregulated and 541 downregulated differentially expressed genes (DEGs) visualized in volcano plots and heatmaps.
Overall analytical workflow: The workflow proceeded in three sequential phases. First, candidate hypoxia-related genes were identified by intersecting hub module genes from WGCNA with 200 hypoxia hallmark genes from the Molecular Signatures Database (MSigDB). Second, LASSO regression was applied to the resulting 21 candidate genes to select a minimal, non-redundant diagnostic gene set. Third, Cox regression modeling was used to build a prognostic score, which was then validated in an independent GEO dataset and tested against clinical variables including metastasis status. Experimental validation of gene expression was performed using quantitative real-time PCR (qRT-PCR) in OS cell lines (HOS, 143B, Saos-2, U-2OS) and normal osteoblast cells (hFOB1.19).
This train-test-external validation architecture is notable for its rigor relative to many single-dataset bioinformatics studies. By using separate discovery, internal validation, and external validation cohorts with distinct patient populations, the authors substantially reduced the risk of overfitting that plagues many published genomic biomarker studies in rare cancers.
Weighted Gene Co-Expression Network Analysis (WGCNA) is a systems biology approach that identifies groups of genes (modules) that are coordinately expressed across samples. Rather than treating each gene as independent, WGCNA builds a weighted adjacency matrix based on pairwise Pearson correlations, then applies a soft-thresholding power to approximate a scale-free network topology. The approach groups genes into modules with similar expression patterns and then correlates each module to clinical traits of interest, in this case OS disease status (tumor vs. normal).
Module construction and selection: The GSE99671 dataset's top 50% of most variably expressed genes were input into WGCNA. The soft-thresholding power was set to 7, chosen as the lowest power achieving scale-free network properties based on scale independence and mean connectivity plots. Dynamic tree cutting identified 26 co-expression modules from this analysis. Module-trait relationship analysis at P less than 0.05 identified modules significantly correlated with OS disease status. Two modules, the brown and pink modules, were selected as hub modules based on their module-trait correlation coefficient and were used for downstream analysis. Genes in these two modules with module membership (MM) greater than 0.5 and gene significance (GS) greater than 0.2 were retained as hub genes, resulting in a combined pool of 752 genes across both modules.
GSEA and hypoxia pathway activation: GSEA on the GSE99671 dataset confirmed that the hypoxia hallmark pathway (from MSigDB) was significantly enriched in OS samples compared to normal samples (P less than 0.001), providing biological justification for focusing on hypoxia-related mechanisms. This activation signal set the stage for the key filtering step: the 752 hub module genes were intersected with the 200 MSigDB hypoxia hallmark genes using a Venn diagram approach. This intersection yielded 21 candidate genes that are both significantly co-expressed in OS disease modules and members of the canonical hypoxia hallmark gene set.
The convergence of co-expression network analysis and pathway enrichment analysis is methodologically sound because it requires candidate genes to satisfy two independent criteria: coordinated expression with other disease-relevant genes (from WGCNA) and membership in a biologically meaningful pathway (from GSEA). Genes satisfying both criteria are more likely to represent genuine mechanistic contributors to OS pathobiology rather than statistical artifacts of any single analysis.
From the 21 candidate hypoxia-related genes identified by the WGCNA-GSEA intersection, LASSO (Least Absolute Shrinkage and Selection Operator) regression was applied to select the minimal set of genes with maximal diagnostic utility. LASSO is a penalized regression technique that shrinks coefficients of less informative variables to zero, performing simultaneous variable selection and regularization. This prevents overfitting in high-dimensional settings where the number of candidate variables exceeds the number of samples, which is a common problem in genomic biomarker studies. The GSE99671 dataset was randomly divided 70/30 into a training set and test set for this analysis.
Model construction and nomogram: LASSO regression in the training set selected exactly two genes, Stanniocalcin 2 (STC2) and Transmembrane Protein 45A (TMEM45A), as the characteristic diagnostic genes for OS. A logistic regression model was then built using these two variables and quantified in a nomogram for clinical application. The nomogram translates the continuous expression values of STC2 and TMEM45A into a combined probability score that clinicians could theoretically use to assess the likelihood that a given tissue sample represents osteosarcoma.
Diagnostic performance across datasets: ROC curve analysis demonstrated excellent discriminatory performance of the two-gene model. In the training set, the area under the ROC curve (AUC) was 0.995, indicating near-perfect separation of OS from normal samples in the data the model was trained on. In the independent test set (30% held-out portion of GSE99671), the AUC was 0.958, representing robust performance on unseen data from the same cohort. Most notably, in the completely external validation cohort GSE16088 (a separate GEO dataset with 4 normal and 16 OS samples), the model achieved an AUC of 1.0, indicating perfect discrimination on an independent patient population. Decision curve analysis (DCA) confirmed net clinical benefit of the model across a range of threshold probabilities in both the training and test sets, while calibration curves demonstrated good agreement between predicted probabilities and observed outcomes.
Biological plausibility: Both selected genes have established roles in hypoxia biology and cancer. TMEM45A, a transmembrane protein of 275 amino acids, is upregulated under hypoxia conditions and has been shown to mediate chemotherapy resistance in breast and liver cancer cells. Prior work demonstrated that TMEM45A upregulation by hypoxia contributes to cisplatin resistance. STC2 (Stanniocalcin 2), a 302-amino acid glycosylated peptide hormone, is directly regulated by HIF-1 transcription factors and plays roles in calcium regulation, glucose metabolism, and angiogenesis. Upregulated STC2 has been reported as a poor prognostic marker in breast, colorectal, gastric, prostate, hepatocellular, and lung cancers.
While the diagnostic model addresses the question of whether a tissue sample is malignant, the prognostic model addresses a different clinical question: among confirmed OS patients, which individuals are at highest risk of death? For this analysis, the authors used the TARGET-OS cohort, which provides OS patient gene expression data linked to overall survival and clinical variables including age, sex, and metastasis status. The cohort was again split 70/30 into training and test sets.
Cox regression pipeline: The 21 hypoxia-related candidate genes were first evaluated by LASSO Cox regression to identify survival-associated genes in the training set. Univariate Cox regression analysis identified both STC2 expression and metastasis status as significantly correlated with overall survival (P less than 0.05). Multivariate Cox regression confirmed that STC2 and metastasis status remained independently prognostic after mutual adjustment. A risk score was calculated for each patient using the weighted sum of STC2 expression multiplied by its multivariate Cox regression coefficient. Patients were then stratified into high-risk and low-risk groups based on the median risk score cutoff, and Kaplan-Meier survival analysis confirmed significantly inferior overall survival in the high-risk group (P less than 0.05 in the TARGET cohort).
Time-dependent AUC performance: ROC analysis using time-dependent AUC methodology assessed the model's ability to predict survival at specific time points. AUC values were 0.96 for 1-year survival, 0.82 for 3-year survival, and 0.80 for 5-year survival in the primary TARGET cohort. These values compare favorably with published prognostic models in other cancers and suggest that STC2 expression carries particularly strong short-term prognostic information, which may reflect its role in driving rapid disease progression through angiogenesis promotion and hypoxia-mediated resistance to early treatment.
External validation in GSE21257: The prognostic model was validated in the independent GSE21257 dataset comprising 53 OS patients. This external cohort confirmed the model's ability to stratify patients into significantly different survival groups (P less than 0.05 by Kaplan-Meier analysis). Time-dependent AUC values in GSE21257 were consistent with the discovery cohort. Calibration curves for 1-, 3-, and 5-year survival probability demonstrated adequate agreement between model predictions and observed outcomes, and decision curve analysis confirmed that using the nomogram for clinical decision-making would provide net benefit across a meaningful range of threshold probabilities.
Bioinformatics models are prone to false discovery if not supported by independent experimental evidence. To validate the expression findings computationally observed in GEO datasets, the authors cultured four OS cell lines (HOS, 143B, Saos-2, and U-2OS) alongside the normal human osteoblast cell line hFOB1.19 and measured STC2 and TMEM45A expression levels by quantitative RT-PCR (qRT-PCR). Expression levels were normalized to ACTB (beta-actin) as an internal reference, and relative expression was calculated using the 2-delta-delta-Ct method. Results showed that OS cell lines expressed significantly higher levels of both STC2 and TMEM45A compared to the normal osteoblast cell line, directionally consistent with the bioinformatics findings across both GSE99671 and GSE225588 datasets. This in vitro confirmation adds mechanistic credibility to the computational observations and supports the hypothesis that these genes are not merely statistically associated with OS but are actively upregulated in OS cells.
Immune cell infiltration analysis: The tumor microenvironment (TME) is increasingly recognized as a major determinant of OS prognosis and treatment response, particularly as immunotherapy gains traction in sarcoma. Patients in the TARGET cohort were divided into high-risk and low-risk groups based on their prognostic risk scores, and immune cell infiltration was estimated using three independent computational methods: EPIC, MCPCOUNTER, and XCELL. These algorithms use gene expression signatures to deconvolve the cellular composition of bulk tumor RNA data.
Cancer-associated fibroblasts and STC2: A consistent finding across all three deconvolution methods was that cancer-associated fibroblast (CAF) abundance was significantly downregulated in the high-risk group relative to the low-risk group. Furthermore, Pearson correlation analysis directly linking STC2 expression to CAF infiltration scores revealed a significant negative association (P less than 0.05), meaning that higher STC2 expression in OS tumors correlates with lower CAF infiltration. The biological interpretation is complex, as CAFs can have dual roles in cancer, sometimes promoting tumor growth through matrix remodeling and sometimes restricting tumor spread. The observation that high STC2 expression associates with low CAF infiltration and poor prognosis raises the hypothesis that STC2 may modulate stromal cell recruitment or differentiation in the OS microenvironment, a relationship warranting dedicated mechanistic study.
To contextualize the role of STC2 beyond osteosarcoma, the authors performed a pan-cancer expression and prognosis analysis using two established online platforms: GEPIA (Gene Expression Profiling Interactive Analysis, hosted by Peking University Cancer Institute) for expression comparisons across TCGA cancer types, and Kaplan-Meier Plotter for survival association analysis. These tools leverage TCGA and GEO data to rapidly assess whether a gene of interest shows consistent dysregulation and prognostic relevance across the full spectrum of human cancers.
Expression patterns across cancer types: Pan-cancer expression analysis confirmed that STC2 is significantly upregulated in multiple cancer types compared to matched normal tissues. The analysis highlighted particularly elevated STC2 expression in esophageal squamous cell carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), clear cell renal cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), and stomach adenocarcinoma (STAD). This pattern of overexpression across multiple epithelial malignancies suggests that STC2 upregulation is not an OS-specific phenomenon but rather a broadly shared feature of hypoxia-driven tumor biology.
Prognostic relevance across cancer types: Kaplan-Meier survival analysis stratified by STC2 expression level (high vs. low) demonstrated that high STC2 expression was significantly associated with poor overall survival in ESCA, HNSC, KIRC, LUSC, and STAD. This cross-cancer prognostic consistency substantially strengthens the argument that STC2 is a functionally important oncogenic driver rather than a passenger event that is incidentally co-expressed with disease-promoting genes. The fact that STC2's prognostic relevance holds across anatomically and histologically diverse cancer types implies that its mechanism of action, likely through HIF-1-mediated transcriptional programs driving angiogenesis and metabolic adaptation, operates independently of tissue-of-origin context.
For OS specifically, the pan-cancer data provides a broader biological framework for interpreting the mechanistic experiments: the chemoresistance and pro-tumorigenic functions of STC2 established in breast, colorectal, and gastric cancer likely apply in osteosarcoma as well, making STC2 an attractive candidate for therapeutic targeting beyond biomarker utility alone.
Database-only discovery cohorts: All gene expression datasets used for biomarker discovery and initial validation were retrospectively downloaded from public databases (GEO and TCGA/TARGET). While this approach is now standard in bioinformatics research, publicly available datasets have known limitations including batch effects between studies, variation in tissue processing and RNA extraction protocols, differences in patient demographics and treatment histories, and potential ascertainment biases in which patients were enrolled in the original studies. The authors acknowledge that prospective clinical studies specifically designed to test the diagnostic and prognostic models will be necessary to confirm real-world applicability.
Functional validation gaps: The study validated that STC2 and TMEM45A are differentially expressed between OS and normal cells in vitro, but functional experiments demonstrating the causal roles of these genes in OS biology were not performed. Specifically, knockdown or knockout of STC2 or TMEM45A in OS cell lines with assessment of proliferation, invasion, migration, apoptosis, and drug sensitivity would establish whether these genes are drivers of OS pathobiology or merely correlated markers. The authors explicitly note that further experiments investigating the functional mechanisms and therapeutic targetability of STC2 and TMEM45A are needed.
Sample size constraints: GSE225588, used as one of the validation datasets, contained only 6 normal and 6 OS samples, which is insufficient to provide statistical power for independent validation. GSE99671 (the discovery dataset after quality control) comprised only 33 samples. While the TARGET-OS prognostic cohort is larger, the relatively small sizes of several component datasets introduce uncertainty about the generalizability of findings. Rare cancer bioinformatics research is inherently constrained by sample availability, but this limitation should be acknowledged when interpreting the reported AUC values, particularly the perfect AUC of 1.0 in the small GSE16088 external validation set (4 normal, 16 OS).
Clinical translation considerations: For the diagnostic model to enter clinical practice, it would need to demonstrate value over and above existing diagnostic tools, including imaging, biopsy histology, and immunohistochemistry. The prognostic model would need prospective multicenter validation demonstrating that risk stratification using STC2 expression improves treatment decisions compared to established prognostic factors such as metastasis status and histological response to neoadjuvant chemotherapy. Integration of these genomic signatures with clinical variables and imaging biomarkers within multimodal prognostic frameworks represents the most promising path toward actionable clinical tools.