Osteosarcoma (OS) is the most common primary bone malignancy in children and young adults, accounting for roughly 5% of all pediatric malignancies and 20% of bone-related cancers. Despite the availability of surgery, chemotherapy, and radiotherapy, osteosarcoma continues to carry a poor prognosis because of its aggressive biological behavior and high metastatic potential. Five-year survival for patients with localized disease is approximately 60-70%, but this drops to below 30% in patients who present with or develop metastatic disease. The molecular heterogeneity of OS has long complicated efforts to identify actionable therapeutic targets and reliable prognostic biomarkers.
Post-translational modifications (PTMs): PTMs are chemical changes that occur to proteins after they are synthesized, fundamentally altering their function, stability, localization, and interaction with other molecules. The 11 types of PTMs investigated in this study include ubiquitination, methylation, phosphorylation, glycosylation, acetylation, SUMOylation, citrullination, neddylation, palmitoylation, ADP-ribosylation, and succinylation. A growing body of evidence links PTM dysregulation to OS pathogenesis: H3K27 acetylation activates COL6A1 to promote lung metastasis, ALKBH1-regulated N6-methyladenosine modification of TRAF1 mRNA enhances OS proliferation and invasion, and knockdown of SUMO-specific protease 1 suppresses OS cell viability and invasive capacity.
The knowledge gap: Despite the biological importance of PTMs in OS, their roles at single-cell resolution and within the spatial architecture of the tumor microenvironment (TME) had not been comprehensively mapped. The cellular heterogeneity of OS tumors means that bulk RNA sequencing averages signals across diverse cell types, obscuring which specific cell populations are most strongly shaped by PTM activity and how those cells interact with neighboring immune and stromal cells.
This study integrates three complementary high-resolution approaches, including single-cell RNA sequencing (scRNA-seq), spatial transcriptomics (ST), and consensus machine learning on bulk transcriptomic data, to systematically characterize PTM activity across osteosarcoma cell types and build a prognostic model with direct clinical utility for immunotherapy stratification.
The scRNA-seq component of this study drew from nine publicly available GEO datasets (GSE152048, GSE200529, GSE198896, GSE270231, GSE162454, GSE250015, GSE234187, GSE169396, and GSE217792), collectively yielding 69,342 cells from 37 OS tumor samples and 13 normal bone tissue samples. Data processing used the Seurat R package, with log-normalization applied uniformly across all datasets. Batch effects were mitigated using the Harmony algorithm, and dimensionality reduction was performed via UMAP (uniform manifold approximation and projection). Cell clustering used a resolution parameter of 1.5, producing 49 distinct cell clusters subsequently annotated into seven major cell populations using established marker genes.
PTM activity scoring: PTM activity in each cell was quantified using the AddModuleScore function in Seurat, which calculates an aggregate expression score for a predefined gene set relative to randomly selected background genes. The PTM gene catalog was assembled from KEGG database entries, GSEA gene sets, and published literature, ultimately encompassing 3,401 genes spanning 11 distinct PTM categories. This scoring approach assigns a continuous PTM activity value to each individual cell, enabling classification of tumor cells by their PTM activity level.
Spatial transcriptomics analysis: Spatial transcriptomic data from OS tissue sections were processed with Seurat following log normalization. Spot clustering used the FindClusters and FindNeighbors functions, with initial cluster identification aided by H&E staining and unsupervised analysis. The mistyR package (v1.6.1) evaluated spatial interactions across three distance scales: single spots (intra-region), adjacent regions within a radius of 5 spots (juxta-region), and longer-range neighborhoods within a radius of 15 spots (para-region). The robust cell type decomposition (RCTD) algorithm mapped scRNA-seq-derived cell type annotations onto spatial coordinates, enabling direct comparison of PTM expression patterns across physical tissue sections.
Bulk transcriptomic cohorts: For prognostic model construction, gene expression data were retrieved from TARGET-OS (87 patients, designated as the training cohort) and three GEO datasets (GSE39055, GSE16091, and GSE21257, combined into a validation cohort of 88 patients). The limma and sva R packages were used for batch correction and cohort merging. Pseudotime trajectory analysis was performed using Monocle2 to order osteoblastic cell subpopulations along developmental trajectories.
Across the seven major cell populations identified in the scRNA-seq data (myeloid cells, fibroblasts, osteoblastic cells, T cells, NKT cells, endothelial cells, and B cells), PTM activity was most significantly enriched in osteoblastic cells from tumor tissues relative to normal bone. This preferential enrichment establishes osteoblastic cells as the primary cellular context in which PTM dysregulation drives OS pathology. Based on PTM expression levels, osteoblastic cells were divided into two subgroups: PTMs-high osteoblastic cells (PTMshighos) and PTMs-low osteoblastic cells (PTMslowos).
Pseudotime trajectory analysis: Monocle2-based pseudotime analysis placed osteoblastic cells along a developmental continuum, revealing that molecular markers including LTF, PADI4, and VIM showed elevated expression during the neoplastic (late differentiation) phase. Trajectory analysis identified a transitional pattern in which clusters 7 and 17 (dominant in PTMslowos) progressed toward clusters 6 and 8 (enriched in PTMshighos), indicating that elevated PTM activity correlates with a more terminally differentiated and potentially more aggressive cellular state. This trajectory suggests that PTM accumulation is not a static feature but a dynamic process associated with OS tumor progression.
Pathway enrichment in PTMshighos: GSEA of PTMshighos versus PTMslowos cells revealed that the Notch, Hedgehog, and Wnt signaling pathways, which are well-established drivers of cancer stem cell maintenance and tumor progression, were markedly enriched in the PTMshighos subpopulation. The activation of these three developmental pathways in PTM-high cells suggests that PTM dysregulation may contribute to stemness, resistance to differentiation, and metastatic capacity in OS. The association between PTMshighos cells and a terminally differentiated but highly malignant phenotype is not paradoxical given that in OS biology, terminal differentiation does not necessarily correlate with reduced malignancy.
Reduced immune interaction: PTMshighos cells showed reduced interaction signals with immune cell populations including B cells and T cells in both scRNA-seq and spatial transcriptomics analyses. This finding suggests that elevated PTM activity may contribute to immune exclusion within the tumor microenvironment, potentially suppressing antigen presentation and weakening antitumor immune responses. The implication is that PTM-driven tumor biology is not only cell-intrinsic but also shapes the surrounding immune landscape in ways that may render these tumors less responsive to immunotherapy.
To translate single-cell findings into the physical tissue context, the RCTD (robust cell type decomposition) algorithm was applied to deconvolve spatial transcriptomic data, assigning scRNA-seq-derived cell type identities to individual tissue spots within OS sections. All osteoblastic cells were again classified into PTMshighos and PTMslowos subgroups based on PTM expression, and the spatial distribution of eight distinct cell types was mapped across tumor tissue sections. This analysis confirmed that the patterns observed in scRNA-seq data are not artifacts of tissue dissociation but are spatially organized features of OS tissue architecture.
Spatial interaction analysis using mistyR: The mistyR package evaluated interactions at three spatial scales: intra-regional (within a single spot), juxta-regional (within a 5-spot radius), and para-regional (within a 15-spot radius). Cell-type composition analysis highlighted the intra and para_15 interaction regions as particularly important for modulating cellular interactions and spatial organization within OS tumors. PTMshighos osteoblastic cells showed pronounced co-localization and interaction with fibroblasts at both juxta- and para-regional scales, confirming a spatial relationship between high-PTM tumor cells and cancer-associated fibroblasts (CAFs) that had been suggested by the scRNA-seq interaction analysis.
CAF interactions and immunosuppression: The spatial enrichment of CAF-related signals in the vicinity of PTMshighos cells is biologically significant. Cancer-associated fibroblasts are well-established contributors to an immunosuppressive TME, secreting TGF-beta, CXCL12, and other immunomodulatory factors that exclude cytotoxic T cells and promote tumor immune evasion. The spatial confirmation that PTMshighos cells preferentially neighbor fibroblasts strengthens the model in which PTM dysregulation operates partly through remodeling of the stromal and immune microenvironment rather than exclusively through cell-intrinsic oncogenic signaling.
Concordance between modalities: The convergence of scRNA-seq interaction data and spatial transcriptomics findings provides robust, multi-modal evidence for the PTMshighos-fibroblast axis. The spatial data also confirmed that VIM, one of the ten hub genes in the eventual prognostic model, exhibited broad expression across six distinct cell types in the spatial analysis, consistent with its high expression in five of seven cell types in the scRNA-seq heatmap analysis. This cross-modal validation strengthens confidence in VIM as a functionally relevant prognostic target.
To translate the single-cell PTM findings into a clinically usable prognostic tool, the authors applied a comprehensive machine learning framework to bulk transcriptomic data. Univariate Cox proportional hazards regression identified 97 PTM-related genes significantly associated with OS overall survival across the TARGET-OS and GEO-OS cohorts. These 97 genes were then used as input features for the construction of the consensus machine learning-derived post-translational modification gene signature (CMDPTMS), using the ML.Dev.Prog.Sig function from the Mime R package.
The 117-algorithm ensemble: The Mime package systematically evaluates 117 different algorithmic combinations spanning established machine learning frameworks including StepCox (forward and backward selection), Ridge regression, LASSO, elastic net, random survival forest, gradient boosting machines, and others. Each model is evaluated using K-fold cross-validation in the TARGET-OS training cohort, with performance measured by the mean concordance index (C-index). This exhaustive evaluation avoids the selection bias inherent in testing only one or a few algorithms. The optimal model is selected by maximizing the area under the ROC curve (AUC) using the cal_AUC_ml_res function.
Model selection and performance: Among the top five models by C-index, the StepCox[forward] + Ridge combination achieved the highest mean AUC at both 1-year and 3-year time points, as confirmed by time-dependent ROC curve analysis. The final CMDPTMS model was therefore defined as the StepCox[forward] + Ridge model. It incorporates ten hub genes: FZD8, IBTK, PPP1CC, VPS25, BAG2, TRIM35, UBE3B, USP51, BRD4, and VIM. The model demonstrated a mean C-index of 0.74 across training and validation cohorts, with AUC values at 1-year and 3-year survival endpoints confirming robust predictive performance across external validation datasets.
Benchmarking against existing signatures: In a systematic comparison against published OS prognostic signatures from 2020 to 2025, CMDPTMS demonstrated superior C-index values in both the TARGET-OS and GEO-OS cohorts relative to nearly all reference models. Independent prognostic analysis using a comprehensive nomogram incorporating CMDPTMS risk score alongside clinical variables confirmed the independent prognostic contribution of the risk score, with good calibration curves and decision curve analysis (DCA) indicating net clinical benefit.
Patients were stratified into high-CMDPTMS and low-CMDPTMS groups using the median risk score as a cutoff. Kaplan-Meier survival analysis confirmed that the high-CMDPTMS group had significantly inferior overall survival outcomes. To understand the biological underpinning of this prognostic divergence, the IOBR R package was used to compute enrichment scores for a comprehensive set of published immune-related gene signatures across both patient groups.
Immune cell infiltration differences: High-CMDPTMS patients showed markedly elevated infiltration levels of T cell exhaustion markers, plasma cells, fibroblasts, and myeloid-derived suppressor cells (MDSCs). These are hallmarks of an immunosuppressive tumor microenvironment. By contrast, immunosuppressive molecular markers and immune checkpoint molecules were paradoxically enriched in the low-CMDPTMS group, while immunoexclusion markers including cancer-associated fibroblast (CAF) signatures and TGF-beta family members were predominantly elevated in the high-CMDPTMS group. This pattern is consistent with low-CMDPTMS tumors being "hot" in terms of immune activation potential and high-CMDPTMS tumors being "cold" due to active immune exclusion.
Tumor mutational burden: The low-CMDPTMS group displayed elevated tumor mutational burden (TMB) levels, indicating greater immunogenicity. Higher TMB is generally associated with more neoantigens, broader T cell recognition, and better response to immune checkpoint inhibition. Survival analysis showed that CMDPTMS could complement TMB as a stratification tool: patients with high CMDPTMS combined with high TMB tended to have the most favorable survival outcomes, suggesting that the two biomarkers capture partially independent aspects of immune responsiveness.
Immunotherapy prediction algorithms: Three independent algorithms consistently predicted superior immunotherapy response in the low-CMDPTMS group. TIDE (tumor immune dysfunction and exclusion) analysis predicted lower likelihood of immune dysfunction and exclusion in low-CMDPTMS patients. TIP (tumor immunophenotype) scoring showed marked differences particularly at step 4 (B cell recruiting) and step 6, aligning with the immune cell infiltration findings. Subclass mapping using a melanoma cohort treated with PD-1 blockade confirmed that reduced CMDPTMS levels correlated with more effective PD-1 therapy response. Validation in the IMvigor bladder cancer immunotherapy cohort further confirmed that the low-CMDPTMS group had superior post-treatment survival outcomes, with the divergence becoming apparent after three months of treatment.
Anticancer drug sensitivity was evaluated across CMDPTMS risk groups using IC50 values estimated by the oncoPredict package, drawing on drug response data from the CTRP v2.0 and PRISM Repurposing (19Q4) datasets. The analysis identified differential sensitivity to 37 compounds total: 30 drugs showed sensitivity correlating with low CMDPTMS, while 7 drugs showed sensitivity correlating with high CMDPTMS. This asymmetry suggests that high-CMDPTMS patients, who face worse prognosis and poor immunotherapy response, may nevertheless have a distinct set of therapeutic vulnerabilities that could be exploited.
Identified drug candidates for high-CMDPTMS patients: From the PRISM repurposing dataset, ten compounds showed differential activity in the high-CMDPTMS group: barasertib, BMS-986020, D-4476, GSK2110183, GZD824, masitinib, norfloxacin, ponatinib, SR-27897, and TG100-115. From the CTRP dataset, ten additional compounds were identified including axitinib, AZD8055, PF-4708671, erlotinib, dabrafenib, ribociclib, and ulixertinib. These include kinase inhibitors (barasertib targets Aurora B, ponatinib targets BCR-ABL and FGFR, erlotinib targets EGFR), a CDK4/6 inhibitor (ribociclib), and an mTOR inhibitor (AZD8055), representing multiple mechanistically distinct therapeutic options for patients in the unfavorable risk group.
VIM functional experiments: Among the ten CMDPTMS hub genes, VIM (vimentin) emerged as the most broadly expressed across cell types in both scRNA-seq and spatial transcriptomics analyses. Vimentin is a type III intermediate filament protein classically associated with epithelial-to-mesenchymal transition and cellular motility. To confirm its functional relevance, siRNA-mediated knockdown of VIM was performed in two OS cell lines (MG-63 and U-2 OS) using Lipofectamine 3000 transfection. CCK-8 viability assays measured cell proliferation over 72 hours at 24-hour intervals. Transwell invasion assays assessed migratory capacity by seeding 1 x 10^4 cells in the upper chamber with DMEM plus 15% FBS as a chemoattractant in the lower chamber.
VIM knockdown results: VIM knockdown significantly inhibited both growth and migration in MG-63 and U-2 OS cell lines, validating vimentin as a functional driver of OS cell aggressiveness rather than merely a prognostic correlate. RT-qPCR confirmed elevated VIM mRNA expression in both OS cell lines relative to normal osteoblast cells (hFOB1.19). Fluorescence staining in U-2 OS cells showed cytoplasmic localization of vimentin, consistent with its role as an intermediate filament structural protein. These laboratory findings align with prior reports of elevated VIM and CD63 in plasma exosomes from OS patients compared to healthy individuals.
Sample size limitations: The training cohort consisted of 87 patients from the TARGET-OS database, and the validation cohort comprised 88 patients pooled from three GEO datasets. While the use of multiple independent GEO cohorts strengthens validation relative to a single-center study, these remain relatively small sample sizes for a complex multi-feature machine learning model. The authors explicitly acknowledge that small sample size restricts generalizability, particularly for rare OS molecular subtypes that may be underrepresented in the combined cohort.
OS subtype focus: The study's PTM analysis focused on osteoblastic cells within OS, with findings directly applicable to the osteoblastic subtype. Osteosarcoma includes several histological subtypes (osteoblastic, chondroblastic, fibroblastic, telangiectatic, small cell), and it is not clear whether the CMDPTMS model and the PTMshighos-fibroblast interaction axis generalize across all OS variants. This limits the scope of conclusions that can be drawn about PTM biology in non-osteoblastic OS subtypes, let alone other sarcoma types.
Drug sensitivity prediction limitations: The pharmacological screening identifying 20 candidate compounds relied on in vitro IC50 data from cancer cell line panels (CTRP and PRISM datasets), which are known to have limited predictive value for in vivo drug response in human patients. The tissue context, drug metabolism, tumor heterogeneity, and pharmacokinetics of OS in patients differ fundamentally from cell line experiments. Comprehensive clinical validation of these drug predictions remains imperative and has not been attempted in this study.
Mechanism and clinical translation: The mechanisms by which each of the ten CMDPTMS hub genes contributes to OS prognosis and PTM biology require further mechanistic exploration. The in vitro functional experiments validated only VIM among the ten hub genes. For FZD8, BRD4, UBE3B, and the remaining genes, functional evidence in OS-specific cellular contexts is lacking. The study also notes that practical clinical relevance must be confirmed through extensive multicenter prospective studies before CMDPTMS can be considered for routine clinical application. Integration into clinical decision-making for immunotherapy selection will require prospective trials with uniform treatment protocols and standardized biomarker assessment.