Osteosarcoma is the most common primary malignant bone tumor, with peak incidence during adolescent growth spurts and ranking eighth among all pediatric cancers by total incidence. The tumor preferentially arises at the ends of long bones, particularly the distal femur, proximal tibia, and proximal humerus, where rapid bone remodeling during growth creates a permissive environment for oncogenic transformation. Despite its relatively low absolute incidence, osteosarcoma carries a disproportionately high mortality burden because of its aggressive growth pattern and strong propensity for early hematogenous metastasis, predominantly to the lungs.
The metastasis problem: Approximately 10-20% of patients present with detectable metastatic lesions at diagnosis, and an additional 30-40% of patients who initially appear to have localized disease will develop metastases during or after treatment. The five-year survival rate for localized osteosarcoma is approximately 60-70%, but drops to below 30% in the metastatic setting. This stark survival difference makes early identification of patients at high metastatic risk a central clinical and research priority.
Treatment landscape: Standard treatment has remained largely unchanged for four decades, consisting of surgical resection combined with multi-agent chemotherapy using doxorubicin, cisplatin, and high-dose methotrexate. While this regimen controls primary tumor growth effectively, it has not meaningfully improved outcomes for metastatic patients. The stagnation in treatment progress reflects, in part, a fundamental gap in understanding the molecular mechanisms that drive metastatic progression, particularly at the level of tumor evolution and clonal dynamics.
This 2025 study in the Journal of Translational Medicine takes a genomic evolutionary approach to this problem, using whole-exome sequencing from the publicly available TARGET database to reconstruct how osteosarcoma tumors evolve, which mutations accumulate early versus late, and whether the pattern of clonal evolution can predict metastatic behavior. The goal is not just biological insight but a deployable machine learning classifier that stratifies patients by metastatic risk at diagnosis.
The study is a reanalysis of publicly available genomic data rather than a prospective trial. The primary dataset was drawn from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database, accessed via the Genomic Data Commons (GDC) platform. From the full TARGET osteosarcoma cohort, the authors assembled a working set of 61 patients who had matching somatic mutation data, copy number variation (CNV) data, and clinical information. Of these 61, 19 presented with metastasis at diagnosis and 42 did not.
Variant calling pipeline: For 13 patients with available raw whole-exome sequencing data (published by Xu et al.), the authors performed de novo variant calling. Sequencing reads were aligned to the human reference genome hg38 using the Burrows-Wheeler Aligner (BWA, v.0.7.17). Somatic mutations were called using Mutect2 following GATK (v.4.3.0.0) best practices, with duplicate marking via MarkDuplicates and base quality score recalibration applied beforehand. Copy number variants were detected with CNVkit (v.0.9.10), and all mutations were annotated using GATK Funcotator.
Tumor purity estimation: A significant methodological challenge in the TARGET dataset is that existing tumor purity estimation tools such as ABSOLUTE and FACETS cannot be directly applied because the database lacks the required data formats. The authors addressed this by implementing a method from Locallo et al. that estimates tumor purity from the variant allele frequency (VAF) distribution of single-nucleotide variants (SNVs) in copy number-neutral genomic regions. The DBSCAN clustering algorithm was used to identify the dominant VAF cluster, with the k-distance method for setting the neighborhood radius parameter and minPts empirically set to the feature dimension plus one.
Cancer cell fraction and clonal reconstruction: For each mutation, the cancer cell fraction (CCF) was calculated by combining VAF with tumor purity and local somatic copy number. A mutation with CCF close to 1 was classified as a founder clone mutation present in nearly all tumor cells; CCF below 1 indicated a subclonal mutation present in only a subset. Subclonal architectures and evolutionary trees were reconstructed using the R package RETCHER, which employs a beta mixture model with Bayesian variational inference. A maximum of 10 clusters was allowed, with a minimum of three mutations per cluster. The resulting clonal trees were inferred using sum and cross rule constraints on cluster proportions.
Mutational signatures: De novo mutational signatures were identified using the non-negative matrix factorization module in maftools (v.2.14.0), then mapped to known signatures in the COSMIC database by cosine similarity. Signatures with cosine similarity greater than 0.5 were considered interpretable.
The genomic characterization of the 61 TARGET-OS patients revealed a mutational landscape consistent with the known biology of osteosarcoma as a predominantly structural variant-driven malignancy with relatively low point mutation burden. The median number of non-synonymous mutations was 19 (range 7-60) in the metastatic group and 18 (range 2-77) in the non-metastatic group. Tumor mutational burden (TMB) was estimated at approximately 0.45 mutations per megabase, which is notably low compared to carcinomas and consistent with the characterization of osteosarcoma as a "cold" tumor immunologically.
Mutational signatures: Using non-negative matrix factorization combined with cosine similarity mapping to COSMIC signatures, two robust mutational signatures were identified across the cohort: SBS5 and SBS40. Both are classified as "unknown" in COSMIC in terms of etiology, though SBS5 is considered a clock-like signature associated with the accumulation of mutations over time. A key supporting observation came from a large study of 27 pediatric cancer types by Venu et al., which found SBS5 present in 96.6% of pediatric tumor samples and significantly correlated with age at diagnosis. SBS40 was found to be active in nine pediatric cancer types including osteosarcoma, suggesting a shared oncogenic mechanism across pediatric bone and soft tissue tumors.
Most frequent mutations: The top recurrently mutated genes across the cohort were dominated by known osteosarcoma drivers. TP53 mutations, long established as key driver events in osteosarcoma, appeared among the most frequent somatic mutations. RB1 mutations were also prevalent, consistent with the tumor's frequent origin in rapidly dividing osteoblastic cells during skeletal growth. Structural genomic features including copy number amplifications and deletions were heterogeneous across patients, with the CNV heatmap showing substantial inter-patient variation rather than a shared pattern of focal amplifications or deletions.
The broad distribution of mutation counts (ranging from 2 to 77 non-synonymous mutations) highlights the high degree of inter-patient genomic heterogeneity in osteosarcoma, a characteristic that makes identifying universally actionable driver mutations challenging. The relatively low TMB also suggests that immune checkpoint inhibitor approaches, which perform best in high-TMB cancers, face an inherent biological barrier in osteosarcoma.
The clonal evolution analysis using RETCHER revealed that 62% (38/61) of osteosarcoma patients exhibited a simple linear evolutionary pattern, in which successive mutations accumulate in a single dominant clone that expands sequentially. The remaining 38% showed branching evolutionary patterns or contained only a single clone cluster. The linear pattern was more common in non-metastatic tumors, while metastatic tumors showed a higher prevalence of branching evolution, consistent with the idea that greater subclonal diversity provides more opportunities for a subclone to acquire metastatic capability.
CCF distribution and clonal complexity: Across the cohort, most mutations had high CCF values, indicating that they occurred early in tumor development and were present in the majority of tumor cells. Only a few samples had a median CCF below 0.8. Analysis of subclonal cluster counts showed that 55% (12/40) of non-metastatic patients had only one or two subclonal clusters, while this proportion was lower in the metastatic group, consistent with greater clonal complexity in tumors that had spread. This quantitative difference in clonal architecture, rather than simply the presence or absence of specific mutations, carries prognostic information.
ATRX: a key early metastatic driver: The most functionally significant finding from the clonal analysis was the contrasting behavior of ATRX mutations across metastatic and non-metastatic tumors. In a non-metastatic representative sample (PALKDP), the ATRX mutation had a CCF of 0.35, placing it in a subclone that represented only about one-third of tumor cells. In contrast, in a metastatic representative sample (PATMXR), the ATRX mutation had a CCF of 0.93, placing it in the founder clone present in virtually all tumor cells. This finding suggests that when ATRX mutations occur early, as clonal rather than subclonal events, they confer a strong selective advantage that may drive metastatic progression.
ATRX biology: ATRX is a chromatin remodeling factor whose functional loss is closely associated with genomic instability. ATRX mutations activate the alternative lengthening of telomeres (ALT) mechanism, allowing tumor cells to maintain telomere length without telomerase, which is associated with more aggressive behavior in osteosarcoma and other pediatric cancers. Loss of ATRX also promotes NF-kB signaling activation and enhanced integrin-mediated cell adhesion, both of which facilitate invasion and metastasis. The study's evolutionary analysis provides in vivo evidence that ATRX mutations, when they occur early enough to achieve high clonal prevalence, reshape the entire tumor's clonal dynamics in a direction favoring metastatic spread.
The classifier development began with a differential mutation analysis comparing the 19 metastatic and 42 non-metastatic TARGET-OS patients using the mafCompare function in maftools. This analysis considered only non-synonymous mutations present in at least two samples to avoid noise from private mutations. The result was 62 differentially mutated genes. The CCF values of these 62 genes were then used as feature vectors, with metastatic status as the label, in a LASSO (Least Absolute Shrinkage and Selection Operator) regression model to perform feature selection and reduce dimensionality.
Feature selection and the eight key mutations: LASSO regularization identified eight key genes whose CCF values were most predictive of metastatic status: ATRX, RB1, DNAH9, RYR2, ADH1A, ARHGEF15, CARMIL1, and ZMAT1. These eight features were then standardized to ensure balanced scale and stability across classifiers. The biological roles of these genes vary: ATRX and RB1 are established osteosarcoma drivers; DNAH9 encodes a dynein axonemal heavy chain involved in cytoskeletal organization; RYR2 is a ryanodine receptor with emerging roles in calcium signaling in cancer; ADH1A is an alcohol dehydrogenase implicated in metabolic reprogramming; and ARHGEF15, CARMIL1, and ZMAT1 represent less-characterized genes with potential roles in cytoskeletal dynamics and gene regulation.
Model comparison and performance: Four classical machine learning algorithms were compared: logistic regression (LR), decision tree (DT), support vector machine (SVM), and random forest (RF). All models were optimized using random search and grid search for hyperparameter tuning, and evaluated by five-fold cross-validation. LR achieved the highest accuracy (0.83, SD 0.10), recall (0.58, SD 0.49), and F1-score (0.67, SD 0.29) with AUC of 0.72. SVM achieved the highest AUC (0.76, SD 0.25) with accuracy 0.80. DT performed worst with accuracy 0.72 and AUC 0.58. Random forest was intermediate at accuracy 0.75 and AUC 0.65.
Final model evaluation: LR was selected as the final model and evaluated on a held-out test set using a 70:30 train-test split (43 training cases, 18 test cases). The LR model achieved 83% accuracy on the internal test set, confirming that the cross-validation results were not simply artifacts of the small training set. The eight-feature CCF-based logistic regression classifier thus represents the primary deliverable of the study: a biologically grounded, machine learning-validated tool for predicting metastatic status at diagnosis.
External validation used whole-exome sequencing data from 13 osteosarcoma patients published by Xu et al., which was not part of the TARGET-OS training or test sets. This cohort had an important distinguishing characteristic: all 13 patients were confirmed to have developed lung metastasis, making it a uniformly high-risk validation set rather than a mixed metastatic/non-metastatic cohort. The genomic profile of this external cohort differed substantially from TARGET-OS: the median number of non-synonymous mutations was 1,729, ranging from 10 to 5,313, an order of magnitude higher than the TARGET dataset, reflecting either biological differences or differences in sequencing depth and analytical pipelines.
Genomic differences in the validation cohort: Missense mutations were the most frequent mutation type in the external cohort, and none of the common high-frequency mutations from the TARGET-OS analysis appeared in the top 10 most frequent mutations of the Xu et al. dataset. This lack of direct overlap in dominant mutations is notable and highlights the substantial inter-cohort heterogeneity in osteosarcoma genomics. Despite this, the CNV heatmap of the external cohort showed recurrent patterns of chromosomal copy number alterations, including regions of consistent amplification and deletion across patients.
Subclonal reconstruction in the validation cohort: RETCHER was applied to one representative sample (OS12T) from the external cohort to demonstrate that the clonal reconstruction methodology generalized to data from a different source. The ATRX mutation was identified and localized to a specific clonal cluster in this sample, and the evolutionary analysis of the 13-patient external cohort produced coherent clonal trees consistent with the patterns observed in TARGET-OS data, supporting the biological reproducibility of the analytical framework.
Model performance on external data: The logistic regression classifier, trained on TARGET-OS data, achieved 75% accuracy on the external validation cohort. This 8-percentage-point drop from the 83% internal test accuracy is modest and within expectations for a small-sample genomic model applied to a cohort with different data collection methods and genomic characteristics. The external validation result importantly also demonstrates that the model can correctly classify cases outside the original dataset, a prerequisite for any clinical tool aspiring to deployment beyond a single institution.
Beyond the classifier, the study employed Suppes' probabilistic theory of causality to infer the temporal ordering of shared high-frequency mutations at the cohort level. This approach is well-suited to cross-sectional genomic data where longitudinal tumor sampling is unavailable. The method defines a causal relationship between mutation A and mutation B if A occurs before B in evolutionary time (i.e., A has higher CCF than B) and if the occurrence of A increases the probability of B occurring. Mutations present in at least three patients were defined as high-frequency; 12 such mutations were shared between the metastatic and non-metastatic groups.
Non-metastatic evolutionary pathway: In the non-metastatic group (n=23 for this analysis), the evolutionary pathway was centered on TP53 and RB1 as initiating mutations. Five downstream mutated genes, including ATRX, were identified as descendants of TP53 in the inferred evolutionary tree. PCDH15 was identified as a descendant of RB1. This topology, with TP53 and RB1 as early drivers and other mutations accumulating downstream, reflects the classical tumor suppressor gene-driven model of osteosarcoma progression, consistent with the established oncogenic role of TP53 inactivation in initiating osteosarcoma.
Metastatic evolutionary pathway: The metastatic group (n=14 for this analysis) showed a strikingly different topology. Rather than TP53 and RB1 as the primary initiating events, the metastatic evolutionary pathway began with ATRX mutations as an early, initiating event, with other mutations positioned downstream. This inversion of the evolutionary hierarchy, where ATRX moves from a downstream consequence of TP53 mutation (in non-metastatic tumors) to an early clonal driver (in metastatic tumors), is the most striking finding of the cohort-level analysis. PCDH15 and ZFHX3 were observed to have subclonal features specifically in metastatic tumors, suggesting they may contribute to metastatic progression during later stages of tumor evolution.
These evolutionary pathway differences between metastatic and non-metastatic groups provide mechanistic context for the classifier. The eight CCF-based features capture not just which mutations are present, but implicitly encode information about when those mutations occurred in the tumor's evolutionary history, which is precisely the information that distinguishes the two groups.
Sample size constraints: The study's most significant limitation is its small sample size. The primary cohort of 61 patients, with only 19 in the metastatic group, creates a class imbalance that challenges machine learning model training and makes it difficult to reliably estimate performance metrics. The standard deviation values on the cross-validation performance metrics reflect this instability: recall had a SD of 0.49 for the logistic regression model, indicating that performance varied substantially across folds. The external validation dataset of 13 patients, all with confirmed metastasis, provides directional confirmation of generalizability but does not constitute the kind of balanced independent validation required for clinical tool development.
Cross-sectional data limitation: The study relies entirely on cross-sectional data: single biopsies taken at diagnosis, without longitudinal sampling at disease progression, metastasis development, or post-treatment recurrence. This design prevents evaluation of how the tumor's clonal landscape changes over time in response to chemotherapy or how early subclonal mutations evolve into dominant clones during metastatic progression. Longitudinal whole-exome or whole-genome sequencing at multiple timepoints, analogous to studies performed in AML and other leukemias, would substantially strengthen the evolutionary conclusions.
Lack of transcriptomic and multi-omics integration: The current model relies solely on somatic mutation CCF values. It does not incorporate gene expression data, DNA methylation, protein abundance, or metabolomic measurements, all of which carry complementary prognostic information. In particular, transcriptomic data would allow direct measurement of the functional consequences of mutations like ATRX loss, rather than inferring these effects from genomic data alone. Integration of single-cell RNA sequencing would enable mapping of subclonal mutations onto specific cell states within the tumor microenvironment.
Future directions: The authors identify several priority areas for follow-up. Expanding the cohort by integrating data from multiple osteosarcoma genomic repositories, including the ICGC and additional institutional sequencing efforts, is the most immediate need. Combining the mutational CCF classifier with imaging-based biomarkers, particularly MRI radiomics and metabolic PET imaging features, could improve predictive accuracy by capturing complementary spatial information about tumor heterogeneity. The integration of liquid biopsy approaches, particularly circulating tumor DNA (ctDNA) from serial blood draws, could ultimately enable non-invasive, dynamic risk stratification as treatment proceeds, capturing the real-time evolution of tumor clones in response to chemotherapy.