Four major malignant diseases converge on the skeleton: osteosarcoma (OS), Ewing's sarcoma (EW), multiple myeloma (MM), and cancer bone metastases (BM), primarily from prostate cancer. Despite arising from distinct cell lineages with very different clinical behaviors, these four conditions share strikingly similar histological appearances on standard biopsy. They all produce lesions in or adjacent to bone, all cause overlapping radiological signs such as cortical destruction and periosteal reaction, and all can present with localized pain, swelling, and elevated inflammatory markers. This morphological overlap creates a genuine diagnostic hazard in routine clinical practice.
The scale of misdiagnosis: The consequences of getting the diagnosis wrong are severe, since each disease demands a completely different treatment strategy, from chemotherapy combinations for OS and EW, to plasma cell-directed agents for MM, to systemic treatment of the primary tumor site for BM. Published data underscore how frequently errors occur: imaging-based misdiagnosis of osteosarcoma has been reported at rates of 23-43%, particularly in older patients where OS is less common and therefore less clinically suspected. For Ewing's sarcoma, Wurtz et al. found that the average diagnostic delay was approximately 10 months, and a separate study by Widhe et al. documented a misdiagnosis rate of 80.77%. Bone metastasis from an unknown primary tumor is similarly difficult, usually requiring PET-CT plus targeted biopsy confirmation.
The immune microenvironment as a discriminator: Each of these four tumors carries a distinct immune infiltration profile that reflects its underlying biology. Multiple myeloma is fundamentally a B cell malignancy, with bone marrow infiltrated by malignant plasma cells and abundant reactive B lymphocytes. Bone metastatic carcinoma is associated with extensive immunosuppression, with markedly reduced T cell, B cell, and dendritic cell infiltration. Ewing's sarcoma presents with extensive inflammation and high expression of cytokines such as IL-6, often mimicking osteomyelitis radiologically. Osteosarcoma originates from the bone marrow mesenchyme and is dominated by monocyte-macrophage lineage cells. These differences are measurable at the gene expression level, raising the possibility of a computational diagnostic approach based on immune profiling rather than subjective morphology.
This 2021 study from Frontiers in Cell and Developmental Biology tests precisely that hypothesis: can an immune infiltration signature, quantified computationally from gene expression data, be used to construct a machine learning model that accurately differentiates between OS, EW, MM, and BM? The authors also ask a second question: does the immune microenvironment within each tumor type predict patient prognosis and clinical outcomes?
The study assembled 22 publicly available gene expression datasets from two major repositories: the Gene Expression Omnibus (GEO) and the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. These datasets collectively span osteosarcoma, Ewing's sarcoma, multiple myeloma, and prostate cancer bone metastases samples. Careful sample curation was applied: cell line samples, purified tumor cell preparations, duplicated samples, normal tissue samples, and metastatic samples from sites other than bone were all excluded, retaining only primary bone tumors and cancer bone metastasis samples. After cleaning, the total dataset comprised 1,459 samples, split into a discovery cohort of 1,385 samples and an independent external validation cohort of 74 samples.
Batch effect correction: Because the 22 datasets were generated across different laboratories, sequencing platforms, and microarray technologies, cross-platform batch effects were a major methodological concern. Batch effects were removed using the ComBat method implemented in R. Gene expression values were harmonized using Entrez IDs, with the highest mean signal selected when multiple probe sets mapped to the same gene. All data were log2-normalized per dataset before batch effect removal. Within-dataset normalization of ssGSEA scores was deliberately omitted, since published evidence shows that unnormalized ssGSEA scores are more robust for estimating pathway enrichment across independently collected datasets.
ssGSEA immune profiling: Immune infiltration was quantified using the single-sample gene set enrichment analysis (ssGSEA) algorithm implemented via the R package GSVA. This non-parametric method computes an enrichment score (ES score) for each of 28 immune cell types in each individual sample, based on curated immune signature gene sets. Unlike bulk deconvolution methods such as CIBERSORT, ssGSEA does not require a reference matrix and can be applied to any gene expression platform, making it well-suited to integrating datasets across technologies. The 28 cell types profiled include major lymphoid populations (B cells, activated B cells, immature B cells, CD4 T cells, CD8 T cells, effector memory CD8 T cells), NK cell subsets (CD56dim NK cells), myeloid populations (macrophages, monocytes, dendritic cells, MDSCs, neutrophils), and regulatory populations (Treg cells, mast cells, type 2 T helper cells, plasmacytoid dendritic cells, and others).
Random forest classification: The random forest algorithm was chosen to construct the differential diagnostic model. Random forest operates as an ensemble of decision trees, each trained on a random bootstrap sample of the training data using a random subset of features at each split. Final class assignment uses majority voting across all trees. This approach is robust to sample imbalance and overfitting, both of which are concerns when integrating multi-disease datasets of different sizes. The discovery cohort was divided 70/30 into a training set (n = 998) and a testing set (n = 387) via stratified sampling. Model hyperparameters were tuned by varying the number of trees (ntree from 10 to 500) and the number of features sampled per split (mtry from 0 to 15). Optimal performance was reached at ntree = 150 and mtry = 5, with an out-of-bag (OOB) error rate of 2.3%. Further increasing ntree produced no meaningful improvement; increasing mtry reduced accuracy.
Before constructing the diagnostic model, the study characterized the immune landscape of each disease type using unsupervised clustering and pairwise comparisons. Unsupervised hierarchical clustering of the 1,459 samples using Euclidean distances on their 28 ES scores produced a heatmap showing clear groupings consistent with disease type, confirming that immune infiltration profiles carry strong disease-discriminating information. ANOVA testing was used to formally compare ES scores across the four tumor types, and all comparisons were highly statistically significant (p less than 0.0001).
Multiple myeloma: MM showed dramatically elevated B cell enrichment scores, consistent with its biology as a malignant plasma cell disorder arising from the B cell lineage. Activated B cell and immature B cell ES scores were so distinctively elevated in MM that they alone could reliably separate MM from the other three tumor types in pairwise scatter plots. The specificity for MM using activated B cells alone reached 1.0 in both the test and validation datasets, confirming B cell abundance as a near-perfect positive marker for MM in this context.
Osteosarcoma: OS showed the highest enrichment of macrophages and CD56dim NK cells among the four tumor types. Treg cell infiltration was also more prominent in OS than in EW or BM. The correlation matrices within each disease revealed that OS and EW both have generally positive correlations between different immune cell types, meaning that when one immune population is high, others tend to be high as well. This suggests coordinated immune activation in both sarcoma subtypes, in contrast to MM (where immune correlations were weak or absent) and BM (which showed a complex, partly negative immune cell correlation network consistent with an immunosuppressive environment).
Bone metastases: BM from prostate cancer showed consistently reduced infiltration of T cells, B cells, and plasmacytoid dendritic cells compared to the primary bone tumors, consistent with the established concept of an immunosuppressive bone metastatic niche. This reduced infiltration pattern was visible in both the heatmap clustering and the pair plot analysis. Ewing's sarcoma: EW did not show a single dominant immune cell subpopulation that could serve as a standalone diagnostic marker, distinguishing it from both MM (B cells) and OS (macrophages). Instead, EW showed broadly elevated infiltration across multiple immune cell types, consistent with its clinical presentation as a highly inflammatory tumor. This finding implies that single-variable EW diagnosis is not feasible from the immune microenvironment alone, and that a multi-variable model such as random forest is necessary.
The random forest model was evaluated on three separate datasets: the training set, the internal test set, and the independent external validation set. Variable importance analysis identified monocytes, CD56dim NK cells, and activated B cells as the three most important features for classification accuracy, while activated B cells, immature B cells, and CD56dim NK cells contributed most to reducing sample variance (Gini impurity) across the decision trees. MDS (multidimensional scaling) plots of the 28 ES scores confirmed that the four tumor types occupy clearly separable regions of feature space, providing visual confirmation that the classification problem was learnable.
Internal test set performance (n = 387): The overall diagnostic accuracy of the model on the internal test set was 97.42% (95% CI: 95.3-98.75%, p less than 2.2 x 10 to the -9). Per-class sensitivity values were 0.828 for BM, 0.954 for EW, 1.000 for MM, and 0.992 for OS. Per-class specificity values were 1.000 for BM, 0.997 for EW, 0.996 for MM, and 0.970 for OS. F1 scores, which balance precision and recall, were 0.906 for BM, 0.976 for EW, 1.000 for MM, and 0.964 for OS on the test dataset. BM had the lowest sensitivity (0.828), reflecting the greater challenge of distinguishing immunosuppressed bone metastasis samples from the other groups when infiltration levels are generally low.
External validation performance (n = 74): On the completely independent external validation cohort, overall accuracy improved to 98.65% (95% CI: 92.7-99.97%, p less than 2.2 x 10 to the -16). Sensitivity values reached 1.000 for BM, 1.000 for EW, 0.941 for MM, and 1.000 for OS. Specificity reached 1.000 for BM, 0.984 for EW, 1.000 for MM, and 1.000 for OS. F1 scores were 1.000 for BM, 0.952 for EW, 0.970 for MM, and 1.000 for OS. The fact that external validation performance matched or exceeded internal test performance is a meaningful finding, suggesting that the model generalizes rather than simply overfitting to the training distribution.
Comparison with single-variable models: ROC curves compared the full random forest model against single ES score diagnostic models built from each of the five most important variables. For each of the four tumor classes, the random forest model achieved substantially better AUC than any single immune cell ES score alone. This confirms that no single immune marker is sufficient for reliable differential diagnosis and that the ensemble model captures complementary information from multiple immune features simultaneously.
Having established the immune infiltration landscape, the study then examined whether ES scores predict patient prognosis within each disease. Four OS datasets (TARGET OS, GSE21257, GSE16091, GSE29055) and two EW datasets (GSE17618, GSE63157) contained clinical overall survival data with sufficient follow-up for prognostic analysis. Cox univariate regression was performed for each immune cell type across each dataset separately, followed by Kaplan-Meier (KM) survival analysis using the per-dataset median ES score as the cutoff for high versus low infiltration.
Macrophage infiltration in osteosarcoma: Macrophages showed a consistent positive trend in all four OS datasets, with high macrophage ES score associated with better patient survival. However, when evaluated per-dataset using each dataset's own median as the cutoff, the prognostic effect was statistically significant in only some datasets, reflecting heterogeneity in cohort composition and follow-up time. To overcome this instability, the authors pooled all four OS datasets and applied the global median macrophage ES score (computed across all 1,459 samples of all four diseases) as a fixed, unbiased cutoff. In the combined OS analysis, patients with low macrophage infiltration had significantly shorter overall survival (p = 0.009), supporting macrophages as a favorable prognostic marker in OS.
B cell infiltration in osteosarcoma: Activated B cells demonstrated statistically significant prognostic value in the GSE39055 dataset when analyzed per-dataset, and remained significant in the combined OS dataset analysis (p = 0.032), with higher B cell infiltration associated with better survival. This is a notable finding because B cells have historically received very little attention in osteosarcoma research, which has focused predominantly on T cell and macrophage responses. Recent oncology literature has identified B cell markers and tertiary lymphoid structures (organized B cell-rich lymphoid aggregates within tumors) as strong predictors of immunotherapy response in multiple cancer types, including soft tissue sarcoma. The present study extends this observation to osteosarcoma at the molecular immune microenvironment level.
Clinical associations with immune infiltration in OS: Beyond overall survival, immune ES scores were correlated with specific OS clinical characteristics using t-tests. CD4 T cell infiltration differed significantly between male and female patients, though other immune cell types did not show sex-based differences. Using age 18 as a cutoff, younger OS patients showed generally lower immune infiltration in 26 of 28 cell types, suggesting that immune microenvironment maturity may differ by developmental stage. Patients with tumor metastasis had significantly lower ES scores for activated B cells, immature B cells, activated dendritic cells, effector memory CD8 T cells, MDSCs, NK T cells, neutrophils, and type 2 T helper cells, while plasmacytoid dendritic cell ES score was elevated, pointing toward specific immune escape mechanisms associated with OS metastasis.
The prognostic relationship between immune infiltration and outcomes in Ewing's sarcoma was strikingly inverse compared to osteosarcoma. In the combined EW dataset analysis (two datasets pooled, global median as cutoff), 17 of the 28 immune cell types showed a negative correlation with patient survival, meaning that higher immune infiltration was associated with shorter overall survival. This is the opposite of the pattern seen in OS and of the general expectation that immune infiltration is favorable in solid tumors. The CD56dim NK cell ES score had the most consistent negative prognostic association across both EW datasets.
Inflammatory EW biology: EW frequently presents with fever, elevated ESR, and radiological changes that mimic osteomyelitis, and its gene expression profile includes high levels of pro-inflammatory cytokines such as IL-6. The clinical literature on EW immunobiology has noted that increased infiltration of M2-like macrophages, CD8 T cells, NK cells, mast cells, B cells, and type 2 T helper cells in the EW microenvironment correlates with poor patient prognosis. The present study's finding that 17 of 28 immune cell types negatively predict EW survival is consistent with this body of literature and suggests that the extensive inflammatory tumor microenvironment in EW may actively promote rather than restrain tumor progression.
Mechanistic implications: The paradoxical finding that immune infiltration is favorable in OS but unfavorable in EW likely reflects fundamental differences in how these two sarcomas interact with immune cells. In OS, macrophages appear to play a protective role (supported by clinical evidence that mifamurtide, a macrophage activator, improves OS disease-free survival, and that high CD163-positive macrophage infiltration correlates with better OS outcomes). In EW, the pervasive inflammatory environment, potentially driven by the EWS-FLI1 fusion oncogene, may drive tumor growth and immune evasion through different mechanisms, including cytokine-mediated suppression of cytotoxic immune responses. These opposing patterns highlight that the immune microenvironment cannot be interpreted in a tumor-type-agnostic manner; the same immune cell population may have opposing effects in different sarcoma subtypes.
Chemotherapy responsiveness: In the osteosarcoma dataset where Huvos grading (a pathological measure of chemotherapy-induced tumor necrosis) was available, patients with high Huvos grade (greater than or equal to 90% necrosis, indicating excellent chemotherapy response) had significantly more CD8 T cell, mast cell, and type 2 T helper cell infiltration compared to low-grade responders. This links pre-treatment immune infiltration to chemotherapy efficacy, raising the possibility that immune microenvironment profiling could help predict which OS patients will respond well to neoadjuvant chemotherapy.
Retrospective public data only: The entire study is based on publicly available retrospective gene expression datasets collected for purposes other than this analysis. No prospective or clinically controlled data collection was performed. The 22 datasets span a wide range of collection periods, geographic origins, patient selection criteria, and gene expression platforms. While the ComBat batch correction method is widely used and validated, it cannot perfectly remove all sources of technical and biological variation introduced by cross-platform integration. Some residual batch effects may persist and could artificially inflate or deflate immune ES scores for specific cell types.
Limited sample size for bone metastasis: The BM cohort comprised prostate cancer bone metastases only, with 109 BM samples in the training/testing cohort. This limits the generalizability of the model to other primary tumors that commonly metastasize to bone, including breast cancer, lung cancer, and renal cell carcinoma. The prognostic analyses for BM and MM were not performed because suitable clinical outcome datasets were not available in the public repositories for these disease categories. The external validation cohort of 74 samples is also small, and while performance on this set was strong, larger independent cohorts would provide more reliable confidence intervals.
No histological validation: The diagnostic model operates entirely on gene expression-derived immune infiltration scores. It was not validated against expert hematopathologist review, IHC panels, or cytogenetics in a head-to-head comparison, so it is not possible to directly assess whether it would outperform or complement current standard-of-care diagnostic tools. For clinical deployment, a prospective comparison study with histopathological gold-standard diagnosis as the reference would be required.
Computational nature of immune profiling: ssGSEA immune infiltration scores are derived from bulk gene expression data, meaning they represent average infiltration across all cells captured in a tissue sample. Single-cell RNA sequencing would provide far higher resolution on the spatial distribution, activation states, and functional phenotypes of immune cells within these tumors. It is possible that the prognostic associations identified reflect the abundance of immune cell types without capturing whether those cells are functionally active, exhausted, or phenotypically distinct in ways relevant to patient outcomes.
Expanding to a diagnostic complement: The authors note that ssGSEA requires only a single sample's gene expression data as input, without needing a full reference matrix. This means the model could in principle be applied to a biopsy from a single patient using microarray or RNA-seq data, functioning as a computer-aided diagnostic tool alongside conventional pathology. In this context, the model would not replace histopathology but could serve as a quantitative second-opinion tool, particularly in ambiguous cases where morphology alone does not allow confident classification. Prospective validation studies embedding the model in a clinical workflow, with blinded comparison against expert panel diagnoses, would be the necessary next step.
Expanding the BM panel: The current model was trained on prostate cancer bone metastases. Future work should incorporate bone metastasis samples from breast cancer, lung cancer, renal cell carcinoma, and other common primaries, since each may carry a distinct immune microenvironment that differs from prostate BM. Building a broader BM training set would improve the model's clinical utility and its sensitivity for BM arising from different primary sites. This would require curating additional public datasets or prospective collection, with careful attention to batch correction across newly added platforms.
Therapeutic targeting of the immune microenvironment: The prognostic data on macrophage and B cell infiltration in OS, and the broadly negative immune prognostic signal in EW, suggest specific therapeutic hypotheses. For OS, stimulating macrophage and B cell infiltration (for example with macrophage activators such as mifamurtide, or with strategies promoting tertiary lymphoid structure formation) could represent a novel adjuvant approach. For EW, the extensive inflammatory microenvironment associated with poor prognosis suggests that targeting specific pro-tumorigenic immune populations (such as M2-like macrophages or suppressive NK cell subsets) rather than broadly promoting immune infiltration may be needed. The study provides molecular data to prioritize these hypotheses.
Integration with single-cell and spatial technologies: Single-cell RNA sequencing and spatial transcriptomics would allow the immune infiltration observations made at bulk resolution to be refined to the level of individual cell types, activation states, and tissue compartments. Particularly for osteosarcoma, identifying whether the favorable macrophage and B cell signals originate from specific spatial niches (such as tertiary lymphoid structures or the tumor-bone interface) would deepen the biological understanding and potentially identify more precise therapeutic targets. The ssGSEA-based approach provides a scalable first layer of immune characterization that can guide subsequent single-cell studies toward the most biologically informative populations.