Osteosarcoma (OS) is the most common primary malignant bone tumor in adolescents and young adults, with peak incidence between ages 15 and 19. The standard treatment follows a sandwich protocol: neoadjuvant chemotherapy, surgical resection, then adjuvant chemotherapy, using agents such as adriamycin, cisplatin, ifosfamide, and high-dose methotrexate. This regimen has pushed 5-year overall survival above 60% for localized disease. However, outcomes for patients with metastatic or relapsed disease remain dismal, hovering at approximately 20%, a figure that has not improved meaningfully in several decades despite incremental changes in chemotherapy dosing and surgical technique.
The molecular complexity problem: Osteosarcoma has unusually high intra-tumor genomic heterogeneity. Unlike many solid tumors where a small set of driver mutations recur frequently and can be targeted, OS shows a paucity of high-frequency recurrent mutations at the intra-tumor level. This genomic instability makes it difficult to identify actionable targets, and attempts to apply immunotherapy with immune checkpoint inhibitors have yielded disappointing results in unselected OS populations. One reason may be that the tumor microenvironment (TME) in OS is broadly immunosuppressive, and identifying which patients harbor a more permissive immune landscape, one that might respond to immunotherapy or chemotherapy differently, requires better molecular characterization.
The role of the tumour microenvironment: The TME is composed of tumor cells, stromal cells, immune cells (including tumor-infiltrating lymphocytes and macrophages), fibroblasts, vasculature, and extracellular matrix (ECM). In OS, TME components do not passively surround tumor cells; they actively regulate proliferation, migration, metastasis, drug resistance, and immune escape. The balance between pro-inflammatory M1-type macrophages (anti-tumor) and immunosuppressive M2-type macrophages (pro-tumor) is particularly relevant, as is the PD-1/PD-L1 axis that allows tumor cells to evade T-cell killing. A meta-analysis cited in this paper found that PD-L1 overexpression predicted poor overall survival (HR 1.45, 95% CI 1.11-1.90), worse metastasis-free survival (HR 1.58), and worse event-free survival (HR 2.82) in OS, as well as a higher rate of tumor metastasis (OR 2.95).
Given this background, the study set out to build a data-driven prognostic risk model anchored specifically in TME gene expression profiles, with a particular focus on the adolescent and young adult OS population where age-specific immune landscapes differ from adult OS. The key question: can a small set of immune and stromal genes, selected from a large expression dataset using rigorous statistical methods, reliably stratify OS patients into high-risk and low-risk groups?
The study used two publicly available datasets: TARGET (Therapeutically Applicable Research to Generate Effective Treatments) as the primary training cohort and GSE21257 from the Gene Expression Omnibus (GEO) database as the external validation cohort. From TARGET, 82 patients aged 25 years or younger with complete gene expression data and at least one month of follow-up were enrolled. The GSE21257 dataset contributed 53 OS patients with clinical information including age, gender, histological subtype, Huvos grade (the standard histopathologic grading of tumor necrosis after neoadjuvant chemotherapy), tumor location, metastasis at diagnosis, and survival data. All samples came from primary lesions, and gene expression profiling was performed by microarray. The validation dataset was normalized using the R package "limma."
ESTIMATE algorithm for TME scoring: TME immune and stromal scores were calculated from mRNA expression profiles using the ESTIMATE algorithm, which was originally designed for and validated across multiple solid tumor types including ovarian cancer, bladder cancer, and gastric cancer. Patients were split into high and low immune score groups and high and low stromal score groups based on the median of each score. This stratification revealed that high immune scores were associated with a 5-year overall survival of 82.6% versus 48.7% in the low immune score group (P = 0.03). High stromal scores corresponded to a 5-year OS of 65% versus 47.3% for low stromal scores (P = 0.053). Recurrence-free survival (RFS) showed parallel trends: RFS rates of 65% versus 47.3% for high versus low immune score groups (P = 0.06).
DEG identification and functional enrichment: Differential gene expression analysis between high and low immune score groups and high and low stromal score groups was performed using the limma R package, with filtering criteria of fold-change greater than 1 (or less than -1) and adjusted P less than 0.05. The high immune score group yielded 32 downregulated and 262 upregulated genes; the high stromal score group produced 252 downregulated and 16 upregulated genes. Venn diagram overlap identified 31 commonly upregulated and 12 commonly downregulated genes. Combining these overlap genes with the top differentially expressed genes from each individual comparison yielded 122 TME-related differentially expressed genes (tmDEGs) for further analysis. Gene Ontology (GO) enrichment confirmed these tmDEGs were involved in leukocyte activation, cell-cell adhesion, and antigen binding. KEGG pathway analysis implicated osteoclast differentiation, B-cell receptor signaling, cytokine-cytokine receptor interaction, and JAK-STAT signaling as key enriched pathways.
Three-stage gene selection pipeline: Prognostic gene selection proceeded through three sequential steps. First, univariate Cox regression on the 122 tmDEGs identified 32 genes significantly correlated with overall survival (P less than 0.05). Second, LASSO (least absolute shrinkage and selection operator) regression on those 32 genes identified 13 candidate genes based on their relative regression coefficients. Third, XGBoost machine learning ranked feature importance across all 122 tmDEGs and identified the top 3 most predictive genes. The intersection of the LASSO-selected 13 genes and the XGBoost top features yielded 9 candidate genes. Finally, multivariate Cox regression analysis with Akaike's Information Criterion (AIC) minimization selected the final 3-gene combination with AIC value of 188.9.
The final three-gene risk model (3-GRM) combines expression levels of COCH (Cochlin), MYOM2 (Myomesin 2), and PDE1B (Phosphodiesterase 1B) using coefficients derived from multivariate Cox regression. The scoring formula is: risk score = [0.30 x COCH expression level] + [0.41 x MYOM2 expression level] - [2.10 x PDE1B expression level]. The concordance index (C-index) of the model was 0.77 (SE = 0.042, 95% CI: 0.688-0.852), indicating good discriminative ability for overall survival prediction. Patients in the TARGET training cohort were dichotomized into high and low 3-GRM score groups based on the median score cutoff.
COCH (Cochlin): COCH is highly expressed in sensory organs including the inner ear and eye, and is also found in lymph nodes and spleen. The protein encoded by COCH is the main non-collagenous component of the extracellular matrix and has high affinity for multiple collagen types. In the context of cancer, COCH overexpression has been linked to stage progression in clear cell renal cell carcinoma. The N-terminal LCCL domain of COCH can signal innate immune cells and amplify cytokine responses. Cochlin secreted from follicular dendritic cells in the spleen drives systemic immune responses by inducing IL-1-beta and IL-6 secretion and enhancing recruitment of neutrophils and macrophages. In the 3-GRM, COCH overexpression in the high-score group is associated with poor prognosis, and the authors hypothesize this may reflect COCH-mediated immune microenvironment imbalance and immune escape promotion.
MYOM2 (Myomesin 2): MYOM2 encodes a muscle fiber structural protein that is an essential component of the cytoskeleton and the M-band at the center of the sarcomere, critical for sarcomere contraction stability. Its role in cancer is less established, but it has been identified as a recurrently mutated gene in relapsed or refractory diffuse large B-cell lymphoma, suggesting a connection to oncogenic events. MYOM2 overexpression has also been associated with inflammation: it is the only significantly upregulated gene in localized invasive periodontitis. Research on ankylosing spondylitis found that IgG galactosylation status combined with the MYOM2 rs229466 polymorphism (T allele leads to MYOM2 overexpression) predicts anti-TNF therapy response, indicating that MYOM2 interacts with inflammatory signaling. Like COCH, MYOM2 is overexpressed in the high 3-GRM score group and associated with poor prognosis; the proposed mechanism involves modulation of the TME inflammatory response to create a more pro-tumor microenvironment.
PDE1B (Phosphodiesterase 1B): PDE1B is located on chromosome 12q13 and encodes a member of the cyclic nucleotide phosphodiesterase family, enzymes that degrade intracellular cAMP and cGMP. These second messengers play important roles in immune cell differentiation. Studies show PDE1B regulates differentiation of multiple immune cell types. Granulocyte macrophage colony stimulating factor (GM-CSF) can shift monocyte differentiation toward dendritic cells and also upregulates PDE1B; conversely, IL-4 treatment in the presence of GM-CSF suppresses PDE1B2 upregulation, altering macrophage phenotype toward increased phagocytosis and leukocyte adhesion molecule CD11b expression. In the 3-GRM, PDE1B is downregulated in the high-score group, making it a protective factor. The negative coefficient (-2.10) in the risk formula reflects this: higher PDE1B expression lowers the risk score. The authors speculate that PDE1B upregulation may promote M1-type macrophage differentiation under GM-CSF influence, enhancing anti-tumor immune responses.
In the TARGET training cohort of 82 patients, dividing by the median 3-GRM score produced two groups with sharply different survival outcomes. Kaplan-Meier analysis showed that the high 3-GRM score group had significantly worse overall survival than the low-score group (P less than 0.05). The 5-year survival rates were 51.3% in the high-score group versus 80.1% in the low-score group, a difference that was both statistically significant and clinically meaningful for a disease where such separation translates directly into treatment decision-making.
ROC curve performance: The discriminative ability of the 3-GRM was evaluated using time-dependent ROC curves. The area under the ROC curve (AUC) values for predicting 1-year, 3-year, and 5-year survival were 0.890, 0.822, and 0.773, respectively. These values indicate that the model performs well at near-term prediction (1-year AUC nearly 0.90) and retains clinically useful predictive power at 5 years (AUC 0.77), suggesting the risk stratification is not a transient phenomenon but reflects a durable underlying biological difference between groups.
Independence from clinical covariates: Univariate Cox regression of all available clinical and molecular variables showed that high 3-GRM score, metastasis at diagnosis, and lung metastasis were all statistically significant risk factors for prognosis. Age, gender, race, primary site location, and surgical method were not independently associated with outcome. Because metastasis at diagnosis and lung metastasis are closely interrelated in the TARGET database (93% of patients with metastasis at diagnosis had lung metastasis), these were combined as a single "metastasis at diagnosis" variable in subsequent multivariate analysis. Multivariate Cox regression confirmed that both high 3-GRM score and metastasis at diagnosis were independent prognostic factors. Patients with metastases at diagnosis had a 5-year survival of 31.3% versus 77.4% for those without (P less than 0.01).
The 3-GRM score was also significantly correlated with the underlying TME metrics derived from the ESTIMATE algorithm. Correlation analysis showed that higher 3-GRM score was associated with lower immune score (R = -0.398, P less than 0.01) and lower stromal score (R = -0.523, P less than 0.01), confirming that the risk model captures the immunological underpinnings of the TME rather than an unrelated biological signal. This internal validation of biological coherence strengthens confidence that the model is not simply an overfitted statistical artifact.
Replication in an independent external dataset is the most demanding test for any prognostic model. The GSE21257 validation cohort comprised 53 OS patients from GEO, with clinical annotation including Huvos grade, metastasis at diagnosis, and survival. Applying the same 3-GRM scoring formula to this cohort and splitting patients at the same median threshold resulted in 28 patients in the high-score group and 25 in the low-score group. Kaplan-Meier curves confirmed that patients with high 3-GRM scores had significantly shorter overall survival than those with low scores, directionally consistent with the training cohort results.
Validation AUC values: Time-dependent ROC curves in the validation cohort yielded AUC values of 0.861 at 1 year, 0.710 at 3 years, and 0.694 at 5 years. The 1-year AUC of 0.861 indicates strong near-term discriminative ability, very close to the training cohort performance of 0.890. The 3-year and 5-year AUC values show a modest reduction compared to training (from 0.822 to 0.710 at 3 years, and 0.773 to 0.694 at 5 years), which is a common and expected consequence of model generalization to a smaller, independent sample. Crucially, the AUC values remain above 0.69 even at 5 years, indicating that the model retains meaningful predictive power beyond the training context.
Independent prognostic factor confirmation: Univariate Cox regression in the GSE21257 cohort identified three significant prognostic variables: Huvos grade (histological response grade after chemotherapy), metastasis at diagnosis, and high 3-GRM score. Multivariate Cox regression confirmed all three remained independent prognostic factors in this external dataset. The appearance of Huvos grade as a significant variable in the validation cohort but not in the training cohort likely reflects the more complete clinical annotation available in GSE21257, including post-chemotherapy histopathologic response data that was not uniformly available for the TARGET cohort.
The consistency of the 3-GRM's prognostic value across both datasets, despite differences in cohort size, patient demographics, clinical annotation completeness, and microarray platform, constitutes meaningful evidence of model robustness. The model performed with similar directionality and comparable magnitude in two datasets collected independently, reducing the likelihood that the prognostic signal is driven by training dataset-specific confounders.
To understand why the 3-GRM score stratifies prognosis so effectively, the authors characterized the immune cell composition of the TME in high versus low 3-GRM score groups using two complementary computational deconvolution methods. CIBERSORT, a deconvolution tool based on gene expression signatures, estimated the proportions of 22 distinct tumor-infiltrating immune cell (TIIC) types across the 82 TARGET samples. xCELL, which uses a gene set enrichment approach, calculated enrichment scores for 64 immune and stromal cell types. Using two orthogonal methods on the same samples provides a form of technical replication for the immune infiltration findings.
Dominant immune cell populations: CIBERSORT analysis of the full 82-patient cohort revealed that the most abundant TIICs were M2 macrophages (27.8% plus or minus 1.14%), M0 macrophages (27.63% plus or minus 11.35%), and CD4-positive memory T cells (17.8% plus or minus 7.86%), collectively accounting for more than 70% of immune infiltration. This macrophage-dominant landscape is consistent with other reports describing OS as having a heavily macrophage-infiltrated TME. The high prevalence of M2-type macrophages (classically immunosuppressive) relative to M1-type (classically anti-tumor) macrophages is consistent with the generally poor immunotherapy response seen in unselected OS populations.
Differences between high and low 3-GRM score groups: Among the 22 CIBERSORT cell types, the proportion of M0 macrophages was significantly higher in the high 3-GRM score group (P less than 0.01). Conversely, monocytes (P = 0.05), M1 macrophages (P = 0.04), and M2 macrophages (P = 0.02) were all significantly lower in the high-score group. The shift toward undifferentiated M0 macrophages and away from both M1 and M2 subtypes in high-risk patients suggests a stalled or dysregulated macrophage differentiation state rather than simple M2 polarization, which may reflect deeper immunosuppression in the highest-risk tumors.
xCELL confirmation: Across the 64 cell types evaluated by xCELL, the high 3-GRM score group showed significantly lower enrichment of multiple immune and stromal cell populations including activated dendritic cells (aDC), conventional dendritic cells (cDC), immature dendritic cells (iDC), M1 macrophages, M2 macrophages, endothelial cells, and mesenchymal stem cells (MSC). Notably, "muscle cells" and "skeletal muscle cells" were significantly enriched in the high 3-GRM score group, consistent with the upregulation of MYOM2 (a muscle structural protein) in this group. The convergent finding of broadly suppressed immune infiltration in high-risk patients across both CIBERSORT and xCELL provides mechanistic grounding for why this group has worse outcomes.
Nomograms are graphical tools that integrate multiple prognostic variables into a single visual prediction framework, producing individualized probability estimates for clinical outcomes. Since multivariate Cox regression confirmed that both 3-GRM score and metastasis at diagnosis are independent prognostic factors for OS in osteosarcoma, the authors constructed a nomogram incorporating these two variables to predict 1-year, 3-year, and 5-year overall survival probabilities. The nomogram was built using the "rms" R package and calibrated using bootstrap resampling with 1,000 iterations to evaluate the consistency between predicted and actual survival rates.
Nomogram performance in the training cohort: The nomogram achieved a C-index of 0.825 in the TARGET training cohort, substantially higher than the 3-GRM alone (C-index 0.77). The AUC values for predicting 1-year, 3-year, and 5-year survival reached 0.971, 0.853, and 0.818, respectively. These values represent a meaningful improvement over the 3-GRM alone at every time point (1-year: 0.971 vs 0.890; 3-year: 0.853 vs 0.822; 5-year: 0.818 vs 0.773). The calibration plots showed good agreement between nomogram-predicted survival probabilities and actual observed survival rates across all three time horizons, confirming that the model is well-calibrated rather than systematically over- or underestimating risk.
Nomogram performance in the validation cohort: The nomogram was also validated in the GSE21257 dataset. AUC values in the validation cohort were 0.781 at 1 year, 0.840 at 3 years, and 0.795 at 5 years. The calibration plots in the validation cohort also showed good alignment between predicted and actual survival, and the predicted results were consistent with actual observations in both supplementary calibration plots (Supplementary Figures 3A and 3B). Interestingly, the 3-year and 5-year AUC values in the validation cohort (0.840 and 0.795) are higher than or comparable to those in the training cohort (0.853 and 0.818), suggesting the nomogram generalizes well rather than deteriorating on external data.
Comparison with prior OS prognostic models: The paper benchmarks the nomogram against several published OS prognostic models from the literature, reporting AUC values at each time point. The nomogram model outperforms all comparators in this analysis at each of the three time horizons. This comparison is necessarily limited by the fact that each model was tested on the same dataset (TARGET), which may favor the locally trained model, but the consistency across time points and the validation in GSE21257 provides some reassurance that the gains are real.
Dataset size and clinical annotation gaps: The most significant limitation is the relatively small training cohort of 82 patients, which restricts the statistical power of the modeling pipeline and leaves some clinical covariates (such as tumor location, histological grade, and formal clinical staging) underrepresented in the analysis. These variables were not included in the multivariate Cox models partly because they were not available with sufficient completeness in the TARGET dataset. The absence of staging and grade from the final model means the nomogram cannot currently be applied to individualized staging-adjusted risk assessment in clinical practice without supplemental data.
Mechanistic gaps for the three genes: While each of COCH, MYOM2, and PDE1B has established functions in immune and stromal biology, their specific mechanisms of action in osteosarcoma remain unclear. The authors explicitly acknowledge this gap: none of the three genes has been previously reported in the context of OS prognosis, and the pathways through which they collectively shape the TME and influence clinical outcomes have not been experimentally validated in OS cell lines or animal models. The mechanistic hypotheses proposed in the discussion section (COCH-mediated cytokine imbalance, MYOM2-mediated inflammatory phenotype modulation, PDE1B-driven M1 macrophage differentiation) are plausible based on known biology in other contexts but remain speculative for OS specifically.
Retrospective design and need for prospective validation: Both the training and validation datasets are retrospective cohorts from publicly available gene expression repositories. All expression data were from microarray platforms, and no RNA-seq validation was performed. The prognostic cutoff (median of the training cohort's risk score distribution) was not prospectively defined or validated in a clinically independent prospective study. The authors explicitly call for multi-center clinical trials and prospective studies to validate the 3-GRM before any clinical application, an appropriately cautious position given the constraints of the available data.
Future directions and clinical potential: If validated prospectively, the 3-GRM represents a potentially low-cost complement to standard imaging and pathology assessment for OS prognosis. The three-gene panel could, in principle, be converted from microarray expression data into a quantitative RT-PCR assay deployable on standard clinical tumor samples. The immune cell infiltration findings also open avenues for patient selection in immunotherapy trials: patients with low 3-GRM scores (and correspondingly higher immune infiltration, more M1 macrophages, and intact dendritic cell populations) may represent the subset of OS patients most likely to benefit from PD-1/PD-L1 blockade or macrophage-targeting strategies. Integration of the 3-GRM with existing clinical risk factors such as Huvos grade and metastasis status could form the basis for a more comprehensive OS risk stratification tool with genuine utility in treatment planning.